Overcoming dissipation in the calculation of standard binding free energies by ligand extraction

This article addresses calculations of the standard free energy of binding from molecular simulations in which a bound ligand is extracted from its binding site by steered molecular dynamics (MD) simulations or equilibrium umbrella sampling (US). Host–guest systems are used as test beds to examine the requirements for obtaining the reversible work of ligand extraction. We find that, for both steered MD and US, marked irreversibilities can occur when the guest molecule crosses an energy barrier and suddenly jumps to a new position, causing dissipation of energy stored in the stretched molecule(s). For flexible molecules, this occurs even when a stiff pulling spring is used, and it is difficult to suppress in calculations where the spring is attached to the molecules by single, fixed attachment points. We, therefore, introduce and test a method, fluctuation‐guided pulling, which adaptively adjusts the spring's attachment points based on the guest's atomic fluctuations relative to the host. This adaptive approach is found to substantially improve the reversibility of both steered MD and US calculations for the present systems. The results are then used to estimate standard binding free energies within a comprehensive framework, termed attach‐pull‐release, which recognizes that the standard free energy of binding must include not only the pulling work itself, but also the work of attaching and then releasing the spring, where the release work includes an accounting of the standard concentration to which the ligand is discharged. © 2013 Wiley Periodicals, Inc.     

Molecular dynamics simulation of the ligand binding domain of mGluR1 in response to agonist and antagonist binding

The interdomain movements of the ligand binding domain (LBD) of mGluR1 in response to agonist or antagonist binding are studied by 2 ns molecular dynamics (MD) simulations. Our results indicate that MD is able to reproduce many of the experimentally determined features of the open and closed conformations of LBD. Analysis of the ligand behavior over time allows to delineate some of the molecular determinants responsible for the agonist-induced or antagonist-blocked LBD responses.     

Energetics and dynamics in MbCN: CN--vibrational relaxation from molecular dynamics simulations

The dynamics of the cyanide anion bound to sperm-whale myoglobin is investigated using atomistic simulations. With density-functional theory, a 2D potential energy surface for the cyanide-heme complex is calculated. Two deep minima with a stabilization energy of approximately 50 kcal/mol corresponding to two different binding orientations (Fe-CN and Fe-NC) of the ligand are found. The Fe-CN conformation is favored over Fe-NC by several kcal/mol. Mixed quantum mechanics/molecular mechanics calculations show that the binding orientation affects the bond strength of the ligand, with a significantly different bond length and a 25 cm-1 shift in the fundamental CN-frequency. For the molecular dynamics (MD) simulations, a 3-center fluctuating charge model for the Fe-CN unit is developed that captures polarization and ligand-metal charge transfer. Stability arguments based on the energetics around the active site and the CN- frequency shifts suggest that the Fe-CN conformation with epsilon-protonation of His epsilon 64 are most likely, which is in agreement with experiment. Both equilibrium and nonequilibrium MD simulations are carried out to investigate the relaxation time scale and possible relaxation pathways in bound MbCN. The nonequilibrium MD simulations with a vibrationally excited ligand reveal that vibrational relaxation takes place on a time scale of hundreds of picoseconds within the active site. This finding supports the hypothesis that the experimentally observed relaxation rate (3.6 ps) reflects the repopulation of the electronic ground state.     

Hydroxamate类抑制剂与MMP-2的分子对接研究

通过分子动力学模拟研究了MMP-2和hydroxamate抑制剂之间的作用模式。在分子动力学模拟中,对于催化区的锌离子和其共价结合的配体(包括抑制剂和组氨酸)采用了键合的模型。从模拟的结果可以看到,R^1取代基和MMP-2的S1疏水口袋中的部分残基能形成很好的几何匹配,从而可以产生很强的范德华和疏水相互作用。模拟结果也表明,两个抑制剂和MMP-2之间分别能形成5个和8个氢键,抑制剂B比A活性更高的原因就是能够形成更加有利氢键作用模式。在整个模拟过程中,催化锌都能保持好的五配位形式,配位键的长度也处于稳定的状态,预测得到的MMP-2和其抑制剂的相互作用模式对于全新抑制剂的设计提供了非常重要的结构信息。     

MDMS: Software Facilitating Performing Molecular Dynamics Simulations

Molecular Dynamics Made Simple (MDMS) is software that facilitates performing molecular dynamics (MD) simulations of solvated protein/protein–ligand complexes with Amber, one of the most popular MD codes. It guides users through the whole process of running MD starting with choosing a protein structure, preparing the model, parametrization of the system, establishing parameters for controlling MD, and finally running simulations. By accommodating every step required for running MD, this software ensures that the simulations performed by a user will provide as realistic insight as it is possible. Its sequential structure and a text-based interface ensure ease of use, while the flexibility required for complex cases is still preserved. MDMS also provides a very time-efficient and streamlined method to start MD simulations, which makes it a feasible tool for both novices and experienced computational chemists. © 2019 Wiley Periodicals, Inc.     

Can MM/GBSA calculations be sped up by system truncation?

We have studied whether calculations of the binding free energy of small ligands to a protein by the MM/GBSA approach (molecular mechanics combined with generalized Born and surface area solvation) can be sped up by including only a restricted number of atoms close to the ligand. If the protein is truncated before the molecular dynamics (MD) simulations, quite large changes are observed for the calculated binding energies, for example, 4 kJ/mol average difference for a radius of 19 Å for the binding of nine phenol derivatives to ferritin. The results are improved if no atoms are fixed in the simulations, with average and maximum errors of 2 and 3 kJ/mol at 19 Å and 3 and 6 kJ/mol at 7 Å. Similar results are obtained for two additional proteins, p38α MAP kinase and factor Xa. On the other hand, if energies are calculated on snapshots that are truncated after the MD simulation, all residues more than 8.5 Å from the ligand can be omitted without changing the energies by more than 1 kJ/mol on average (maximum error 1.4 kJ/mol). At the molecular mechanics level, the gain in computer time for such an approach is small. However, it shows what size of system should be used if the energies instead are calculated with a more demanding method, for example, quantum‐mechanics. © 2017 Wiley Periodicals, Inc.     

Effect of ligand binding on the intraminimum dynamics of proteins

Effects of ligand binding on protein dynamics are studied via molecular dynamics (MD) simulations on two different enzymes, dihydrofolate reductase (DHFR) and triosephosphate isomerase (TIM), in their unliganded (free) and liganded states. Domain motions in MD trajectories are analyzed by collectivities and rotation angles along the principal components (PCs). DHFR in the free state has well‐defined domain rotations, whereas rotations are slightly damped in the binary complex with nicotinamide adenine dinucleotide phosphate (NADPH), and remarkably distorted in the presence of NADP⁺, showing that NADP⁺ is solely responsible for the loss of correlation of the domains in DHFR. Although mean square fluctuations of MD simulations in the same PC subspaces are similar for different ligation states, linear stochastic time series models show that backbone flexibility along the first five PCs is decreased upon NADPH and NADP⁺ binding in subpicosecond scale. This shows that mobility of the protein along the PCs is closely related with intraminimum dynamics, and alterations in ligation states may change the intraminimum dynamics significantly. Low vibrational frequencies of the alpha‐carbon atoms of DHFR are determined from the time series models of a larger number of low indexed PCs, and it is found that number of modes in the lowest frequencies is reduced upon ligand binding. A similar result is obtained for TIM in the unliganded and dihydroxyacetone phosphate bound states. We suggest that stochastic time series modeling is a promising method to be used in determining subtle perturbations in protein dynamics. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Conformational transition and energy landscape of ErbB4 activated by neuregulin1β: one microsecond molecular dynamics simulations

ErbB4, a receptor tyrosine kinase of the ErbB family, plays crucial roles in cell growth and differentiation, especially in the development of the heart and nervous system. Ligand binding to its extracellular region could modulate the activation process. To understand the mechanism of ErbB4 activation induced by ligand binding, we performed one microsecond molecular dynamics (MD) simulations on the ErbB4 extracellular region (ECR) with and without its endogenous ligand neuregulin1β (NRG1β). The conformational transition of the ECR-ErbB4/NRG1β complex from a tethered inactive conformation to an extended active-like form has been observed, while such large and function-related conformational change has not been seen in the simulation on the ECR-ErbB4, suggesting that ligand binding is indeed the active inducing force for the conformational transition and further dimerization. On the basis of MD simulations and principal component analysis, we constructed a rough energy landscape for the conformational transition of ECR-ErbB4/NRG1β complex, suggesting that the conformational change from the inactive state to active-like state involves a stable conformation. The energy barrier for the tether opening was estimated as ~2.7 kcal/mol, which is very close to the experimental value (1-2 kcal/mol) reported for ErbB1. On the basis of the simulation results, an atomic mechanism for the ligand-induced activation of ErbB4 was postulated. The present MD simulations provide a new insight into the conformational changes underlying the activation of ErbB4.     

A molecular dynamics investigation of CDK8/CycC and ligand binding: conformational flexibility and implication in drug discovery

Abnormal activity of cyclin-dependent kinase 8 (CDK8) along with its partner protein cyclin C (CycC) is a common feature of many diseases including colorectal cancer. Using molecular dynamics (MD) simulations, this study determined the dynamics of the CDK8-CycC system and we obtained detailed breakdowns of binding energy contributions for four type-I and five type-II CDK8 inhibitors. We revealed system motions and conformational changes that will affect ligand binding, confirmed the essentialness of CycC for inclusion in future computational studies, and provide guidance in development of CDK8 binders. We employed unbiased all-atom MD simulations for 500 ns on twelve CDK8-CycC systems, including apoproteins and protein–ligand complexes, then performed principal component analysis (PCA) and measured the RMSF of key regions to identify protein dynamics. Binding pocket volume analysis identified conformational changes that accompany ligand binding. Next, H-bond analysis, residue-wise interaction calculations, and MM/PBSA were performed to characterize protein–ligand interactions and find the binding energy. We discovered that CycC is vital for maintaining a proper conformation of CDK8 to facilitate ligand binding and that the system exhibits motion that should be carefully considered in future computational work. Surprisingly, we found that motion of the activation loop did not affect ligand binding. Type-I and type-II ligand binding is driven by van der Waals interactions, but electrostatic energy and entropic penalties affect type-II binding as well. Binding of both ligand types affects protein flexibility. Based on this we provide suggestions for development of tighter-binding CDK8 inhibitors and offer insight that can aid future computational studies.     

Picosecond dynamics of ligand interconversion in the primary docking site of heme proteins

We have used femtosecond IR spectroscopy to probe interconversion dynamics of ligand in the primary docking site of heme proteins under physiological conditions. The docking site, fashioned with highly conserved amino acid residues, modulates ligand-binding activity by mediating the passage of ligand to and from the active binding site. Ligands in two states of the docking site interconvert on the picosecond time scale, and the rates are about 4 times slower in hemoglobin than that in myoglobin. The accurate interconversion rates on the time scale readily accessible by MD simulations can be used to refine computer simulations, which could in turn provide a detailed mechanistic picture of ligand binding in heme proteins.     

Cooperative motions of protein and hydration water molecules: molecular dynamics study of scytalone dehydratase

Two molecular dynamics (MD) simulations totaling 25 ns of simulation time of monomeric scytalone dehydratase (SD) were performed. The enzyme has a ligand-binding pocket containing a cone-shaped alpha+beta barrel, and the C-terminal region covers the binding pocket. Our simulations clarified the difference in protein dynamics and conformation between the liganded protein and the unliganded protein. The liganded protein held the ligand molecule tightly and the initial structure was maintained during the simulation. The unliganded protein, on the other hand, fluctuated dynamically and its structure changed largely from the initial structure. In the equilibrium state, the binding pocket was fully solvated by opening of the C-terminal region, and the protein dynamics was connected with hydration water molecules entry into and release from the binding pocket. In addition, the cooperative motions of the unliganded protein and the hydration water molecules produced the path through the protein interior for ligand binding.     

Computation of binding free energy with molecular dynamics and grand canonical Monte Carlo simulations

The binding of a ligand to a receptor is often associated with the displacement of a number of bound water molecules. When the binding site is exposed to the bulk region, this process may be sampled adequately by standard unbiased molecular dynamics trajectories. However, when the binding site is deeply buried and the exchange of water molecules with the bulk region may be difficult to sample, the convergence and accuracy in free energy perturbation (FEP) calculations can be severely compromised. These problems are further compounded when a reduced system including only the region surrounding the binding site is simulated. To address these issues, we couple molecular dynamics (MD) with grand canonical Monte Carlo (GCMC) simulations to allow the number of water to fluctuate during an alchemical FEP calculation. The atoms in a spherical inner region around the binding pocket are treated explicitly while the influence of the outer region is approximated using the generalized solvent boundary potential (GSBP). At each step during thermodynamic integration, the number of water in the inner region is equilibrated with GCMC and energy data generated with MD is collected. Free energy calculations on camphor binding to a deeply buried pocket in cytochrome P450cam, which causes about seven water molecules to be expelled, are used to test the method. It concluded that solvation free energy calculations with the GCMC/MD method can greatly improve the accuracy of the computed binding free energy compared to simulations with fixed number of water.     

Evaluation of the effect of the chiral centers of Taxol on binding to β-tubulin: A docking and molecular dynamics simulation study

Taxol is one of the most important anti-cancer drugs. The interaction between different variants of Taxol, by altering one of its chiral centers at a time, with β-tubulin protein has been investigated. To achieve such goal, docking and molecular dynamics (MD) simulation studies have been performed. In docking studies, the preferred conformers have been selected to further study by MD method based on the binding energies reported by the AutoDock program. The best result of docking study which shows the highest affinity between ligand and protein has been used as the starting point of the MD simulations. All of the complexes have shown acceptable stability during the simulation process, based on the RMSDs of the backbone of the protein structure. Finally, MM-GBSA calculations have been carried out to select the best ligand, considering the binding energy criteria. The results predict that two of the structures have better affinity toward the mentioned protein, in comparison with Taxol. Three of the structures have affinity similar to that of the Taxol toward the β-tubulin.     

Significant enhancement of docking sensitivity using implicit ligand sampling

The efficient and accurate quantification of protein-ligand interactions using computational methods is still a challenging task. Two factors strongly contribute to the failure of docking methods to predict free energies of binding accurately: the insufficient incorporation of protein flexibility coupled to ligand binding and the neglected dynamics of the protein-ligand complex in current scoring schemes. We have developed a new methodology, named the 'ligand-model' concept, to sample protein conformations that are relevant for binding structurally diverse sets of ligands. In the ligand-model concept, molecular-dynamics (MD) simulations are performed with a virtual ligand, represented by a collection of functional groups that binds to the protein and dynamically changes its shape and properties during the simulation. The ligand model essentially represents a large ensemble of different chemical species binding to the same target protein. Representative protein structures were obtained from the MD simulation, and docking was performed into this ensemble of protein conformation. Similar binding poses were clustered, and the averaged score was utilized to rerank the poses. We demonstrate that the ligand-model approach yields significant improvements in predicting native-like binding poses and quantifying binding affinities compared to static docking and ensemble docking simulations into protein structures generated from an apo MD simulation.     

Theoretical studies on the structural rearrangement of ligand binding pocket in human vitamin D receptor

The structural rearrangement of the ligand binding domain (LBD) of human Vitamin D receptor (hVDR) complexed with 1α, 25‐dihydroxyvitamin D₃ (natural ligand) and its analogues (denoted as b and c ) was studied by molecular dynamics (MD) simulations. MD simulations revealed that these ligands could induce different structural changes of LBD, in which 1α, 25‐dihydroxyvitamin D₃ only led to a minute change, suggesting that LBD adopted its canonical active conformation upon binding the natural ligand, while b and c could provoke a clear structural rearrangement of the LBD. In complex of hVDR‐LBD/ b , it is found that helix 6 (H6) and subsequent loop 6‐7 shift outward and the last turn of H11 shifts away from H12, which generate a new cavity at the bottom of binding pocket to accommodate the extra butyl group on the side chain of ligand b . As for hVDR‐LBD/ c , the steric exclusion of the second side chain of ligand c makes the N‐terminal of H7 move outsides and C‐terminal of H11 close to H12, expanding the bottom of the pocket. These calculation results agree well the experimental observations. © 2010 Wiley Periodicals, Inc. Int J Quantum Chem, 2011     

Oxoanion binding by guanidiniocarbonylpyrrole cations in water: a combined DFT and MD investigation

Structures and properties of nonbonding interactions involving guanidinium-functionalized hosts and carboxylate substrates were investigated by a combination of ab initio and molecular dynamics approaches. The systems under study are on one hand intended to be a model of the arginine-anion bond, so often observed in proteins and nucleic acids, and on the other to provide an opportunity to investigate the influence of molecular structure on the formation of supramolecular complexes in detail. Use of DFT calculations, including extended basis sets and implicit water treatment, allowed us to determine minimum-energy structures and binding enthalpies that compared well with experimental data. Intermolecular forces were found to be mostly due to electrostatic interactions through three hydrogen bonds, one of which is bifurcate, and are sufficiently strong to induce a conformational change in the ligand consisting of a rotation of about 180 degrees around the guanidiniocarbonylpyrrole axis. Free binding energies of the complexes were evaluated through MD simulations performed in the presence of explicit water molecules by use of the molecular mechanics Poisson-Boltzmann solvent accessible surface area (MM-PBSA) and linear interaction energy (LIE) approaches. LIE energies were in quantitative agreement with experimental data. A detailed analysis of the MD simulations revealed that the complexes cannot be described in terms of a single binding structure, but that they are characterized by a significant internal mobility responsible for several low-energy metastable structures.     

Conformational selection or induced fit for Brinker and DNA recognition

Brinker is the key target protein of the Drosophila Decapentaplegic morphogen signalling pathway. Brinker is widely expressed and can bind with DNA. NMR spectra suggest that apo-Brinker is intrinsically unstructured and undergoes a folding transition upon DNA-binding. However, the coupled mechanism of binding and folding is poorly understood. Here, we performed molecular dynamics (MD) simulations for both bound and apo-Brinker to study the mechanism. Room-temperature MD simulations suggest that Brinker becomes more rigid and stable upon DNA-binding. Kinetic analysis of high-temperature MD simulations shows that both bound and apo-Brinker unfold via a two-state process. The time scale of tertiary unfolding is significantly different between bound and apo-Brinker. The predicted Φ-values suggest that there are more residues with native-like transition state ensembles (TSEs) for bound Brinker than for apo-Brinker. The average RMSD differences between bound and apo-Brinker and Kolmogorov-Smirnov (KS) test analysis illustrate that Brinker folding upon DNA-binding might obey induced-fit mechanism based on MD simulations. These methods can be used for the research of other biomolecular folding upon ligand-binding.     

Prediction of protein–protein complexes using replica exchange with repulsive scaling

The realistic prediction of protein–protein complex structures is import to ultimately model the interaction of all proteins in a cell and for the design of new protein–protein interactions. In principle, molecular dynamics (MD) simulations allow one to follow the association process under realistic conditions including full partner flexibility and surrounding solvent. However, due to the many local binding energy minima at the surface of protein partners, MD simulations are frequently trapped for long times in transient association states. We have designed a replica-exchange based scheme employing different levels of a repulsive biasing between partners in each replica simulation. The bias acts only on intermolecular interactions based on an increase in effective pairwise van der Waals radii (repulsive scaling (RS)-REMD) without affecting interactions within each protein or with the solvent. For a set of five protein test cases (out of six) the RS-REMD technique allowed the sampling of near-native complex structures even when starting from the opposide site with respect to the native binding site for one partner. Using the same start structures and same computational demand regular MD simulations sampled near native complex structures only for one case. The method showed also improved results for the refinement of docked structures in the vicinity of the native binding geometry compared to regular MD refinement.