Insights into the dynamics of HIV-1 protease: a kinetic network model constructed from atomistic simulations

Contingency and Entrenchment of Drug-Resistance Mutations in HIV Viral Proteins

The ability of HIV-1 to rapidly mutate leads to antiretroviral therapy (ART) failure among infected patients. Drug-resistance mutations (DRMs), which cause a fitness penalty to intrinsic viral fitness, are compensated by accessory mutations with favorable epistatic interactions which cause an evolutionary trapping effect, but the kinetics of this overall process has not been well characterized. Here, using a Potts Hamiltonian model describing epistasis combined with kinetic Monte Carlo simulations of evolutionary trajectories, we explore how epistasis modulates the evolutionary dynamics of HIV DRMs. We show how the occurrence of a drug-resistance mutation is contingent on favorable epistatic interactions with many other residues of the sequence background and that subsequent mutations entrench DRMs. We measure the time-autocorrelation of fluctuations in the likelihood of DRMs due to epistatic coupling with the sequence background, which reveals the presence of two evolutionary processes controlling DRM kinetics with two distinct time scales. Further analysis of waiting times for the evolutionary trapping effect to reverse reveals that the sequences which entrench (trap) a DRM are responsible for the slower time scale. We also quantify the overall strength of epistatic effects on the evolutionary kinetics for different mutations and show these are much larger for DRM positions than polymorphic positions, and we also show that trapping of a DRM is often caused by the collective effect of many accessory mutations, rather than a few strongly coupled ones, suggesting the importance of multiresidue sequence variations in HIV evolution. The analysis presented here provides a framework to explore the kinetic pathways through which viral proteins like HIV evolve under drug-selection pressure.

Structure-based virtual screening workflow to identify antivirals targeting HIV-1 capsid

We have identified novel HIV-1 capsid inhibitors targeting the PF74 binding site. Acting as the building block of the HIV-1 capsid core, the HIV-1 capsid protein plays an important role in the viral life cycle and is an attractive target for antiviral development. A structure-based virtual screening workflow for hit identification was employed, which includes docking 1.6 million commercially-available drug-like compounds from the ZINC database to the capsid dimer, followed by applying two absolute binding free energy (ABFE) filters on the 500 top-ranked molecules from docking. The first employs the Binding Energy Distribution Analysis Method (BEDAM) in implicit solvent. The top-ranked compounds are then refined using the Double Decoupling method in explicit solvent. Both docking and BEDAM refinement were carried out on the IBM World Community Grid as part of the FightAIDS@Home project. Using this virtual screening workflow, we identified 24 molecules with calculated binding free energies between − 6 and − 12 kcal/mol. We performed thermal shift assays on these molecules to examine their potential effects on the stability of HIV-1 capsid hexamer and found that two compounds, ZINC520357473 and ZINC4119064 increased the melting point of the latter by 14.8 °C and 33 °C, respectively. These results support the conclusion that the two ZINC compounds are primary hits targeting the capsid dimer interface. Our simulations also suggest that the two hit molecules may bind at the capsid dimer interface by occupying a new sub-pocket that has not been exploited by existing CA inhibitors. The possible causes for why other top-scored compounds suggested by ABFE filters failed to show measurable activity are discussed.

Multiscale Approach for Computing Gated Ligand Binding from Molecular Dynamics and Brownian Dynamics Simulations

We develop an approach to characterize the effects of gating by a multiconformation protein consisting of macrostate conformations that are either accessible or inaccessible to ligand binding. We first construct a Markov state model of the apo-protein from atomistic molecular dynamics simulations from which we identify macrostates and their conformations, compute their relative macrostate populations and interchange kinetics, and structurally characterize them in terms of ligand accessibility. We insert the calculated first-order rate constants for conformational transitions into a multistate gating theory from which we derive a gating factor γ that quantifies the degree of conformational gating. Applied to HIV-1 protease, our approach yields a kinetic network of three accessible (semi-open, open, and wide-open) and two inaccessible (closed and a newly identified, "parted") macrostate conformations. The parted conformation sterically partitions the active site, suggesting a possible role in product release. We find that the binding kinetics of drugs and drug-like inhibitors to HIV-1 protease falls in the slow gating regime. However, because γ = 0.75, conformational gating only modestly slows ligand binding. Brownian dynamics simulations of the diffusional association of eight inhibitors to the protease─having a wide range of experimental association constants (∼10⁴-10¹⁰ M^-1 s^-1)─yields gated rate constants in the range of ∼0.5-5.7 × 10⁸ M^-1 s^-1. This indicates that, whereas the association rate of some inhibitors could be described by the model, for many inhibitors either subsequent conformational transitions or alternate binding mechanisms may be rate-limiting. For systems known to be modulated by conformational gating, the approach could be scaled computationally efficiently to screen association kinetics for a large number of ligands.

Classification and Design of HIV-1 Integrase Inhibitors Based on Machine Learning

A key enzyme in human immunodeficiency virus type 1 (HIV-1) life cycle, integrase (IN) aids the integration of viral DNA into the host DNA, which has become an ideal target for the development of anti-HIV drugs. A total of 1785 potential HIV-1 IN inhibitors were collected from the databases of ChEMBL, Binding Database, DrugBank, and PubMed, as well as from 40 references. The database was divided into the training set and test set by random sampling. By exploring the correlation between molecular descriptors and inhibitory activity, it is found that the classification and specific activity data of inhibitors can be more accurately predicted by the combination of molecular descriptors and molecular fingerprints. The calculation of molecular fingerprint descriptor provides the additional substructure information to improve the prediction ability. Based on the training set, two machine learning methods, the recursive partition (RP) and naive Bayes (NB) models, were used to build the classifiers of HIV-1 IN inhibitors. Through the test set verification, the RP technique accurately predicted 82.5% inhibitors and 86.3% noninhibitors. The NB model predicted 88.3% inhibitors and 87.2% noninhibitors with correlation coefficient of 85.2%. The results show that the prediction performance of NB model is slightly better than that of RP, and the key molecular segments are also obtained. Additionally, CoMFA and CoMSIA models with good activity prediction ability both were constructed by exploring the structure-activity relationship, which is helpful for the design and optimization of HIV-1 IN inhibitors.

Increasing the Sampling Efficiency of Protein Conformational Change by Combining a Modified Replica Exchange Molecular Dynamics and Normal Mode Analysis

Understanding conformational change at an atomic level is significant when determining a protein functional mechanism. Replica exchange molecular dynamics (REMD) is a widely used enhanced sampling method to explore protein conformational space. However, REMD with an explicit solvent model requires huge computational resources, immensely limiting its application. In this study, a variation of parallel tempering metadynamics (PTMetaD) with the omission of solvent-solvent interactions in exchange attempts and the use of low-frequency modes calculated by normal-mode analysis (NMA) as collective variables (CVs), namely ossPTMetaD, is proposed with the aim to accelerate MD simulations simultaneously in temperature and geometrical spaces. For testing the performance of ossPTMetaD, five protein systems with diverse biological functions and motion patterns were selected, including large-scale domain motion (AdK), flap movement (HIV-1 protease and BACE1), and DFG-motif flip in kinases (p38α and c-Abl). The simulation results showed that ossPTMetaD requires much fewer numbers of replicas than temperature REMD (T-REMD) with a reduction of ∼70% to achieve a similar exchange ratio. Although it does not obey the detailed balance condition, ossPTMetaD provides consistent results with T-REMD and experimental data. The high accessibility of the large conformational change of protein systems by ossPTMetaD, especially in simulating the very challenging DFG-motif flip of protein kinases, demonstrated its high efficiency and robustness in the characterization of the large-scale protein conformational change pathway and associated free energy profile.

Drug Re-positioning Studies for Novel HIV-1 Inhibitors Using Binary QSAR Models and Multi-target-driven In Silico Studies

Current antiretroviral therapies against HIV involve the usage of at least two drugs that target different stages of HIV life cycle. However, potential drug interactions and side effects pose a problem. A promising concept for complex disease treatment is 'one molecule-multiple target' approach to overcome undesired effects of multiple drugs. Additionally, it is beneficial to consider drug re-purposing due to the cost of taking a drug into the market. Taking these into account, here potential anti-HIV compounds are suggested by virtually screening small approved drug molecules and clinical candidates. Initially, binary QSAR models are used to predict the therapeutic activity of around 7900 compounds against HIV and to predict the toxicity of molecules with high therapeutic activities. Selected compounds are considered for molecular docking studies against two targets, HIV-1 protease enzyme, and chemokine co-receptor CCR5. The top docking poses for all 549 molecules are then subjected to short (1 ns) individual molecular dynamics (MD) simulations and they are ranked based on their calculated relative binding free energies. Finally, 25 molecules are selected for long (200 ns) MD simulations, and 5 molecules are suggested as promising multi-target HIV agents. The results of this study may open new avenues for the designing of new dual HIV-1 inhibitor scaffolds.

Absolute Protein Binding Free Energy Simulations for Ligands with Multiple Poses, a Thermodynamic Path That Avoids Exhaustive Enumeration of the Poses

We propose a free energy calculation method for receptor-ligand binding, which have multiple binding poses that avoids exhaustive enumeration of the poses. For systems with multiple binding poses, the standard procedure is to enumerate orientations of the binding poses, restrain the ligand to each orientation, and then, calculate the binding free energies for each binding pose. In this study, we modify a part of the thermodynamic cycle in order to sample a broader conformational space of the ligand in the binding site. This modification leads to more accurate free energy calculation without performing separate free energy simulations for each binding pose. We applied our modification to simple model host-guest systems as a test, which have only two binding poses, by using a single decoupling method (SDM) in implicit solvent. The results showed that the binding free energies obtained from our method without knowing the two binding poses were in good agreement with the benchmark results obtained by explicit enumeration of the binding poses. Our method is applicable to other alchemical binding free energy calculation methods such as the double decoupling method (DDM) in explicit solvent. We performed a calculation for a protein-ligand system with explicit solvent using our modified thermodynamic path. The results of the free energy simulation along our modified path were in good agreement with the results of conventional DDM, which requires a separate binding free energy calculation for each of the binding poses of the example of phenol binding to T4 lysozyme in explicit solvent. © 2019 Wiley Periodicals, Inc.

On the removal of initial state bias from simulation data

Classical atomistic simulations of biomolecules play an increasingly important role in molecular life science. The structure of current computing architectures favors methods that run multiple trajectories at once without requiring extensive communication between them. Many advanced sampling strategies in the field fit this mold. These approaches often rely on an adaptive logic and create ensembles of comparatively short trajectories whose starting points are not distributed according to the correct Boltzmann weights. This type of bias is notoriously difficult to remove, and Markov state models (MSMs) are one of the few strategies available for recovering the correct kinetics and thermodynamics from these ensembles of trajectories. In this contribution, we analyze the performance of MSMs in the thermodynamic reweighting task for a hierarchical set of systems. We show that MSMs can be rigorous tools to recover the correct equilibrium distribution for systems of sufficiently low dimensionality. This is conditional upon not tampering with local flux imbalances found in the data. For a real-world application, we find that a pure likelihood-based inference of the transition matrix produces the best results. The removal of the bias is incomplete, however, and for this system, all tested MSMs are outperformed by an alternative albeit less general approach rooted in the ideas of statistical resampling. We conclude by formulating some recommendations for how to address the reweighting issue in practice.

Conformational and dynamical basis for cross-reactivity observed between anti HIV-1 protease antibody with protease and an epitope peptide from it

F11.2.32 is a monoclonal antibody raised against HIV-1 protease and it inhibits protease activity. While the structure of the epitope peptide in complex with the antibody is known, how protease interacts with the antibody is not known. In this study, we model the conformational features of the free and bound epitope peptide and protease-antibody interactions. We find through our simulations, that the free epitope peptide P36-46 samples conformations akin to the bound conformation of the peptide in complex with the Ab, with a β-turn conformation sampled by the ³⁸LPGR⁴¹ sequence highlighting the role of inherent conformational preferences of the peptide. Further, to determine the interactions present between the protease and antibody, we docked the protease in its conformation observed in the crystal structure, onto the antibody and simulated the dynamics of the complex in explicit water. We have identified the key residues involved in hydrogen-bond interactions and salt-bridges in Ag-Ab complex and examined the role of CDR flexibility in binding different conformations of the same epitope sequence in peptide and protein antigens. Thus, our results provide the basis for understanding the cross-reactivity observed between the antibody with protease and the epitope peptide from it.

All-atom virus simulations

The constant threat of viral disease can be combated by the development of novel vaccines and therapeutics designed to disrupt key features of virus structure or infection cycle processes. Such development relies on high-resolution characterization of viruses and their dynamical behaviors, which are often challenging to obtain solely by experiment. In response, all-atom molecular dynamics simulations are widely leveraged to study the structural components of viruses, leading to some of the largest simulation endeavors undertaken to date. The present work reviews exemplary all-atom simulation work on viruses, as well as progress toward simulating entire virions.

Committor of elementary reactions on multistate systems

In our study, we extend the committor concept on multi-minima systems, where more than one reaction may proceed, but the feasible data evaluation needs the projection onto partial reactions. The elementary reaction committor and the corresponding probability density of the reactive trajectories are defined and calculated on a three-hole two-dimensional model system explored by single-particle Langevin dynamics. We propose a method to visualize more elementary reaction committor functions or probability densities of reactive trajectories on a single plot that helps to identify the most important reaction channels and the nonreactive domains simultaneously. We suggest a weighting for the energy-committor plots that correctly shows the limits of both the minimal energy path and the average energy concepts. The methods also performed well on the analysis of molecular dynamics trajectories of 2-chlorobutane, where an elementary reaction committor, the probability densities, the potential energy/committor, and the free-energy/committor curves are presented.

Self-Assembly Pathways of β-Sheet-Rich Amyloid-β(1-40) Dimers: Markov State Model Analysis on Millisecond Hybrid-Resolution Simulations

Early oligomerization during amyloid-β (Aβ) aggregation is essential for Aβ neurotoxicity. Understanding how unstructured Aβs assemble into oligomers, especially those rich in β-sheets, is essential but remains challenging as the assembly process is too transient for experimental characterization and too slow for molecular dynamics simulations. So far, atomic simulations are limited only to studies of either oligomer structures or assembly pathways for short Aβ segments. To overcome the computational challenge, we combine in this study a hybrid-resolution model and adaptive sampling techniques to perform over 2.7 ms of simulations of formation of full-length Aβ40 dimers that are the earliest toxic oligomeric species. The Markov state model is further employed to characterize the transition pathways and associated kinetics. Our results show that for two major forms of β-sheet-rich structures reported experimentally, the corresponding assembly mechanisms are markedly different. Hairpin-containing structures are formed by direct binding of soluble Aβ in β-hairpin-like conformations. Formation of parallel, in-register structures resembling fibrils occurs ∼100-fold more slowly and involves a rapid encounter of Aβ in arbitrary conformations followed by a slow structural conversion. The structural conversion proceeds via diverse pathways but always requires transient unfolding of encounter complexes. We find that the transition kinetics could be affected differently by intra-/intermolecular interactions involving individual residues in a conformation-dependent manner. In particular, the interactions involving Aβ's N-terminal part promote the assembly into hairpin-containing structures but delay the formation of fibril-like structures, thus explaining puzzling observations reported previously regarding the roles of this region in the early assembly process.

Stratified UWHAM and Its Stochastic Approximation for Multicanonical Simulations Which Are Far from Equilibrium

We describe a new analysis tool called Stratified unbinned Weighted Histogram Analysis Method (Stratified-UWHAM), which can be used to compute free energies and expectations from a multicanonical ensemble when a subset of the parallel simulations is far from being equilibrated because of barriers between free energy basins which are only rarely (or never) crossed at some states. The Stratified-UWHAM equations can be obtained in the form of UWHAM equations but with an expanded set of states. We also provide a stochastic solver, Stratified RE-SWHAM, for Stratified-UWHAM to remove its computational bottleneck. Stratified-UWHAM and Stratified RE-SWHAM are applied to study three test topics: the free energy landscape of alanine dipeptide, the binding affinity of a host-guest binding complex, and path sampling for a two-dimensional double well potential. The examples show that when some of the parallel simulations are only locally equilibrated, the estimates of free energies and equilibrium distributions provided by the conventional UWHAM (or MBAR) solutions exhibit considerable biases, but the estimates provided by Stratified-UWHAM and Stratified RE-SWHAM agree with the benchmark very well. Lastly, we discuss features of the Stratified-UWHAM approach which is based on coarse-graining in relation to two other maximum likelihood-based methods which were proposed recently, that also coarse-grain the multicanonical data.

Impact of Genetic Variations in HIV-1 Tat on LTR-Mediated Transcription via TAR RNA Interaction

HIV-1 evades host defense through mutations and recombination events, generating numerous variants in an infected patient. These variants with an undiminished virulence can multiply rapidly in order to progress to AIDS. One of the targets to intervene in HIV-1 replication is the trans-activator of transcription (Tat), a major regulatory protein that transactivates the long terminal repeat promoter through its interaction with trans-activation response (TAR) RNA. In this study, HIV-1 infected patients (n = 120) from North India revealed Ser46Phe (20%) and Ser61Arg (2%) mutations in the Tat variants with a strong interaction toward TAR leading to enhanced transactivation activities. Molecular dynamics simulation data verified that the variants with this mutation had a higher binding affinity for TAR than both the wild-type Tat and other variants that lacked Ser46Phe and Ser61Arg. Other mutations in Tat conferred varying affinities for TAR interaction leading to differential transactivation abilities. This is the first report from North India with a clinical validation of CD4 counts to demonstrate the influence of Tat genetic variations affecting the stability of Tat and its interaction with TAR. This study highlights the co-evolution pattern of Tat and predominant nucleotides for Tat activity, facilitating the identification of genetic determinants for the attenuation of viral gene expression.

Mechanism of the Association Pathways for a Pair of Fast and Slow Binding Ligands of HIV-1 Protease

Equilibrium constants, together with kinetic rate constants of binding, are key factors in the efficacy and safety of drug compounds, informing drug design. However, the association pathways of protein-ligand binding, which contribute to their kinetic behaviors, are little understood. In this work, we used unbiased all-atom molecular dynamics (MD) simulations with an explicit solvent model to study the association processes of protein-ligand binding. Using the HIV protease (HIVp)-xk263 and HIVp-ritonavir protein-ligand systems as cases, we observed that ligand association is a multistep process involving diffusion, localization, and conformational rearrangements of the protein, ligand, and water molecules. Moreover, these two ligands preferred different routes of binding, which reflect two well-known binding mechanisms: induced-fit and conformation selection models. Our study shows that xk263 has a stronger capacity for desolvating surrounding water molecules, thereby inducing a semiopen conformation of the HIVp flaps (induced-fit model). In contrast, the slow dehydration characteristic of ritonavir allows for gradual association with the binding pocket of HIVp when the protein's flap conformation is fully open (conformation selection model). By studying the mechanism of ligand association and understanding the role of solvent molecules during the binding event, we can obtain a different perspective on the mechanism of macromolecule recognition, providing insights into drug discovery.

Structural insights into HIV-1 protease flap opening processes and key intermediates

The study provided an integrated view of the transition pathway of the flap opening of HIV-1 protease using MD simulation.

Effect of allosteric molecules on structure and drug affinity of HIV-1 protease by molecular dynamics simulations

Recent experiments show that small molecules can bind onto the allosteric sites of HIV-1 protease (PR), which provides a starting point for developing allosteric inhibitors. However, the knowledge of the effect of such binding on the structural dynamics and binding free energy of the active site inhibitor and PR is still lacking. Here, we report 200ns long molecular dynamics simulation results to gain insight into the influences of two allosteric molecules (1H-indole-6-carboxylic acid, 1F1 and 2-methylcyclohexano, 4D9). The simulations demonstrate that both allosteric molecules change the PR conformation and stabilize the structures of PR and the inhibitor; the residues of the flaps are sensitive to the allosteric molecules and the flexibility of the residues is pronouncedly suppressed; the additions of the small molecules to the allosteric sites strengthen the binding affinities of 3TL-PR by about 12-15kal/mol in the binding free energy, which mainly arises from electrostatic term. Interestingly, it is found that the action mechanisms of 1F1 and 4D9 are different, the former behaviors like a doorman that keeps the inhibitor from escape and makes the flaps (door) partially open; the latter is like a wedge that expands the allosteric space and meanwhile closes the flaps. Our data provide a theoretical support for designing the allosteric inhibitor.

Exploration of gated ligand binding recognizes an allosteric site for blocking FABP4-protein interaction

Fatty acid binding protein 4 (FABP4), reversibly binding to fatty acids and other lipids with high affinities, is a potential target for treatment of cancers. The binding site of FABP4 is buried in an interior cavity and thereby ligand binding/unbinding is coupled with opening/closing of FABP4. It is a difficult task both experimentally and computationally to illuminate the entry or exit pathway, especially with the conformational gating. In this report we combine extensive computer simulations, clustering analysis, and the Markov state model to investigate the binding mechanism of FABP4 and troglitazone. Our simulations capture spontaneous binding and unbinding events as well as the conformational transition of FABP4 between the open and closed states. An allosteric binding site on the protein surface is recognized for the development of novel FABP4 inhibitors. The binding affinity is calculated and compared with the experimental value. The kinetic analysis suggests that ligand residence on the protein surface may delay the binding process. Overall, our results provide a comprehensive picture of ligand diffusion on the protein surface, ligand migration into the buried cavity, and the conformational change of FABP4 at an atomic level.

Transition path theory analysis of c-Src kinase activation

Nonreceptor tyrosine kinases of the Src family are large multidomain allosteric proteins that are crucial to cellular signaling pathways. In a previous study, we generated a Markov state model (MSM) to simulate the activation of c-Src catalytic domain, used as a prototypical tyrosine kinase. The long-time kinetics of transition predicted by the MSM was in agreement with experimental observations. In the present study, we apply the framework of transition path theory (TPT) to the previously constructed MSM to characterize the main features of the activation pathway. The analysis indicates that the activating transition, in which the activation loop first opens up followed by an inward rotation of the αC-helix, takes place via a dense set of intermediate microstates distributed within a fairly broad "transition tube" in a multidimensional conformational subspace connecting the two end-point conformations. Multiple microstates with negligible equilibrium probabilities carry a large transition flux associated with the activating transition, which explains why extensive conformational sampling is necessary to accurately determine the kinetics of activation. Our results suggest that the combination of MSM with TPT provides an effective framework to represent conformational transitions in complex biomolecular systems.

Connecting free energy surfaces in implicit and explicit solvent: an efficient method to compute conformational and solvation free energies

The ability to accurately model solvent effects on free energy surfaces is important for understanding many biophysical processes including protein folding and misfolding, allosteric transitions, and protein–ligand binding. Although all-atom simulations in explicit solvent can provide an accurate model for biomolecules in solution, explicit solvent simulations are hampered by the slow equilibration on rugged landscapes containing multiple basins separated by barriers. In many cases, implicit solvent models can be used to significantly speed up the conformational sampling; however, implicit solvent simulations do not fully capture the effects of a molecular solvent, and this can lead to loss of accuracy in the estimated free energies. Here we introduce a new approach to compute free energy changes in which the molecular details of explicit solvent simulations are retained while also taking advantage of the speed of the implicit solvent simulations. In this approach, the slow equilibration in explicit solvent, due to the long waiting times before barrier crossing, is avoided by using a thermodynamic cycle which connects the free energy basins in implicit solvent and explicit solvent using a localized decoupling scheme. We test this method by computing conformational free energy differences and solvation free energies of the model system alanine dipeptide in water. The free energy changes between basins in explicit solvent calculated using fully explicit solvent paths agree with the corresponding free energy differences obtained using the implicit/explicit thermodynamic cycle to within 0.3 kcal/mol out of ∼3 kcal/mol at only ∼8% of the computational cost. We note that WHAM methods can be used to further improve the efficiency and accuracy of the implicit/explicit thermodynamic cycle.

Investigation of Structural Dynamics of Enzymes and Protonation States of Substrates Using Computational Tools

This review discusses the use of molecular modeling tools, together with existing experimental findings, to provide a complete atomic-level description of enzyme dynamics and function. We focus on functionally relevant conformational dynamics of enzymes and the protonation states of substrates. The conformational fluctuations of enzymes usually play a crucial role in substrate recognition and catalysis. Protein dynamics can be altered by a tiny change in a molecular system such as different protonation states of various intermediates or by a significant perturbation such as a ligand association. Here we review recent advances in applying atomistic molecular dynamics (MD) simulations to investigate allosteric and network regulation of tryptophan synthase (TRPS) and protonation states of its intermediates and catalysis. In addition, we review studies using quantum mechanics/molecular mechanics (QM/MM) methods to investigate the protonation states of catalytic residues of β-Ketoacyl ACP synthase I (KasA). We also discuss modeling of large-scale protein motions for HIV-1 protease with coarse-grained Brownian dynamics (BD) simulations.

Role of Conformational Motions in Enzyme Function: Selected Methodologies and Case Studies

It is now common knowledge that enzymes are mobile entities relying on complex atomic-scale dynamics and coordinated conformational events for proper ligand recognition and catalysis. However, the exact role of protein dynamics in enzyme function remains either poorly understood or difficult to interpret. This mini-review intends to reconcile biophysical observations and biological significance by first describing a number of common experimental and computational methodologies employed to characterize atomic-scale residue motions on various timescales in enzymes, and second by illustrating how the knowledge of these motions can be used to describe the functional behavior of enzymes and even act upon it. Two biologically relevant examples will be highlighted, namely the HIV-1 protease and DNA polymerase β enzyme systems.

Flap Dynamics in Aspartic Proteases: A Computational Perspective

Recent advances in biochemistry and drug design have placed proteases as one of the critical target groups for developing novel small-molecule inhibitors. Among all proteases, aspartic proteases have gained significant attention due to their role in HIV/AIDS, malaria, Alzheimer's disease, etc. The binding cleft is covered by one or two β-hairpins (flaps) which need to be opened before a ligand can bind. After binding, the flaps close to retain the ligand in the active site. Development of computational tools has improved our understanding of flap dynamics and its role in ligand recognition. In the past decade, several computational approaches, for example molecular dynamics (MD) simulations, coarse-grained simulations, replica-exchange molecular dynamics (REMD) and metadynamics, have been used to understand flap dynamics and conformational motions associated with flap movements. This review is intended to summarize the computational progress towards understanding the flap dynamics of proteases and to be a reference for future studies in this field.

A stochastic spatio-temporal (SST) model to study cell-to-cell variability in HIV-1 infection

Although HIV viremia in infected patients proceeds in a manner that may be accounted for by deterministic mathematical models, single virus-cell encounters following initial HIV exposure result in a variety of outcomes, only one of which results in a productive infection. The development of single molecule tracking techniques in living cells allows studies of intracellular transport of HIV, but it remains less clear what its impact may be on viral integration efficiency. Here, we present a stochastic intracellular mathematical model of HIV replication that incorporates microtubule transport of viral components. Using this model, we could study single round infections and observe how viruses entering cells reach one of three potential fates - degradation of the viral RNA genome, formation of LTR circles, or successful integration and establishment of a provirus. Our model predicts global trafficking properties, such as the probability and the mean time for a HIV viral particle to reach the nuclear pore. Interestingly, our model predicts that trafficking determines neither the probability or time of provirus establishment - instead, they are a function of vRNA degradation and reverse transcription reactions. Thus, our spatio-temporal model provides novel insights into the HIV infection process and may constitute a useful tool for the identification of promising drug targets.

Unique Flap Conformation in an HIV-1 Protease with High-Level Darunavir Resistance

Darunavir (DRV) is one of the most powerful protease inhibitors (PIs) for treating human immunodeficiency virus type-1 (HIV-1) infection and presents a high genetic barrier to the generation of resistant viruses. However, DRV-resistant HIV-1 infrequently emerges from viruses exhibiting resistance to other protease inhibitors. To address this resistance, researchers have gathered genetic information on DRV resistance. In contrast, few structural insights into the mechanism underlying DRV resistance are available. To elucidate this mechanism, we determined the crystal structure of the ligand-free state of a protease with high-level DRV resistance and six DRV resistance-associated mutations (including I47V and I50V), which we generated by in vitro selection. This crystal structure showed a unique curling conformation at the flap regions that was not found in the previously reported ligand-free protease structures. Molecular dynamics simulations indicated that the curled flap conformation altered the flap dynamics. These results suggest that the preference for a unique flap conformation influences DRV binding. These results provide new structural insights into elucidating the molecular mechanism of DRV resistance and aid to develop PIs effective against DRV-resistant viruses.

Enhanced conformational sampling technique provides an energy landscape view of large-scale protein conformational transitions

A novel in silico approach (NMA–ITS) is introduced to rapidly and effectively sample the configuration space and give quantitative data for exploring the conformational changes of proteins.

Extracting intrinsic dynamic parameters of biomolecular folding from single-molecule force spectroscopy experiments

Single-molecule studies in which a mechanical force is transmitted to the molecule of interest and the molecular extension or position is monitored as a function of time are versatile tools for probing the dynamics of protein folding, stepping of molecular motors, and other biomolecular processes involving activated barrier crossing. One complication in interpreting such studies, however, is the fact that the typical size of a force probe (e.g., a dielectric bead in optical tweezers or the atomic force microscope tip/cantilever assembly) is much larger than the molecule itself, and so the observed molecular motion is affected by the hydrodynamic drag on the probe. This presents the experimenter with a nontrivial task of deconvolving the intrinsic molecular parameters, such as the intrinsic free energy barrier and the effective diffusion coefficient exhibited while crossing the barrier from the experimental signal. Here we focus on the dynamical aspect of this task and show how the intrinsic diffusion coefficient along the molecular reaction coordinate can be inferred from single-molecule measurements of the rates of biomolecular folding and unfolding. We show that the feasibility of accomplishing this task is strongly dependent on the relationship between the intrinsic molecular elasticity and that of the linker connecting the molecule to the force probe and identify the optimal range of instrumental parameters allowing determination of instrument-free molecular dynamics.

Switches of hydrogen bonds during ligand-protein association processes determine binding kinetics

Revealing the processes of ligand-protein associations deepens our understanding of molecular recognition and binding kinetics. Hydrogen bonds (H-bonds) play a crucial role in optimizing ligand-protein interactions and ligand specificity. In addition to the formation of stable H-bonds in the final bound state, the formation of transient H-bonds during binding processes contributes binding kinetics that define a ligand as a fast or slow binder, which also affects drug action. However, the effect of forming the transient H-bonds on the kinetic properties is little understood. Guided by results from coarse-grained Brownian dynamics simulations, we used classical molecular dynamics simulations in an implicit solvent model and accelerated molecular dynamics simulations in explicit waters to show that the position and distribution of the H-bond donor or acceptor of a drug result in switching intermolecular and intramolecular H-bond pairs during ligand recognition processes. We studied two major types of HIV-1 protease ligands: a fast binder, xk263, and a slow binder, ritonavir. The slow association rate in ritonavir can be attributed to increased flexibility of ritonavir, which yields multistep transitions and stepwise entering patterns and the formation and breaking of complex H-bond pairs during the binding process. This model suggests the importance of conversions of spatiotemporal H-bonds during the association of ligands and proteins, which helps in designing inhibitors with preferred binding kinetics.

Molecular dynamics, monte carlo simulations, and langevin dynamics: a computational review

Macromolecular structures, such as neuraminidases, hemagglutinins, and monoclonal antibodies, are not rigid entities. Rather, they are characterised by their flexibility, which is the result of the interaction and collective motion of their constituent atoms. This conformational diversity has a significant impact on their physicochemical and biological properties. Among these are their structural stability, the transport of ions through the M2 channel, drug resistance, macromolecular docking, binding energy, and rational epitope design. To assess these properties and to calculate the associated thermodynamical observables, the conformational space must be efficiently sampled and the dynamic of the constituent atoms must be simulated. This paper presents algorithms and techniques that address the abovementioned issues. To this end, a computational review of molecular dynamics, Monte Carlo simulations, Langevin dynamics, and free energy calculation is presented. The exposition is made from first principles to promote a better understanding of the potentialities, limitations, applications, and interrelations of these computational methods.

Distinguishing binders from false positives by free energy calculations: fragment screening against the flap site of HIV protease

Molecular docking is a powerful tool used in drug discovery and structural biology for predicting the structures of ligand–receptor complexes. However, the accuracy of docking calculations can be limited by factors such as the neglect of protein reorganization in the scoring function; as a result, ligand screening can produce a high rate of false positive hits. Although absolute binding free energy methods still have difficulty in accurately rank-ordering binders, we believe that they can be fruitfully employed to distinguish binders from nonbinders and reduce the false positive rate. Here we study a set of ligands that dock favorably to a newly discovered, potentially allosteric site on the flap of HIV-1 protease. Fragment binding to this site stabilizes a closed form of protease, which could be exploited for the design of allosteric inhibitors. Twenty-three top-ranked protein–ligand complexes from AutoDock were subject to the free energy screening using two methods, the recently developed binding energy analysis method (BEDAM) and the standard double decoupling method (DDM). Free energy calculations correctly identified most of the false positives (≥83%) and recovered all the confirmed binders. The results show a gap averaging ≥3.7 kcal/mol, separating the binders and the false positives. We present a formula that decomposes the binding free energy into contributions from the receptor conformational macrostates, which provides insights into the roles of different binding modes. Our binding free energy component analysis further suggests that improving the treatment for the desolvation penalty associated with the unfulfilled polar groups could reduce the rate of false positive hits in docking. The current study demonstrates that the combination of docking with free energy methods can be very useful for more accurate ligand screening against valuable drug targets.

Fibril elongation by Aβ(17-42): kinetic network analysis of hybrid-resolution molecular dynamics simulations

A critical step of β-amyloid fibril formation is fibril elongation in which amyloid-β monomers undergo structural transitions to fibrillar structures upon their binding to fibril tips. The atomic detail of the structural transitions remains poorly understood. Computational characterization of the structural transitions is limited so far to short Aβ segments (5–10 aa) owing to the long time scale of Aβ fibril elongation. To overcome the computational time scale limit, we combined a hybrid-resolution model with umbrella sampling and replica exchange molecular dynamics and performed altogether ∼1.3 ms of molecular dynamics simulations of fibril elongation for Aβ_17–42. Kinetic network analysis of biased simulations resulted in a kinetic model that encompasses all Aβ segments essential for fibril formation. The model not only reproduces key properties of fibril elongation measured in experiments, including Aβ binding affinity, activation enthalpy of Aβ structural transitions and a large time scale gap (τ_lock/τ_dock = 10³–10⁴) between Aβ binding and its structural transitions, but also reveals detailed pathways involving structural transitions not seen before, namely, fibril formation both in hydrophobic regions L17-A21 and G37-A42 preceding fibril formation in hydrophilic region E22-A30. Moreover, the model identifies as important kinetic intermediates strand–loop–strand (SLS) structures of Aβ monomers, long suspected to be related to fibril elongation. The kinetic model suggests further that fibril elongation arises faster at the fibril tip with exposed L17-A21, rather than at the other tip, explaining thereby unidirectional fibril growth observed previously in experiments.

An allosteric modulator of HIV-1 protease shows equipotent inhibition of wild-type and drug-resistant proteases

NMR and MD simulations have demonstrated that the flaps of HIV-1 protease (HIV-1p) adopt a range of conformations that are coupled with its enzymatic activity. Previously, a model was created for an allosteric site located between the flap and the core of HIV-1p, called the Eye site (Biopolymers2008, 89, 643−65218381626). Here, results from our first study were combined with a ligand-based, lead-hopping method to identify a novel compound (NIT). NIT inhibits HIV-1p, independent of the presence of an active-site inhibitor such as pepstatin A. Assays showed that NIT acts on an allosteric site other than the dimerization interface. MD simulations of the ligand–protein complex show that NIT stably binds in the Eye site and restricts the flaps. That bound state of NIT is consistent with a crystal structure of similar fragments bound in the Eye site (Chem. Biol. Drug Des.2010, 75, 257−26820659109). Most importantly, NIT is equally potent against wild-type and a multidrug-resistant mutant of HIV-1p, which highlights the promise of allosteric inhibitors circumventing existing clinical resistance.

Quantifying the Sources of Kinetic Frustration in Folding Simulations of Small Proteins

Experiments and atomistic simulations of polypeptides have revealed structural intermediates that promote or inhibit conformational transitions to the native state during folding. We invoke a concept of “kinetic frustration” to quantify the prevalence and impact of these behaviors on folding rates within a large set of atomistic simulation data for 10 fast-folding proteins, where each protein’s conformational space is represented as a Markov state model of conformational transitions. Our graph theoretic approach addresses what conformational features correlate with folding inhibition and therefore permits comparison among features within a single protein network and also more generally between proteins. Nonnative contacts and nonnative secondary structure formation can thus be quantitatively implicated in inhibiting folding for several of the tested peptides.

Molecular dynamics of the proline switch and its role in Crk signaling

The Crk adaptor proteins play a central role as a molecular timer for the formation of protein complexes including various growth and differentiation factors. The loss of regulation of Crk results in many kinds of cancers. A self-regulatory mechanism for Crk was recently proposed, which involves domain–domain rearrangement. It is initiated by a cis–trans isomerization of a specific proline residue (Pro238 in chicken Crk II) and can be accelerated by Cyclophilin A. To understand how the proline switch controls the autoinhibition at the molecular level, we performed large-scale molecular dynamics and metadynamics simulations in the context of short peptides and multidomain constructs of chicken Crk II. We found that the equilibrium and kinetic properties of the macrostates are regulated not only by the local environments of specified prolines but also by the global organization of multiple domains. We observe the two macrostates (cis closed/autoinhibited and trans open/uninhibited) consistent with NMR experiments and predict barriers. We also propose an intermediate state, the trans closed state, which interestingly was reported to be a prevalent state in human Crk II. The existence of this macrostate suggests that the rate of switching off the autoinhibition by Cyp A may be limited by the relaxation rate of this intermediate state.

Activation pathway of Src kinase reveals intermediate states as targets for drug design

Unregulated activation of Src kinases leads to aberrant signaling, uncontrolled growth, and differentiation of cancerous cells. Reaching a complete mechanistic understanding of large scale conformational transformations underlying the activation of kinases could greatly help in the development of therapeutic drugs for the treatment of these pathologies. In principle, the nature of conformational transition could be modeled in silico via atomistic molecular dynamics simulations, although this is very challenging due to the long activation timescales. Here, we employ a computational paradigm that couples transition pathway techniques and Markov state model-based massively distributed simulations for mapping the conformational landscape of c-src tyrosine kinase. The computations provide the thermodynamics and kinetics of kinase activation for the first time, and help identify key structural intermediates. Furthermore, the presence of a novel allosteric site in an intermediate state of c-src that could be potentially utilized for drug design is predicted.

Virtual screening of integrase inhibitors by large scale binding free energy calculations: the SAMPL4 challenge

As part of the SAMPL4 blind challenge, filtered AutoDock Vina ligand docking predictions and large scale binding energy distribution analysis method binding free energy calculations have been applied to the virtual screening of a focused library of candidate binders to the LEDGF site of the HIV integrase protein. The computational protocol leveraged docking and high level atomistic models to improve enrichment. The enrichment factor of our blind predictions ranked best among all of the computational submissions, and second best overall. This work represents to our knowledge the first example of the application of an all-atom physics-based binding free energy model to large scale virtual screening. A total of 285 parallel Hamiltonian replica exchange molecular dynamics absolute protein-ligand binding free energy simulations were conducted starting from docked poses. The setup of the simulations was fully automated, calculations were distributed on multiple computing resources and were completed in a 6-weeks period. The accuracy of the docked poses and the inclusion of intramolecular strain and entropic losses in the binding free energy estimates were the major factors behind the success of the method. Lack of sufficient time and computing resources to investigate additional protonation states of the ligands was a major cause of mispredictions. The experiment demonstrated the applicability of binding free energy modeling to improve hit rates in challenging virtual screening of focused ligand libraries during lead optimization.

How kinetics within the unfolded state affects protein folding: an analysis based on markov state models and an ultra-long MD trajectory

Understanding how kinetics in the unfolded state affects protein folding is a fundamentally important yet less well-understood issue. Here we employ three different models to analyze the unfolded landscape and folding kinetics of the miniprotein Trp-cage. The first is a 208 μs explicit solvent molecular dynamics (MD) simulation from D. E. Shaw Research containing tens of folding events. The second is a Markov state model (MSM-MD) constructed from the same ultralong MD simulation; MSM-MD can be used to generate thousands of folding events. The third is a Markov state model built from temperature replica exchange MD simulations in implicit solvent (MSM-REMD). All the models exhibit multiple folding pathways, and there is a good correspondence between the folding pathways from direct MD and those computed from the MSMs. The unfolded populations interconvert rapidly between extended and collapsed conformations on time scales ≤40 ns, compared with the folding time of ∼5 μs. The folding rates are independent of where the folding is initiated from within the unfolded ensemble. About 90% of the unfolded states are sampled within the first 40 μs of the ultralong MD trajectory, which on average explores ∼27% of the unfolded state ensemble between consecutive folding events. We clustered the folding pathways according to structural similarity into "tubes", and kinetically partitioned the unfolded state into populations that fold along different tubes. From our analysis of the simulations and a simple kinetic model, we find that, when the mixing within the unfolded state is comparable to or faster than folding, the folding waiting times for all the folding tubes are similar and the folding kinetics is essentially single exponential despite the presence of heterogeneous folding paths with nonuniform barriers. When the mixing is much slower than folding, different unfolded populations fold independently, leading to nonexponential kinetics. A kinetic partition of the Trp-cage unfolded state is constructed which reveals that different unfolded populations have almost the same probability to fold along any of the multiple folding paths. We are investigating whether the results for the kinetics in the unfolded state of the 20-residue Trp-cage is representative of larger single domain proteins.

NMR relaxation in proteins with fast internal motions and slow conformational exchange: model-free framework and Markov state simulations

Calculating NMR relaxation effects for proteins with dynamics on multiple time scales generally requires very long trajectories based on conventional molecular dynamics simulations. In this report, we have built Markov state models from multiple MD trajectories and used the resulting MSM to capture the very fast internal motions of the protein within a free energy basin on a time scale up to hundreds of picoseconds and the more than 3 orders of magnitude slower conformational exchange between macrostates. To interpret the relaxation data, we derive new equations using the model-free framework which includes two slowly exchanging macrostates, each of which also exhibits fast local motions. Using simulations of HIV-1 protease as an example, we show how the populations of slowly exchanging conformational states as well as order parameters for the different states can be determined from the NMR relaxation data.

How conformational changes can affect catalysis, inhibition and drug resistance of enzymes with induced-fit binding mechanism such as the HIV-1 protease

A central question is how the conformational changes of proteins affect their function and the inhibition of this function by drug molecules. Many enzymes change from an open to a closed conformation upon binding of substrate or inhibitor molecules. These conformational changes have been suggested to follow an induced-fit mechanism in which the molecules first bind in the open conformation in those cases where binding in the closed conformation appears to be sterically obstructed such as for the HIV-1 protease. In this article, we present a general model for the catalysis and inhibition of enzymes with induced-fit binding mechanism. We derive general expressions that specify how the overall catalytic rate of the enzymes depends on the rates for binding, for the conformational changes, and for the chemical reaction. Based on these expressions, we analyze the effect of mutations that mainly shift the conformational equilibrium on catalysis and inhibition. If the overall catalytic rate is limited by product unbinding, we find that mutations that destabilize the closed conformation relative to the open conformation increase the catalytic rate in the presence of inhibitors by a factor exp(ΔΔGC/RT) where ΔΔGC is the mutation-induced shift of the free-energy difference between the conformations. This increase in the catalytic rate due to changes in the conformational equilibrium is independent of the inhibitor molecule and, thus, may help to understand how non-active-site mutations can contribute to the multi-drug-resistance that has been observed for the HIV-1 protease. A comparison to experimental data for the non-active-site mutation L90M of the HIV-1 protease indicates that the mutation slightly destabilizes the closed conformation of the enzyme. This article is part of a Special Issue entitled: The emerging dynamic view of proteins: Protein plasticity in allostery, evolution and self-assembly.

Event detection and sub-state discovery from biomolecular simulations using higher-order statistics: application to enzyme adenylate kinase

Biomolecular simulations at millisecond and longer time-scales can provide vital insights into functional mechanisms. Because post-simulation analyses of such large trajectory datasets can be a limiting factor in obtaining biological insights, there is an emerging need to identify key dynamical events and relating these events to the biological function online, that is, as simulations are progressing. Recently, we have introduced a novel computational technique, quasi-anharmonic analysis (QAA) (Ramanathan et al., PLoS One 2011;6:e15827), for partitioning the conformational landscape into a hierarchy of functionally relevant sub-states. The unique capabilities of QAA are enabled by exploiting anharmonicity in the form of fourth-order statistics for characterizing atomic fluctuations. In this article, we extend QAA for analyzing long time-scale simulations online. In particular, we present HOST4MD--a higher-order statistical toolbox for molecular dynamics simulations, which (1) identifies key dynamical events as simulations are in progress, (2) explores potential sub-states, and (3) identifies conformational transitions that enable the protein to access those sub-states. We demonstrate HOST4MD on microsecond timescale simulations of the enzyme adenylate kinase in its apo state. HOST4MD identifies several conformational events in these simulations, revealing how the intrinsic coupling between the three subdomains (LID, CORE, and NMP) changes during the simulations. Further, it also identifies an inherent asymmetry in the opening/closing of the two binding sites. We anticipate that HOST4MD will provide a powerful and extensible framework for detecting biophysically relevant conformational coordinates from long time-scale simulations.

Molecular dynamics simulation in virus research

Virus replication in the host proceeds by chains of interactions between viral and host proteins. The interactions are deeply influenced by host immune molecules and anti-viral compounds, as well as by mutations in viral proteins. To understand how these interactions proceed mechanically and how they are influenced by mutations, one needs to know the structures and dynamics of the proteins. Molecular dynamics (MD) simulation is a powerful computational method for delineating motions of proteins at an atomic-scale via theoretical and empirical principles in physical chemistry. Recent advances in the hardware and software for biomolecular simulation have rapidly improved the precision and performance of this technique. Consequently, MD simulation is quickly extending the range of applications in biology, helping to reveal unique features of protein structures that would be hard to obtain by experimental methods alone. In this review, we summarize the recent advances in MD simulations in the study of virus–host interactions and evolution, and present future perspectives on this technique.

Structural and Entropic Allosteric Signal Transduction Strength via Correlated Motions

Allosteric signal transduction in biomacromolecules can play an essential role in their function. Internal motional correlations in proteins provide a possible communication mechanism, but the quantitative relationship between statistical correlations and allostery is unknown. Quantitative relationships between internal motional correlations and the efficiency of propagation of allosteric structural and entropic effects are introduced and validated against conformational ensembles obtained from molecular dynamics simulations. This framework can explain a range of phenomena, such as the occurrence of an allosteric entropy change in the absence of any noticeable structural change.

Quasi-anharmonic analysis reveals intermediate states in the nuclear co-activator receptor binding domain ensemble

The molten globule nuclear receptor co-activator binding domain (NCBD) of CREB binding protein (CBP) selectively recruits transcription co-activators (TCAs) during the formation of the transcription preinitiation complex. NCBD:TCA interactions have been implicated in several cancers, however, the mechanisms of NCBD:TCA recognition remain uncharacterized. NCBD:TCA intermolecular recognition has challenged traditional investigation as both NCBD and several of its corresponding TCAs are intrinsically disordered. Using 40μs of explicit solvent molecular dynamics simulations, we relate the conformational diversity of ligand-free NCBD to its bound configurations. We introduce two novel techniques to quantify the conformational heterogeneity of ligand-free NCBD, dihedral quasi-anharmonic analysis (dQAA) and hierarchical graph-based diffusive clustering. With this integrated approach we find that three of four ligand-bound states are natively accessible to the ligand-free NCBD simulations with root-mean squared deviation (RMSD) less than 2Å These conformations are accessible via diverse pathways while a rate-limiting barrier must be crossed in order to access the fourth bound state.

Conformational Transitions and Convergence of Absolute Binding Free Energy Calculations

The Binding Energy Distribution Analysis Method (BEDAM) is employed to compute the standard binding free energies of a series of ligands to a FK506 binding protein (FKBP12) with implicit solvation. Binding free energy estimates are in reasonably good agreement with experimental affinities. The conformations of the complexes identified by the simulations are in good agreement with crystallographic data, which was not used to restrain ligand orientations. The BEDAM method is based on λ -hopping Hamiltonian parallel Replica Exchange (HREM) molecular dynamics conformational sampling, the OPLS-AA/AGBNP2 effective potential, and multi-state free energy estimators (MBAR). Achieving converged and accurate results depends on all of these elements of the calculation. Convergence of the binding free energy is tied to the level of convergence of binding energy distributions at critical intermediate states where bound and unbound states are at equilibrium, and where the rate of binding/unbinding conformational transitions is maximal. This finding mirrors similar observations in the context of order/disorder transitions as for example in protein folding. Insights concerning the physical mechanism of ligand binding and unbinding are obtained. Convergence for the largest FK506 ligand is achieved only after imposing strict conformational restraints, which however require accurate prior structural knowledge of the structure of the complex. The analytical AGBNP2 model is found to underestimate the magnitude of the hydrophobic driving force towards binding in these systems characterized by loosely packed protein-ligand binding interfaces. Rescoring of the binding energies using a numerical surface area model corrects this deficiency. This study illustrates the complex interplay between energy models, exploration of conformational space, and free energy estimators needed to obtain robust estimates from binding free energy calculations.

Elucidating the energetics of entropically driven protein-ligand association: calculations of absolute binding free energy and entropy

The binding of proteins and ligands is generally associated with the loss of translational, rotational, and conformational entropy. In many cases, however, the net entropy change due to binding is positive. To develop a deeper understanding of the energetics of entropically driven protein-ligand binding, we calculated the absolute binding free energies and binding entropies for two HIV-1 protease inhibitors Nelfinavir and Amprenavir using the double-decoupling method with molecular dynamics simulations in explicit solvent. For both ligands, the calculated absolute binding free energies are in general agreement with experiments. The statistical error in the computed ΔG(bind) due to convergence problem is estimated to be ≥2 kcal/mol. The decomposition of free energies indicates that, although the binding of Nelfinavir is driven by nonpolar interaction, Amprenavir binding benefits from both nonpolar and electrostatic interactions. The calculated absolute binding entropies show that (1) Nelfinavir binding is driven by large entropy change and (2) the entropy of Amprenavir binding is much less favorable compared with that of Nelfinavir. Both results are consistent with experiments. To obtain qualitative insights into the entropic effects, we decomposed the absolute binding entropy into different contributions based on the temperature dependence of free energies along different legs of the thermodynamic pathway. The results suggest that the favorable entropic contribution to binding is dominated by the ligand desolvation entropy. The entropy gain due to solvent release from binding site appears to be more than offset by the reduction of rotational and vibrational entropies upon binding.