Plain Molecular Dynamics Research Articles

Computational Analysis of Heat Capacity Effects in Protein-Ligand Binding.

Heat capacity effects in protein-ligand binding as measured by calorimetric experiments have recently attracted considerable attention, particularly in the field of enzyme inhibitor design. A significant negative heat capacity change upon ligand binding implies a marked temperature dependence of the binding enthalpy, which is of high relevance for attempts to optimize protein-ligand interactions. In this work, we address the question of how well such heat capacity changes can be predicted by computer simulations. We examine a series of human thrombin inhibitors that all bind with ΔCp values of about -0.4 kcal/mol/K and calculate heat capacity changes from plain molecular dynamics simulations of the bound and free states of the enzyme and ligand. The results show that accurate ΔCp estimates within a few tenths of a kcal/mol/K of the experimental values can be obtained with this approach. This allows us to address the structural and energetic origin of the negative heat capacity changes for the thrombin inhibitors, and it is found that conformational equilibria of the free ligands in solution make a major contribution to the observed effect.

Molecular Dynamics Simulations in Protein-Protein Docking.

Concerted interactions between all the cell components form the basis of biological processes. Protein-protein interactions (PPIs) constitute a tremendous part of this interaction network. Deeper insight into PPIs can help us better understand numerous diseases and lead to the development of new diagnostic and therapeutic strategies. PPI interfaces, until recently, were considered undruggable. However, it is now believed that the interfaces contain "hot spots," which could be targeted by small molecules. Such a strategy would require high-quality structural data of PPIs, which are difficult to obtain experimentally. Therefore, in silico modeling can complement or be an alternative to in vitro approaches. There are several computational methods for analyzing the structural data of the binding partners and modeling of the protein-protein dimer/oligomer structure. The major problem with in silico structure prediction of protein assemblies is obtaining sufficient sampling of protein dynamics. One of the methods that can take protein flexibility and the effects of the environment into account is Molecular Dynamics (MD). While sampling of the whole protein-protein association process with plain MD would be computationally expensive, there are several strategies to harness the method to PPI studies while maintaining reasonable use of resources. This chapter reviews known applications of MD in the PPI investigation workflows.

Decrypting Allostery in Membrane-Bound K-Ras4B Using Complementary In Silico Approaches Based on Unbiased Molecular Dynamics Simulations.

Protein functions are dynamically regulated by allostery, which enables conformational communication even between faraway residues, and expresses itself in many forms, akin to different "languages": allosteric control pathways predominating in an unperturbed protein are often unintuitively reshaped whenever biochemical perturbations arise (e.g., mutations). To accurately model allostery, unbiased molecular dynamics (MD) simulations require integration with a reliable method able to, e.g., detect incipient allosteric changes or likely perturbation pathways; this is because allostery can operate at longer time scales than those accessible by plain MD. Such methods are typically applied singularly, but we here argue their joint application─as a "multilingual" approach─could work significantly better. We successfully prove this through unbiased MD simulations (∼100 μs) of the widely studied, allosterically active oncotarget K-Ras4B, solvated and embedded in a phospholipid membrane, from which we decrypt allostery using four showcase "languages": Distance Fluctuation analysis and the Shortest Path Map capture allosteric hotspots at equilibrium; Anisotropic Thermal Diffusion and Dynamical Non-Equilibrium MD simulations assess perturbations upon, respectively, either superheating or hydrolyzing the GTP that oncogenically activates K-Ras4B. Chosen "languages" work synergistically, providing an articulate, mutually coherent, experimentally consistent picture of K-Ras4B allostery, whereby distinct traits emerge at equilibrium and upon GTP cleavage. At equilibrium, combined evidence confirms prominent allosteric communication from the membrane-embedded hypervariable region, through a hub comprising helix α5 and sheet β5, and up to the active site, encompassing allosteric "switches" I and II (marginally), and two proposed pockets. Upon GTP cleavage, allosteric perturbations mostly accumulate on the switches and documented interfaces.

Calculation of Heat Capacity Changes in Enzyme Catalysis and Ligand Binding.

It has been suggested that heat capacity changes in enzyme catalysis may be the underlying reason for temperature optima that are not related to unfolding of the enzyme. If this were to be a common phenomenon, it would have major implications for our interpretation of enzyme kinetics. In most cases, the support for the possible existence of a nonzero (negative) activation heat capacity, however, only relies on fitting such a kinetic model to experimental data. It is therefore of fundamental interest to try to use computer simulations to address this issue. One way is simply to calculate the temperature dependence of the activation free energy and determine whether the relationship is linear or not. An alternative approach is to calculate the absolute heat capacities of the reactant and transition states from plain molecular dynamics simulations using either the temperature derivative or fluctuation formula for the enthalpy. Here, we examine these different approaches for a designer enzyme with a temperature optimum that is not caused by unfolding. Benchmark calculations for the heat capacity of liquid water are first carried out using different thermostats. It is shown that the derivative formula for the heat capacity is generally the most robust and insensitive to the thermostat used and its parameters. The enzyme calculations using this method give results in agreement with direct calculations of activation free energies and show no sign of a negative activation heat capacity. We also provide a simple scheme for the calculation of binding heat capacity changes, which is of clear interest in ligand design, and demonstrate it for substrate binding to the designer enzyme. Neither in that case do the simulations predict any negative heat capacity change.

Associative bond swaps in molecular dynamics

We implement a three-body potential to model associative bond swaps, and release it as part of the HOOMD-blue software. The use of a three-body potential to model swaps has been proven to be effective and has recently provided useful insights into the mechanics and dynamics of adaptive network materials such as vitrimers. It is elegant because it can be used in plain molecular dynamics simulations without the need for topology-altering Monte Carlo steps, and naturally represents typical physical features such as slip-bond behavior. It is easily tunable with a single parameter to control the average swap rate. Here, we show how associative bond swaps can be used to speed up the equilibration of systems that self-assemble by avoiding traps and pitfalls, corresponding to long-lived metastable configurations. Our results demonstrate the possibilities of these swaps not only for modeling systems that are associative by nature, but also for increasing simulation efficiency in other systems that are modellable in HOOMD-blue.

Transient Two-Way Molecular-Continuum Coupling with OpenFOAM and MaMiCo: A Sensitivity Study

Molecular-continuum methods, as considered in this work, decompose the computational domain into continuum and molecular dynamics (MD) sub-domains. Compared to plain MD simulations, they greatly reduce computational effort. However, the quality of a fully two-way coupled simulation result strongly depends on a variety of system-specific parameters, and the corresponding sensitivity is only rarely addressed in the literature. Using a state-flux molecular-continuum coupling algorithm, we investigated the influences of various parameters, such as the size of the overlapping region, the coupling time step and the quality of ensemble-based sampling of flow velocities, in a Couette flow scenario. In particular, we considered a big setup in terms of domain size and number of time steps, which allowed us to investigate the long-term behavior of the coupling algorithm close to the incompressible regime. While mostly good agreement was reached on short time scales, it was the long-term behavior which differed even with slightly differently parametrized simulations. We demonstrated our findings by measuring the error in velocity, and we summarize our main observations with a few lessons learned.

Molecular basis for the increased affinity of an RNA recognition motif with re-engineered specificity: A molecular dynamics and enhanced sampling simulations study.

The RNA recognition motif (RRM) is the most common RNA binding domain across eukaryotic proteins. It is therefore of great value to engineer its specificity to target RNAs of arbitrary sequence. This was recently achieved for the RRM in Rbfox protein, where four mutations R118D, E147R, N151S, and E152T were designed to target the precursor to the oncogenic miRNA 21. Here, we used a variety of molecular dynamics-based approaches to predict specific interactions at the binding interface. Overall, we have run approximately 50 microseconds of enhanced sampling and plain molecular dynamics simulations on the engineered complex as well as on the wild-type Rbfox·pre-miRNA 20b from which the mutated systems were designed. Comparison with the available NMR data on the wild type molecules (protein, RNA, and their complex) served to establish the accuracy of the calculations.Free energy calculations suggest that further improvements in affinity and selectivity are achieved by the S151T replacement.

All-atom calculation of protein free-energy profiles.

The Bias Functional (BF) approach is a variational method which enables one to efficiently generate ensembles of reactive trajectories for complex biomolecular transitions, using ordinary computer clusters. For example, this scheme was applied to simulate in atomistic detail the folding of proteins consisting of several hundreds of amino acids and with experimental folding time of several minutes. A drawback of the BF approach is that it produces trajectories which do not satisfy microscopic reversibility. Consequently, this method cannot be used to directly compute equilibrium observables, such as free energy landscapes or equilibrium constants. In this work, we develop a statistical analysis which permits us to compute the potential of mean-force (PMF) along an arbitrary collective coordinate, by exploiting the information contained in the reactive trajectories calculated with the BF approach. We assess the accuracy and computational efficiency of this scheme by comparing its results with the PMF obtained for a small protein by means of plain molecular dynamics.

Extracting Conformational Ensembles of Small Molecules from Molecular Dynamics Simulations: Ampicillin as a Test Case

The accurate and exhaustive description of the conformational ensemble sampled by small molecules in solution, possibly at different physiological conditions, is of primary interest in many fields of medicinal chemistry and computational biology. Recently, we have built an on-line database of compounds with antimicrobial properties, where we provide all-atom force-field parameters and a set of molecular properties, including representative structures extracted from cluster analysis over μs-long molecular dynamics (MD) trajectories. In the present work, we used a medium-sized antibiotic from our sample, namely ampicillin, to assess the quality of the conformational ensemble. To this aim, we compared the conformational landscape extracted from previous unbiased MD simulations to those obtained by means of Replica Exchange MD (REMD) and those originating from three freely-available conformer generation tools widely adopted in computer-aided drug-design. In addition, for different charge/protonation states of ampicillin, we made available force-field parameters and static/dynamic properties derived from both Density Functional Theory and MD calculations. For the specific system investigated here, we found that: (i) the conformational statistics extracted from plain MD simulations is consistent with that obtained from REMD simulations; (ii) overall, our MD-based approach performs slightly better than any of the conformer generator tools if one takes into account both the diversity of the generated conformational set and the ability to reproduce experimentally-determined structures.

Ratcheted molecular-dynamics simulations identify efficiently the transition state of protein folding

The atomistic characterization of the transition state (TS) is a fundamental step to improve the understanding of the folding mechanism and the function of proteins. From a computational point of view, the identification of the conformations that build out the transition state is particularly cumbersome, mainly because of the large computational cost of generating a statistically sound set of folding trajectories. Here we show that a biasing algorithm, based on the physics of the ratchet-and-pawl, can be used to approximate efficiently the transition state. The basic idea is that the algorithmic ratchet exerts a force on the protein when it is climbing the free-energy barrier, while it is inactive when it is descending. The transition state can be identified as the point of the trajectory where the ratchet changes regime. Besides discussing this strategy in general terms, we test it within a protein model whose transition state can be studied independently by plain molecular dynamics simulations. Finally, we show its power in explicit-solvent simulations, obtaining and characterizing a set of transition-state conformations for Acyl-Coenzyme A-Binding Protein (ACBP) and Chymotrypsin Inhibitor 2 (CI2).

Enhanced Conformational Sampling in Molecular Dynamics Simulations of Solvated Peptides: Fragment-Based Local Elevation Umbrella Sampling.

A new method, fragment-based local elevation umbrella sampling (FB-LEUS), is proposed to enhance the conformational sampling in explicit-solvent molecular dynamics (MD) simulations of solvated polymers. The method is derived from the local elevation umbrella sampling (LEUS) method [ Hansen and Hünenberger , J. Comput. Chem. 2010 , 31 , 1 - 23 ], which combines the local elevation (LE) conformational searching and the umbrella sampling (US) conformational sampling approaches into a single scheme. In LEUS, an initial (relatively short) LE build-up (searching) phase is used to construct an optimized (grid-based) biasing potential within a subspace of conformationally relevant degrees of freedom, which is then frozen and used in a (comparatively longer) US sampling phase. This combination dramatically enhances the sampling power of MD simulations but, due to computational and memory costs, is only applicable to relevant subspaces of low dimensionalities. As an attempt to expand the scope of the LEUS approach to solvated polymers with more than a few relevant degrees of freedom, the FB-LEUS scheme involves an US sampling phase that relies on a superposition of low-dimensionality biasing potentials optimized using LEUS at the fragment level. The feasibility of this approach is tested using polyalanine (poly-Ala) and polyvaline (poly-Val) oligopeptides. Two-dimensional biasing potentials are preoptimized at the monopeptide level, and subsequently applied to all dihedral-angle pairs within oligopeptides of 4, 6, 8, or 10 residues. Two types of fragment-based biasing potentials are distinguished: (i) the basin-filling (BF) potentials act so as to "fill" free-energy basins up to a prescribed free-energy level above the global minimum; (ii) the valley-digging (VD) potentials act so as to "dig" valleys between the (four) free-energy minima of the two-dimensional maps, preserving barriers (relative to linearly interpolated free-energy changes) of a prescribed magnitude. The application of these biasing potentials may lead to an impressive enhancement of the searching power (volume of conformational space visited in a given amount of simulation time). However, this increase is largely offset by a deterioration of the statistical efficiency (representativeness of the biased ensemble in terms of the conformational distribution appropriate for the physical ensemble). As a result, it appears difficult to engineer FB-LEUS schemes representing a significant improvement over plain MD, at least for the systems considered here.

Multiscale dynamics of macromolecules using normal mode Langevin.

Proteins and other macromolecules have coupled dynamics over multiple time scales (from femtosecond to millisecond and beyond) that make resolving molecular dynamics challenging. We present an approach based on periodically decomposing the dynamics of a macromolecule into slow and fast modes based on a scalable coarse-grained normal mode analysis. A Langevin equation is used to propagate the slowest degrees of freedom while minimizing the nearly instantaneous degrees of freedom. We present numerical results showing that time steps of up to 1000 fs can be used, with real speedups of up to 200 times over plain molecular dynamics. We present results of successfully folding the Fip35 mutant of WW domain.

Using the local elevation method to construct optimized umbrella sampling potentials: Calculation of the relative free energies and interconversion barriers of glucopyranose ring conformers in water

A method is proposed to combine the local elevation (LE) conformational searching and the umbrella sampling (US) conformational sampling approaches into a single local elevation umbrella sampling (LEUS) scheme for (explicit-solvent) molecular dynamics (MD) simulations. In this approach, an initial (relatively short) LE build-up (searching) phase is used to construct an optimized biasing potential within a subspace of conformationally relevant degrees of freedom, that is then used in a (comparatively longer) US sampling phase. This scheme dramatically enhances (in comparison with plain MD) the sampling power of MD simulations, taking advantage of the fact that the preoptimized biasing potential represents a reasonable approximation to the negative of the free energy surface in the considered conformational subspace. The method is applied to the calculation of the relative free energies of beta-D-glucopyranose ring conformers in water (within the GROMOS 45A4 force field). Different schemes to assign sampled conformational regions to distinct states are also compared. This approach, which bears some analogies with adaptive umbrella sampling and metadynamics (but within a very distinct implementation), is shown to be: (i) efficient (nearly all the computational effort is invested in the actual sampling phase rather than in searching and equilibration); (ii) robust (the method is only weakly sensitive to the details of the build-up protocol, even for relatively short build-up times); (iii) versatile (a LEUS biasing potential database could easily be preoptimized for small molecules and assembled on a fragment basis for larger ones).

Conformation Control of Model Peptides by Metal Ions. A New Type of Turn Structure Found in [(Boc–Cys–Pro–Leu–Cys–Gly–Ala)Hg

Abstract The structure of [(Boc–Cys1–Pro–Leu–Cys4–Gly–Ala–OMe)Hg], 1, in DMF(-d7) and DMSO(-d6) solutions was studied by X-ray absorption fine structure (XAFS), rotating frame nuclear Overhauser effect spectroscopy (ROESY), distance geometry (DG), molecular dynamics (MD), and restrained molecular dynamics (RMD). The XAFS study clarified that the inorganic center of compound 1 adopts a linear coordination with r(HG–S) = 2.33 Å. The NMR experiment revealed 33 atom approximate distances for the 1H-1H pairs of the main chain loop, including the side chains of the cysteinyl residues (ROESY). By a complementary use of DG, plain MD, and RMD for distance information from NMR as well as XAFS, we established that compound 1 adopts a new type of hybridized turn structure. The compound has two hydrogen bonds, Cys1 S-Leu HN and Pro CO-Cys4 HN, among which the former is common to the core sites in proteins, Cysi–X–Y–Cysi+3/M2+ (Zn2+, Fe2+), whereas the latter is common to the mirror image of the γ turn structure in proteins.