Biased Molecular Dynamics Research Articles

Adaptive enhanced sampling with a path-variable for the simulation of protein folding and aggregation.

In this article, we present a novel adaptive enhanced sampling molecular dynamics (MD) method for the accelerated simulation of protein folding and aggregation. We introduce a path-variable L based on the un-biased momenta p and displacements dq for the definition of the bias s applied to the system and derive 3 algorithms: general adaptive bias MD, adaptive path-sampling, and a hybrid method which combines the first 2 methodologies. Through the analysis of the correlations between the bias and the un-biased gradient in the system, we find that the hybrid methodology leads to an improved force correlation and acceleration in the sampling of the phase space. We apply our method on SPC/E water, where we find a conservation of the average water structure. We then use our method to sample dialanine and the folding of TrpCage, where we find a good agreement with simulation data reported in the literature. Finally, we apply our methodologies on the initial stages of aggregation of a hexamer of Alzheimer's amyloid β fragment 25-35 (Aβ 25-35) and find that transitions within the hexameric aggregate are dominated by entropic barriers, while we speculate that especially the conformation entropy plays a major role in the formation of the fibril as a rate limiting factor.

Dynamic Histogram Analysis To Determine Free Energies and Rates from Biased Simulations.

We present an algorithm to calculate free energies and rates from molecular simulations on biased potential energy surfaces. As input, it uses the accumulated times spent in each state or bin of a histogram and counts of transitions between them. Optimal unbiased equilibrium free energies for each of the states/bins are then obtained by maximizing the likelihood of a master equation (i.e., first-order kinetic rate model). The resulting free energies also determine the optimal rate coefficients for transitions between the states or bins on the biased potentials. Unbiased rates can be estimated, e.g., by imposing a linear free energy condition in the likelihood maximization. The resulting "dynamic histogram analysis method extended to detailed balance" (DHAMed) builds on the DHAM method. It is also closely related to the transition-based reweighting analysis method (TRAM) and the discrete TRAM (dTRAM). However, in the continuous-time formulation of DHAMed, the detailed balance constraints are more easily accounted for, resulting in compact expressions amenable to efficient numerical treatment. DHAMed produces accurate free energies in cases where the common weighted-histogram analysis method (WHAM) for umbrella sampling fails because of slow dynamics within the windows. Even in the limit of completely uncorrelated data, where WHAM is optimal in the maximum-likelihood sense, DHAMed results are nearly indistinguishable. We illustrate DHAMed with applications to ion channel conduction, RNA duplex formation, α-helix folding, and rate calculations from accelerated molecular dynamics. DHAMed can also be used to construct Markov state models from biased or replica-exchange molecular dynamics simulations. By using binless WHAM formulated as a numerical minimization problem, the bias factors for the individual states can be determined efficiently in a preprocessing step and, if needed, optimized globally afterward.

H4 Tails Potentially Produce the Diversity in the Orientation of Two Nucleosomes

Histone tails play an important role in internucleosomal interaction and chromatin compaction. To understand how the H4 tails are involved in the internucleosomal interaction, an adaptively biased molecular dynamics simulation of 63 models of two stacked nucleosomes, each with the H4 tails in different locations, was carried out. This simulation generated a variety of orientations of the separated nucleosomes depending on the formation of the H4 tail bridge between the H4 tails and the DNA of the neighboring nucleosomes. For the models that showed distinctive orientations of the two nucleosomes, the free energies of the separation of the nucleosomes were further investigated using umbrella sampling simulations. The attractive force between the nucleosomes was estimated from the free energies; the force when two H4 tail bridges formed varied from 36 to 63 pN, depending on the formation of the H4 tail-bridge and the interfacial interaction, whereas the force reduced to 15–18 pN after either one of the H4 tail bridges had broken, regardless of the conformation of the H4 tail. Additional simulations of the nucleosomes show that when the H4 tail was truncated, the force between the nucleosomes became repulsive (from−3 to −7 pN). We concluded that the H4 tails potentially produce the diversity in the orientation of the two nucleosomes, which would contribute to the polymorphism of the chromatin structure.

Thermodynamics of camphor migration in cytochrome P450cam by atomistic simulations

Understanding the mechanisms of ligand binding to enzymes is of paramount importance for the design of new drugs. Here, we report on the use of a novel biased molecular dynamics (MD) methodology to study the mechanism of camphor binding to cytochrome P450cam. Microsecond-long MD simulations allowed us to observe reaction coordinates characterizing ligand diffusion from the active site of cytochrome P450cam to solvent via three egress routes. These atomistic simulations were used to estimate thermodynamic quantities along the reaction coordinates and indicate diverse binding configurations. The results suggest that the diffusion of camphor along the pathway near the substrate recognition site (SRS) is thermodynamically preferred. In addition, we show that the diffusion near the SRS is triggered by a transition from a heterogeneous collection of closed ligand-bound conformers to the basin comprising the open conformations of cytochrome P450cam. The conformational change accompanying this switch is characterized by the retraction of the F and G helices and the disorder of the B′ helix. These results are corroborated by experimental studies and provide detailed insight into ligand binding and conformational behavior of the cytochrome family. The presented methodology is general and can be applied to other ligand-protein systems.

Neural Network and Nearest Neighbor Algorithms for Enhancing Sampling of Molecular Dynamics.

The free energy calculations of complex chemical and biological systems with molecular dynamics (MD) are inefficient due to multiple local minima separated by high-energy barriers. The minima can be escaped using an enhanced sampling method such as metadynamics, which apply bias (i.e., importance sampling) along a set of collective variables (CV), but the maximum number of CVs (or dimensions) is severely limited. We propose a high-dimensional bias potential method (NN2B) based on two machine learning algorithms: the nearest neighbor density estimator (NNDE) and the artificial neural network (ANN) for the bias potential approximation. The bias potential is constructed iteratively from short biased MD simulations accounting for correlation among CVs. Our method is capable of achieving ergodic sampling and calculating free energy of polypeptides with up to 8-dimensional bias potential.

Improving the Efficiency of Non-equilibrium Sampling in the Aqueous Environment via Implicit-Solvent Simulations.

By means of estimators based on non-equilibrium work, equilibrium free energy differences or potentials of mean force (PMFs) of a system of interest can be computed from biased molecular dynamics (MD) simulations. The approach, however, is often plagued by slow conformational sampling and poor convergence, especially when the solvent effects are taken into account. Here, as a possible way to alleviate the problem, several widely used implicit-solvent models, which are derived from the analytic generalized Born (GB) equation and implemented in the AMBER suite of programs, were employed in free energy calculations based on non-equilibrium work and evaluated for their abilities to emulate explicit water. As a test case, pulling MD simulations were carried out on an alanine polypeptide with different solvent models and protocols, followed by comparisons of the reconstructed PMF profiles along the unfolding coordinate. The results show that when employing the non-equilibrium work method, sampling with an implicit-solvent model is several times faster and, more importantly, converges more rapidly than that with explicit water due to reduction of dissipation. Among the assessed GB models, the Neck variants outperform the OBC and HCT variants in terms of accuracy, whereas their computational costs are comparable. In addition, for the best-performing models, the impact of the solvent-accessible surface area (SASA) dependent nonpolar solvation term was also examined. The present study highlights the advantages of implicit-solvent models for non-equilibrium sampling.

Impact of Water-Dilution on the Solvation Properties of the Ionic Liquid 1-Methyltriethoxy-3-ethylimidazolium Acetate for Model Biomass Molecules.

Many studies have suggested that the processing of lignocellulosic biomass could provide a renewable feedstock to supplant much of the current demand on petroleum sources. Currently, alkyl imidazolium-based ionic liquids (ILs) have shown considerable promise in the pretreatment, solvation, and hydrolysis of lignocellulosic materials although their high cost and unfavorable viscosity has limited their widespread use. Functionalizing these ILs with an oligo(ethoxy) tail has previously been shown through experiment to decrease the IL's viscosity resulting in enhanced mass transport characteristics, in addition to other favorable traits including decreased inhibition of some enzymes. Additionally, the use of cosolvents to mitigate the cost and unfavorable traits of ILs is an area of growing interest with particular attention on water as the presence of water in biomass processes is inevitable. Through the use of biased and unbiased molecular dynamics (MD) simulations, this study provides a molecular-level perspective of the various solvent-solvent and solvent-solute interactions in binary mixtures of water and 1-methyltriethoxy-3-ethylimidazolium acetate ([Me-(OEt)3-Et-IM+] [OAc-]) in the presence of model cellulose compounds (i.e., glucose and cellobiose). It is observed that at ∼75% w/w IL and water a transition in the nanostructure of the solvent occurs between water-like and IL-like solvation characteristics. It is shown that H-bonding interactions between the anion and water are a major driving force that significantly impacts the solvent properties of the IL as well as conformational preferences of the cellulosic model compound. In addition, it is found that the oligo(ethoxy) cation tail is responsible for the reduction in the propensity for tail aggregation as compared to alkyl tails of similar length, which, combined with increased ionic shielding, results in increased diffusion and enhanced water-like solvation characteristics.

Variational Flooding Study of a SN2 Reaction.

We have studied the reaction dynamics of a prototypical organic reaction using a variationally optimized truncated bias to accelerate transitions between educt and product reactant states. The asymmetric SN2 nucleophilic substitution reaction of fluoromethane and chloromethane CH3F + Cl- ⇌ CH3Cl + F- is considered, and many independent biased molecular dynamics simulations have been performed at 600, 900, and 1200 K, collecting several hundred transitions at each temperature. The transition times and relative rate constants have been obtained for both reaction directions. The activation energies extracted from an Arrhenius plot compare well with standard static calculations.

Subdiffusion in Membrane Permeation of Small Molecules.

Within the solubility–diffusion model of passive membrane permeation of small molecules, translocation of the permeant across the biological membrane is traditionally assumed to obey the Smoluchowski diffusion equation, which is germane for classical diffusion on an inhomogeneous free-energy and diffusivity landscape. This equation, however, cannot accommodate subdiffusive regimes, which have long been recognized in lipid bilayer dynamics, notably in the lateral diffusion of individual lipids. Through extensive biased and unbiased molecular dynamics simulations, we show that one-dimensional translocation of methanol across a pure lipid membrane remains subdiffusive on timescales approaching typical permeation times. Analysis of permeant motion within the lipid bilayer reveals that, in the absence of a net force, the mean squared displacement depends on time as t0.7, in stark contrast with the conventional model, which assumes a strictly linear dependence. We further show that an alternate model using a fractional-derivative generalization of the Smoluchowski equation provides a rigorous framework for describing the motion of the permeant molecule on the pico- to nanosecond timescale. The observed subdiffusive behavior appears to emerge from a crossover between small-scale rattling of the permeant around its present position in the membrane and larger-scale displacements precipitated by the formation of transient voids.

On Dioxygen Permeation of MaL Laccase from the Thermophilic Fungus Melanocarpus albomyces: An all-Atom Molecular Dynamics Investigation.

In this article, biased molecular dynamics (MD) simulations of O2 egress from the active center of MaL laccase toward the bulk solvent were described. Parameterization of the set of four Cu(II) ions, in the framework of CHARMM-36 FF, was carried out on a recent dummy-atom model that takes into account the Jahn-Teller effect. By carrying out a number of statistically relevant MD simulations, under a tiny randomly oriented external force applied to the molecule O2 , three preferred gates for O2 egress from the enzyme and a few intermediate binding pockets (BPs) for O2 were visible; all the gates and pockets were located on two of the three domains. This wide distribution of preferred gates notwithstanding the molecule O2 was seen to follow specific pathways, exploiting consistently the interstices created by the enzyme thermal fluctuations. These are features that can be imagined to have evolved to make MaL laccase extremely efficient as catalysts in various reactions that require O2 .

Thermal Unthreading of the Lasso Peptides Astexin-2 and Astexin-3.

Lasso peptides are a class of knot-like polypeptides in which the C-terminal tail of the peptide threads through a ring formed by an isopeptide bond between the N-terminal amine group and a side chain carboxylic acid. The small size (∼20 amino acids) and simple topology of lasso peptides make them a good model system for studying the unthreading of entangled polypeptides, both with experiments and atomistic simulation. Here, we present an in-depth study of the thermal unthreading behavior of two lasso peptides astexin-2 and astexin-3. Quantitative kinetics and energetics of the unthreading process were determined for variants of these peptides using a series of chromatography and mass spectrometry experiments and biased molecular dynamics (MD) simulations. In addition, we show that the Tyr15Phe variant of astexin-3 unthreads via an unprecedented "tail pulling" mechanism. MD simulations on a model ring-thread system coupled with machine learning approaches also led to the discovery of physicochemical descriptors most important for peptide unthreading.

Product Release Pathways in Human and Plasmodium falciparum Phosphoribosyltransferase.

Atomistic molecular dynamics (MD) simulations coupled with the metadynamics technique were carried out to delineate the product (PPi.2Mg and IMP) release mechanisms from the active site of both human (Hs) and Plasmodium falciparum (Pf) hypoxanthine-guanine-(xanthine) phosphoribosyltransferase (HG(X)PRT). An early movement of PPi.2Mg from its binding site has been observed. The swinging motion of the Asp side chain (D134/D145) in the binding pocket facilitates the detachment of IMP, which triggers the opening of flexible loop II, the gateway to the bulk solvent. In PfHGXPRT, PPi.2Mg and IMP are seen to be released via the same path in all of the biased MD simulations. In HsHGPRT too, the product molecules follow similar routes from the active site; however, an alternate but minor escape route for PPi.2Mg has been observed in the human enzyme. Tyr 104 and Phe 186 in HsHGPRT and Tyr 116 and Phe 197 in PfHGXPRT are the key residues that mediate the release of IMP, whereas the motion of PPi.2Mg away from the reaction center is guided by the negatively charged Asp and Glu and a few positively charged residues (Lys and Arg) that line the product release channels. Mutations of a few key residues present in loop II of Trypanosoma cruzi (Tc) HGPRT have been shown to reduce the catalytic efficiency of the enzyme. Herein, in silico mutation of corresponding residues in loop II of HsHGPRT and PfHGXPRT resulted in partial opening of the flexible loop (loop II), thus exposing the active site to bulk water, which offers a rationale for the reduced catalytic activity of these two mutant enzymes. Investigations of the product release from these HsHGPRT and PfHGXPRT mutants delineate the role of these important residues in the enzymatic turnover.

Structured and Unstructured Binding of an Intrinsically Disordered Protein as Revealed by Atomistic Simulations.

Intrinsically disordered proteins (IDPs) are a set of proteins that lack a definite secondary structure in solution. IDPs can acquire tertiary structure when bound to their partners; therefore, the recognition process must also involve protein folding. The nature of the transition state (TS), structured or unstructured, determines the binding mechanism. The characterization of the TS has become a major challenge for experimental techniques and molecular simulations approaches since diffusion, recognition, and binding is coupled to folding. In this work we present atomistic molecular dynamics (MD) simulations that sample the free energy surface of the coupled folding and binding of the transcription factor c-myb to the cotranscription factor CREB binding protein (CBP). This process has been recently studied and became a model to study IDPs. Despite the plethora of available information, we still do not know how c-myb binds to CBP. We performed a set of atomistic biased MD simulations running a total of 15.6 μs. Our results show that c-myb folds very fast upon binding to CBP with no unique pathway for binding. The process can proceed through both structured or unstructured TS's with similar probabilities. This finding reconciles previous seemingly different experimental results. We also performed Go-type coarse-grained MD of several structured and unstructured models that indicate that coupled folding and binding follows a native contact mechanism. To the best of our knowledge, this is the first atomistic MD simulation that samples the free energy surface of the coupled folding and binding processes of IDPs.

Multiensemble Markov models of molecular thermodynamics and kinetics

We introduce the general transition-based reweighting analysis method (TRAM), a statistically optimal approach to integrate both unbiased and biased molecular dynamics simulations, such as umbrella sampling or replica exchange. TRAM estimates a multiensemble Markov model (MEMM) with full thermodynamic and kinetic information at all ensembles. The approach combines the benefits of Markov state models-clustering of high-dimensional spaces and modeling of complex many-state systems-with those of the multistate Bennett acceptance ratio of exploiting biased or high-temperature ensembles to accelerate rare-event sampling. TRAM does not depend on any rate model in addition to the widely used Markov state model approximation, but uses only fundamental relations such as detailed balance and binless reweighting of configurations between ensembles. Previous methods, including the multistate Bennett acceptance ratio, discrete TRAM, and Markov state models are special cases and can be derived from the TRAM equations. TRAM is demonstrated by efficiently computing MEMMs in cases where other estimators break down, including the full thermodynamics and rare-event kinetics from high-dimensional simulation data of an all-atom protein-ligand binding model.

GPCR Handshake in the Spotlight: Exploring the Dimerization Interface of Dopamine D2 Receptors by Simulations at Multiple Resolutions

Dopamine is an important central nervous system neurotransmitter, involved in the reward-response mechanism of the brain. Dopamine receptors belong to Class-A family of G Protein-Coupled Receptor (GPCR) superfamily and share the characteristic fold of seven transmembrane helices with the rest of the GPCR superfamily. It is known that a member of the dopamine receptor subfamily, Dopamine D2 receptor (D2R), forms dimers via transmembrane helix 4 (TM4). These TM4 dimers have been implicated in diverse pharmacological responses.The key to understanding the behavior of these dimers is the molecular basis of their dimerization. Although previous computational studies have investigated GPCR dimerization, a systematic sampling of conformational space has not been explored to find out the dimer interface in D2R. Our work focuses on this dimerization interface, involving TM4 helices of D2R. To begin with, a homology model derived from D3R crystallographic structure was used to develop a coarse-grained membrane model of TM4 dimers. A systematic biased Molecular Dynamics (MD) study was performed on this system to obtain the free-energy surface of dimerization. Our initial results suggest that D2R forms a TM4-TM4 dimer, with a free-energy profile characteristic of weak GPCR dimers. We are further analyzing the behavior of the D2R dimer, with TM4 interface, using all-atom MD simulations.To verify the possibility of other interfaces, we are investigating different dimer configurations by systematically varying the monomer orientation within them. These interfaces, such as the TM1/H8 interface, have been reported in dimers and higher-order oligomers of other GPCRs. A comparison between them and the TM4 interface will provide useful insights into D2R dimerization.

Computational Study of the Impact of the L483Y Mutation on the Desensitization of the AMPA Receptor

Ionotropic glutamate receptors play a crucial role in fast synaptic neurotransmission in the central nervous system. Because of their implication in many fundamental neurological processes such as memory and learning, these ligand-gated ion channels have been studied with a particular interest. In recent years, crystallography of the individual domains, and more recently a full length AMPA receptor structure, has been extremely useful in interpreting years of mutagenesis and physiological data in terms of functional mechanisms. However, many of the precise details about how these receptors function are still unclear, especially at the level of the single molecule and how they move between different conformational states.Our work is focused on a single point mutation (L483Y) located at the ligand-binding domain dimer interface of the AMPA receptor. The presence of a tyrosine at this particular position is known to prevent desensitization. Despite several crystallographic structures available, the specific interactions stabilizing the dimer interface have not been clearly characterized at the atomistic level. Indeed, inspection of crystal structures suggests that the tyrosine appears to be oriented away from the interface.In order to provide a better understanding of how the L483Y mutation affects the desensitization of the AMPA receptor, we performed a set of classical and biased molecular dynamics simulations using the umbrella sampling method. Firstly, our results suggest that the dimer interface is stabilized by a combination of hydrophobic contacts and hydrogen bond interactions. Secondly, the conformational ensembles obtained using different mutants reveal that only a few sidechains play a critical role in the stabilization of the dimer interface. We will discuss these results with respect to putative desensitization pathways.

Numerical analysis of Pickering emulsion stability: insights from ABMD simulations.

The issue of the stability of Pickering emulsions is tackled at a mesoscopic level using dissipative particle dynamics simulations within the Adiabatic Biased Molecular Dynamics framework. We consider the early stage of the coalescence process between two spherical water droplets in a decane solvent. The droplets are stabilized by Janus nanoparticles of different shapes (spherical and ellipsoidal) with different three-phase contact angles. Given a sufficiently dense layer of particles on the droplets, we show that the stabilization mechanism strongly depends on the collision speed. This is consistent with a coalescence mechanism governed by the rheology of the interfacial region. When the system is forced to coalesce sufficiently slowly, we investigate at a mesoscopic level how the ability of the nanoparticles to stabilize Pickering emulsions is discriminated by nanoparticle mobility and the associated caging effect. These properties are both related to the interparticle interaction and the hydrodynamic resistance in the liquid film between the approaching interfaces.

Binding of Capsaicin to the TRPV1 Ion Channel.

Transient receptor potential (TRP) ion channels constitute a notable family of cation channels involved in the ability of an organisms to detect noxious mechanical, thermal, and chemical stimuli that give rise to the perception of pain, taste, and changes in temperature. One of the most experimentally studied agonist of TRP channels is capsaicin, which is responsible for the burning sensation produced when chili pepper is in contact with organic tissues. Thus, understanding how this molecule interacts and regulates TRP channels is essential to high impact pharmacological applications, particularly those related to pain treatment. The recent publication of a three-dimensional structure of the vanilloid receptor 1 (TRPV1) in the absence and presence of capsaicin from single particle electron cryomicroscopy experiments provides the opportunity to explore these questions at the atomic level. In the present work, molecular docking and unbiased and biased molecular dynamics simulations were employed to generate a structural model of the capsaicin-channel complex. In addition, the standard free energy of binding was estimated using alchemical transformations coupled with conformational, translational, and orientational restraints on the ligand. Key binding modes consistent with previous experimental data are identified, and subtle but essential dynamical features of the binding site are characterized. These observations shed some light into how TRPV1 interacts with capsaicin, and may help to refine design parameters for new TRPV1 antagonists, and potentially guide further developments of TRP channel modulators.

Determination of Protein Folding Intermediate Structures Consistent with Data from Oxidative Footprinting Mass Spectrometry

The mapping of folding landscapes remains an important challenge in protein chemistry. Pulsed oxidative labeling of exposed residues and their detection via mass spectrometry provide new means of taking time-resolved “snapshots” of the structural changes that occur during protein folding. However, such experiments have been so far only interpreted qualitatively. Here, we report the detailed structural interpretation of mass spectrometry data from fast photochemical oxidation of proteins (FPOP) experiments at atomic resolution in a biased molecular dynamics approach. We are able to calculate structures of the early folding intermediate of the model system barstar that are fully consistent with FPOP data and Φ values. Furthermore, structures calculated with both FPOP data and Φ values are significantly less compact and have fewer helical residues than intermediate structures calculated with Φ values only. This improves the agreement with the experimental β-Tanford value and CD measurements. The restraints that we introduce facilitate the structural interpretation of FPOP data and provide new means for refined structure calculations of transiently sampled states on protein folding landscapes.

A robust algorithm for optimizing protein structures with NMR chemical shifts.

Over the past decade, a number of methods have been developed to determine the approximate structure of proteins using minimal NMR experimental information such as chemical shifts alone, sparse NOEs alone or a combination of comparative modeling data and chemical shifts. However, there have been relatively few methods that allow these approximate models to be substantively refined or improved using the available NMR chemical shift data. Here, we present a novel method, called Chemical Shift driven Genetic Algorithm for biased Molecular Dynamics (CS-GAMDy), for the robust optimization of protein structures using experimental NMR chemical shifts. The method incorporates knowledge-based scoring functions and structural information derived from NMR chemical shifts via a unique combination of multi-objective MD biasing, a genetic algorithm, and the widely used XPLOR molecular modelling language. Using this approach, we demonstrate that CS-GAMDy is able to refine and/or fold models that are as much as 10 Å (RMSD) away from the correct structure using only NMR chemical shift data. CS-GAMDy is also able to refine of a wide range of approximate or mildly erroneous protein structures to more closely match the known/correct structure and the known/correct chemical shifts. We believe CS-GAMDy will allow protein models generated by sparse restraint or chemical-shift-only methods to achieve sufficiently high quality to be considered fully refined and "PDB worthy". The CS-GAMDy algorithm is explained in detail and its performance is compared over a range of refinement scenarios with several commonly used protein structure refinement protocols. The program has been designed to be easily installed and easily used and is available at http://www.gamdy.ca.