Structure and dynamics in drug discovery

Introduction

Developing new drugs is very expensive and time-consuming^1,2. Generally speaking, the drug discovery and development process usually takes around 10–14 years and more than 1 billion dollars³. In the 1980s, the use of computers was extended from data handling to a more prominent role in drug discovery⁴, and since then, computational approaches in the drug discovery, design, and development process have been adopted rapidly^{2,5,6,7,8,9,10}. It has been estimated that by the use of computer-aided drug discovery (CADD) approaches, the cost of drug discovery and development can be reduced up to 50%¹¹.

CADD is a specialized discipline that uses computational methods to simulate drug-receptor interactions to determine if a given molecule will bind to a target, and if so, what its affinity would be¹². There are mainly two types of CADD techniques: ligand-based drug design (LBDD) and structure-based drug design (SBDD). SBDD can only be used when the three-dimensional structures of target proteins are known, while LBDD design is employed in situations that the structures are unavailable. In recent years, there have been simultaneous advancements in structural biology such as cryo-EM, and computational protein structure prediction as with AlphaFold¹³, allowing for access to many more high-resolution 3D structures of novel drug-receptor complexes^14,15,16,17. Because of the convergence of these breakthroughs, SBDD, as a structurally enabled computational method (Fig. 1), has become a driving force for the discovery of novel small molecule therapeutics.

Structure and dynamics in drug discovery — **Fig. 1: Schematic of structure-based drug discovery.**

Once a well-defined drug target structure is obtained, virtual screening of libraries of compounds is conducted, through molecular docking software. By scoring and ranking a collection of ligands, potential drug candidates are identified. Thus, the rapid expansion of drug-like chemical space, easily accessible for hit and lead discovery, is another key factor driving the advancement of CADD. Just a few years ago, the compound libraries were still limited to several million on-shelf compounds from vendors and in-house screening libraries in pharma. Now, screening can be done with ultra-large virtual libraries and chemical spaces of drug-like compounds, which can be readily made on-demand, rapidly growing beyond billions of compounds¹⁸.

However, there are certain limitations of SBDD despite its many successes. One of the biggest remaining challenges in SBDD is target flexibility. Proteins and ligand molecules possess high flexibility in solution and undergo frequent conformational changes. However, most molecular docking tools allow for high flexibility of the ligand, but the protein is kept fixed or provided with only limited flexibility to the residues present within or near the active site. It is very difficult to provide complete molecular flexibility to the protein as this increases the space and time complexity of the computation dramatically. The second challenge is the difficulty of exploring cryptic pockets, which are not shown in the original structure but will be revealed by protein conformation changes. Those pockets often relate to allosteric regulations, which would offer extra opportunities targeting beyond the primary endogenous binding site of the receptor.

One solution to the above challenges is molecular dynamics simulation, also referred as MD simulation¹⁹. Nowadays, MD simulation has become one of the most influential methods for modeling conformational changes within a ligand-target complex upon binding a small molecule^20,21. Nevertheless, it is commonly believed that normal MD is unable to cross substantial energy barriers within a simulation’s lifespan, preventing it from efficiently traversing the energy landscape of a protein in complex with ligands. By adding a boost potential to smooth the system potential energy surface, accelerated molecular dynamics^22,23 (aMD) methods were developed for helping with this issue. The boost potential decreases the energy barriers and therefore accelerates transitions between the different low-energy states²⁴. With this, aMD is able to sample distinct biomolecular conformations and help with the receptor flexibility and cryptic pockets problems.

In a word, the discovery of new therapeutic drugs has been greatly speeded by computer-aided methods based on molecular structure and dynamics. This perspective provides a brief appraisal of the origin, role, and prospects of these methods, starting from early work based on x-ray crystallographic structures of proteins, through the use of molecular dynamics simulations, to early work drawing upon simulation and machine learning.

Structure-based drug discovery

Increase in available structures

Perhaps the earliest use of protein structures for drug discovery was that for the development of the important inhibitors of angiotensin-converting enzyme (ACE), captopril, and enalopril. The design of these drugs, used to treat high blood pressure and other conditions, benefitted from modeling based on the crystallographic structure of carboxypeptidase A, which has a similar active site featuring a catalytically important zinc ion^25,26.

This field has grown rapidly in the ensuing years, due in part to the elucidation of the structures of thousands of proteins, nucleic acids, and other potential drug targets²⁷. In recent years, due to leaps in structural biology, including automation in crystallography¹⁵, microcrystallography¹⁴, and cryo-electron microscopy technology^16,17,28, the 3D structures for many clinically important targets have been revealed, often in a state relevant to its biological function (Fig. 2). Especially impressive has been the recent structural revolution for G protein-coupled receptors (GPCRs)²⁹, ion channels^30,31 and other membrane proteins that mediate the action of more than half of drugs³², providing excellent targets for ligand screening and lead optimization.

**Fig. 2: The increasing of available structures.**

The number of target structures has increased significantly with the arrival of machine learning tools such as AlphaFold, which reliably predict the atomic structure of proteins for which experimental structures may not be available¹³. Since its launch in 2021, the AlphaFold Protein Structure Database has had released over 214 million unique protein structures³³, compared to around 200,000 PDB structures corresponding to approximately 60,000 unique protein sequences (Fig. 2). This new data set almost covers the complete UniProt database. Additionally, AlphaFold models can cover the entire length of protein sequences compared to the fragmented, often short coverage of PDB entries. Clearly, researchers and pharmaceutical companies can try structure-based drug discovery techniques using these models, presenting unprecedented opportunities for targets without a prior experimental structure.

Improving computational methods

The early years of structure-based drug discovery typically featured computational studies of the binding of small molecules to active sites of enzymes, using structure-based virtual screening or molecular docking, with the test molecules chosen from libraries of compounds or analogs of known binders³⁴. The quality of binding of each test molecule was typically determined using model potential energy functions chosen to balance speed and accuracy³⁴.

Docking molecules of a virtual drug-like compound library into a target receptor structure and predicting its binding score is a major step in a structure-based drug discovery campaign, which plays a key role in any successful application^7,29,35. The predicted candidate ligand sets, produced by such virtual screening, usually show useful hit rates, about 10%-40% in experimental testing³⁶. Some novel hits may also exhibit noteworthy potencies, in the 0.1–10-μM range, for different types of targets³⁶.

Special attention has been devoted to ligand scoring functions, which are supposed to reliably select top binders and to rule out false-positive predictions. This is especially important with the growth of library size. For example, a one-in-a-million rate of false positives in a billion-compound library would result in a thousand false hits, which obviously complicates the selection of hit candidates. Another major challenge is the computation cost. With increasing library sizes, the computational time of docking itself is the main bottleneck in virtual screening processes. Nowadays, screenings on ultra-large virtual libraries that include billions of drug-like compounds are feasible, thanks to the recent availability of cloud computing and graphics processing unit (GPU) computing resources³⁶.

Expansion of accessible chemical space

A successful structure-based drug discovery screening campaign depends critically on diverse ligand libraries that cover a large part of the chemical space that might be of interest. There are at least two reasons why the compound library should be large and diverse. Firstly, such a library will obviously increase the chance of identifying potential hits in the virtual screening process³⁷. This has been demonstrated in many ultra-large virtual screening campaigns, for example in refs. ^{38,39,40,41,42}. Secondly, the diversity of drug candidates will be improved, offering more opportunities for further optimization and modification. A good library can expand the chemical diversity, novelty, and patentability of the hits⁴³. The hit analogs in the same library can help build a meaningful structure-activity relationship, which further facilitates downstream optimization steps.

Several approaches have been developed recently to boost the size and diversity of screening libraries, among which the virtual on-demand libraries are most worth mentioning. In 2017, the readily accessible (REAL) database by Enamine¹⁸ was established and became the first commercially available on-demand library. The REAL library uses carefully selected in-stock building blocks and optimized parallel synthesis protocols, which makes it a fast, reliable source of compounds³⁹. Through years of development, the fully enumerated REAL database has grown from approximately 170 million compounds in 2017 to more than 6.7 billion compounds in 2024. The successful application of the REAL database has been recently documented in several virtual drug screening campaigns⁴⁴, some of which showed exquisite performance with nanomolar and even sub-nanomolar affinities. Some other similar ultra-large virtual libraries are also available, for example, synthetically accessible virtual inventory (SAVI)⁴⁵ developed by the US National Institutes of Health.

Dynamics-based drug discovery

The relaxed complex method

Structure-based screening and docking methods are often able to sample accessible conformations of the small molecules being tested but, apart from modest relaxation, they are typically limited in the extent to which they can sample different conformations of the target molecule. Related to this, the docking methods are limited in their ability to estimate free energies of binding.

The introduction of molecular dynamics (MD) simulations for the study of proteins opened the way to improved methods for drug discovery^46,47. MD simulations allow sampling not only of the conformations of the ligands, but also those of the target molecule. In the target molecule, pre-existing pockets for potential binding vary somewhat in size and shape during its normal dynamics. Importantly, cryptic pockets may also appear, providing new binding sites that might be accessed in docking and modeling. A systematic approach to representing this variation in potential binding sites is the Relaxed Complex Method (Fig. 3), in which representative target conformations, often including novel, cryptic binding sites, are selected from MD simulations for use in docking studies^48,49.

**Fig. 3: The workflow of dynamic-based drug discovery.**

An early example of the effective use of MD in drug discovery is the development of the first FDA-approved inhibitor of HIV integrase. MD simulations⁵⁰ starting with x-ray crystallographic structures of the core domain of this integrase, determined by the Davies group⁵¹, provided early indications of significant flexibility in the active site region. The Davies group subsequently reported the first crystal structure of an inhibitor-bound form of the protein⁵². The Relaxed Complex Method was subsequently applied to explore the binding of this inhibitor to different conformations of the protein generated by MD⁵³. Surprisingly, the results showed that the binding pocket seen in the crystallographic structure expanded to yield a trench in the protein surface that could accommodate two different binding poses of the inhibitor, rotated by about 180 degrees and translationally separated⁵³. This discovery was noted by researchers at Merck & Co⁵⁴, who pursued the two-pose concept in developing the first FDA-approved inhibitor of the integrase for clinical use⁵⁵.

The enhanced sampling methods

Classic MD often shows a debilitating timescale problem, since its integration timestep is restricted to a few femtoseconds. As such, a large number of enhanced sampling methods have been proposed over the decades, to accelerate the dynamics of the system and access much longer timescales^56,57,58. These methods can be broadly classified into two categories.

The first category includes methods based on collective variables. These methods, such as umbrella sampling⁵⁹, metadynamics⁵⁷, and the weighted ensemble method⁶⁰, utilize collective variables to represent the degrees of freedom of interest, thereby reducing the dimensionality of the free energy surface and enabling more comprehensive sampling of these regions. However, selecting effective collective variables can be challenging and often requires prior knowledge of energy basins and their potential exiting pathways. As a result, the applicability of such methods is somewhat constrained. For instance, weighted ensemble methods are particularly effective for sampling transitions between well-defined end states.

The second category encompasses methods that do not rely on collective variables and comprises a variety of distinct approaches. Replica exchange dynamics, for example, employs multiple copies of the same simulation at different temperatures to enhance the probability of observing rare events⁶¹. The mixed Monte Carlo–MD method^56,62 combines Monte Carlo simulations with molecular dynamics to explore extended timescales. Additionally, accelerated MD methods have gained increasing attention in recent drug discovery research, for their advantages in simulating target conformations and identifying cryptic binding sites.

Accelerated MD methods^22,23 are enhanced sampling approaches that work by adding a nonnegative boost potential to smooth the system potential energy surface. No collective variables are employed, so that the results can be directly reweighted to yield the probability of occurrence of conformations and other thermodynamic properties³⁸. This allows the accelerated MD methods to be applied in most situations. A more recent Gaussian Accelerated Molecular Dynamics^63,64 (GaMD) method has been developed in which the boost potential follows a Gaussian distribution and allows for simple and accurate reweighting of the simulations.

Employing these accelerated MD simulation techniques has at least two benefits: better accounting for the target flexibility and increasing the likelihood of finding cryptic binding sites. In Miao, et al.⁶⁵, the M2 muscarinic acetylcholine receptor (mAChR) was used as a GPCR model, and accelerated MD simulations were used to account for the receptor flexibility. Through iterative molecular docking and experimental testing, the authors successfully identified positive and negative allosteric modulators of M2 mAChR, with unprecedented chemical diversity and, remarkably, demonstrated in vivo selectivity of the targeted receptor in a family of similar receptors. In Seitz, etc.⁶⁶, through GaMD, computational screening and in vitro activity testing, the authors reported novel scaffolds for inhibiting Mycobacterium tuberculosis cytochrome bd oxidase, again validated experimentally.

Free energy calculations

An important advantage in the use of MD in drug discovery is that it opens the way to the computation of free energies of binding. The components of free energy of a system, such as enthalpy and entropy, reflect averages over its allowed structures. Free energy studies are a natural next step after library screening and docking studies, providing a way to rank promising lead compounds. Computation of relative free energies of binding, which allows for estimation of the relative affinities of congeners or other pairs of molecules, can help to guide screening or synthetic work. Such calculations can be done using the computational alchemy formalism, first demonstrated for proteins in the binding of different ligands to the enzyme trypsin and a mutant form of the enzyme^67,68. Standard free energies of binding of individual ligands is also possible, as demonstrated in 1988⁶⁹ and then in a more complete form in 1997⁷⁰.

Application in hit optimization

Iterative design-make-test cycles are the essential way to turn an initial hit molecule into a lead drug candidate, which requires properties that are aligned with relevant clinical and pharmaceutical demands. Eventually, the developed drug will need to balance a group of properties like drug efficacy, efficiency, toxicity, human tolerance, etc. The early hit-to-lead stage is focused primarily on binding affinity because high-affinity compounds are likely to be active at lower doses and to have longer residence times in the relevant receptors. Much progress has been made in the past decades for rapidly optimizing binding affinity, most notably by utilizing free energy perturbation (FEP) calculations^{71,72,73,74,75,76,77,78,79,80,81}. Previously, it has been challenging for free energy calculations to achieve the accuracy, efficiency, and reliability required for application to the hit optimization process. Fortunately, recent advances in methodologies and the development of specialized parallel computing devices (GPUs), have enabled free energy calculations to attain a level that is sufficient to drive better synthesis decisions⁸². Such free energy calculations, combined with molecular dynamic simulations, can be used to examine molecular conformations, molecular motion, intermolecular interactions, etc. With this level of accuracy and insight, medicinal chemists can test their ideas based on the FEP calculations before actual synthesis and reduce the percentage of compounds synthesized that do not meet potency requirements.

Application in initial screening

Due to its high computational costs (compared with docking), free energy calculations are not usually applied in the initial screening phase. However, recently both industrial and academic groups are embracing binding free energy calculations as a screening tool, as featured by a few prospective applications. The advantage of free energy calculation is that it is more accurate than the conventional docking process, yet more rapid than what can be done with experiments, which allows drug discovery teams to explore a larger part of chemical space. There has been a growing literature, as we show below, demonstrating the value of free-energy simulations in early drug discovery applications.

In 2016, Merck KGaA started a large initiative to prospectively assess the prediction accuracy of relative binding free energy calculations. In this large-scale assessment concerning 12 targets and 23 chemical series which are from their active drug discovery projects, more than 35,000 free energy perturbation calculations were done, over 400 novel molecules with predicted high binding affinities were identified and subsequently synthesized and tested⁸³. Another example came from Bayer, where a team used a computationally empowered workflow to identify novel covalent allosteric binders for the KRAS G12C isoform⁸⁴. In one of the largest free-energy simulation studies on a GPCR, a collection of 3,4- dihydropyrimidine-2(1H)-ones were enumerated and computationally screened with free-energy simulations against the A2B adenosine receptor⁸⁵. In another study of a viral protein, free-energy simulations were performed on the crystal structure of norovirus RNA-dependent RNA polymerase in a complex with several known binders to determine binding free energies of these molecules relative to the natural nucleotide substrates, through which a virtual nucleotide library containing 121 molecules was screened and two novel molecules were successfully identified with in vitro activity⁸⁶.

Machine learning

Machine learning is a promising new tool for drug discovery. Perhaps the major contribution of machine learning in structure-based drug discovery to date has been through the dramatic expansion of available target structures. In the context of structure and dynamics methods for ligand-receptor recognition, AlphaFold, and related methods are already providing input for protein modeling studies utilizing MD⁸⁷.

AlphaFold

One of the most dramatic developments in recent years in structural biology is the appearance of AlphaFold¹³. To the date of this review, DeepMind released the updated version, AlphaFold 3⁸⁸. Accurate target structures are critical to the mission of structure-based drug discovery, as the starting point of any campaign. Since the introduction of the first AlphaFold model¹³, a revolution in modeling the structure of proteins has started. Enormous progress has been achieved in protein structure prediction with the development of AlphaFold 1, and the field has grown tremendously with a number of later methods that build on the ideas and techniques of AlphaFold 2^{33,89,90,91,92,93,94,95,96}. Clearly, the emergence of AlphaFold models opens unprecedented opportunities for drug discovery with targets that were not available before.

AlphaFold 3 is capable of predicting the joint structure of complexes including proteins, nucleic acids, small molecules, ions, and modified residues. This is promising for structure-based drug discovery, as it is commonly believed that a complex structure with both target and a known binder is a better starting point than a single apo structure. The authors of AlphaFold 3 claim in their publication that the performance of their models on protein-ligand systems is better than classical docking tools, such as Vina^97,98, and greatly outperforms all other blind dockings like RoseTTAFold All-Atom. Their evaluation was done on their PoseBusters benchmark set, which is composed of 428 protein-ligand structures that were not included in their training. The reported accuracy, as the percentage of protein-ligand pairs with pocket-aligned ligand root mean squared deviation of less than 2 Å, is over 90% for their high-confidence group.

AlphaFold and other Machine learning models may be limited in the near future in their ability to predict cryptic pockets, given limitations in training data. It may be possible to create training data by enhanced sampling MD and use this data to train machine learning models to produce additional target structures. Related research is in Lyu, etc.⁹⁹, where the authors followed up on several hundred computational hits and found that there was little to no overlap for the same receptor when starting with the AlphaFold 2 model versus the experimental structure. This indicated that AlphaFold models are already showing some potential on modeling different conformations.

Machine learning docking scoring

One of the attractive directions for applying machine-learning techniques in structure-based drug discovery is the generation of machine-learned docking score functions, by extracting information from both known ligand activities and corresponding protein-ligand 3D structures. The time seems ripe for such models as the databases needed have been established, for example, the PDBbind database¹⁰⁰. Thus, there have already been some attempts in this direction, with various approaches to represent the data and different network architectures, including spatial graph-convolutional models^101,102, 3D deep convolutional neural networks^103,104 or their combinations¹⁰⁵. However, most of these models have suffered severe overtraining problems¹⁰⁶. One possible explanation for this phenomenon is that the PDBbind database does not have an appropriate number of negative cases, i.e., ligands with non-optimal binding interfaces to enrich the training set. One solution to this could be to train the model with the results of physics-based docking, in addition to experimental data. A recent example of such a model is RTCNN¹⁰⁷, although its practical utility remains to be demonstrated.

Since the scoring functions, whether physics-based or machine-learned, are primarily designed or trained to effectively separate potential binders from non-binders, it is probably the best idea to combine both of them and wield their maximum advantages in different situations. If the physics-based and machine-learned scoring functions are relatively independent and emphasize different interactions between the target and ligands, their combination can perhaps reduce the false-positive rates and improve the quality of the hits. For example, in Durrant and McCammon¹⁰⁸, they describe a fast and accurate neural-network-based scoring function that can be used to re-score the docked poses of candidate ligands. Another such example is in the latest 3DR Grand Challenge 4 results for ligand IC50 predictions¹⁰⁹, in which the top methods that used a combination of both physics-based and machine-learned scoring outperformed others. Further, for more accurate potency predictions, the smaller focused library of candidate binders selected by the initial machine learning methods can be sent to more elaborate physics-based analysis, such as free energy perturbation methods^{83,110,111,112}. There is no doubt, with its striking capabilities, machine learning will play an important role in virtual screening campaigns, as well as in lead optimization stages.

Other machine-learning models

Early machine learning algorithms such as support vector machine, and random forest have been used to predict ligand activities and properties, even before deep learning/neural networks came into wider use. For example, quantitative structure-property relationship (QSPR) models which were used to predict pharmacodynamic and pharmacokinetic properties, such as solubility, lipophilicity, etc., have been trained based on large and broad experimental datasets^113,114,115. Machine learning techniques are also implemented for many quantitative structure-activity relationship (QSAR) algorithms¹¹⁶, where the training set and the resulting models will focus on a specific target and one type of chemical scaffold. These machine learning models can greatly help on guiding the hit-to-lead and lead optimization process. Other machine learning models aimed at drug repurposing have also been suggested, based on extensive ligand-target binding datasets, chemical similarity clustering, etc.^117,118.

For QSPR, large public and private databases have been accumulated, with various properties such as solubility, lipophilicity or in vitro proxies for oral bioavailability and brain permeability experimentally measured for many thousands of diverse compounds, allowing prediction of these properties in a broad range of new compounds. However, things are not quite similar for the situation of QSAR models, because, for different target classes, the data availability vary greatly. Probably, the most abundant QSAR databases are those for the kinase superfamily and aminergic GPCRs, as they are popular as drug targets. Recently, an unbiased benchmark of the best machine learning QSAR models was given by a IDG-DREAM Drug-Kinase Binding Prediction Challenge, with the best models achieving a root-mean-square error of 0.95 for the predicted versus experimental pK_d values¹¹⁹. Such accuracy may be acceptable in screenings for the initial hits for less explored kinases and guiding lead optimization. However, as noted before, the kinase family is unique as it is the largest target family with more than 500 structures. The performance and generalizability of such machine learning models for other less rich target families remain unclear.

Future directions

In the last decades, with all the successes and challenges, the transformation from computer-aided drug discovery to computer-driven drug discovery is emerging. The rapid advancements in structural biology along with the bloom of computational protein-structure prediction, are allowing for access to many more high-resolution 3D structures of novel drug-receptor complexes. Dynamic-based drug discovery strategies, coupled with advanced sampling methods, are offering new possibilities for drug design. During the early hit identification stage, the ultra-scale virtual screening approaches, both structure-based and AI-based, are becoming ready in providing fast and cost-effective start points into drug discovery campaigns. The rapid expansion of accessible chemical space, together with the dramatic increases in computational power such as GPUs and cloud computing, are resulting in the ability to virtually screen multi-billions of drug-like chemical space. At the hit-to-lead stage, the more elaborate potency prediction tools such as free energy perturbation and AI-based QSAR often guide rational optimization of ligand potency. In the end, data-driven computational tools are used for multi-parameter optimization, including solubility, permeability, and pharmacokinetic properties, to identify the final drug candidate. It is evident that these breakthroughs are converging towards a new era of computer-driven drug discovery (Fig. 4).

**Fig. 4: Paradigm of future computer-driven drug discovery.**