Optimized explicit-solvent replica exchange molecular dynamics from scratch

Dealing with Induced Fit, Conformational Selection, and Secondary Poses in Molecular Dynamics Simulations for Reliable Free Energy Predictions

In this study, we have tested the performance of standard molecular dynamics (MD) simulations, replicates of shorter standard MD simulations, and Hamiltonian Replica Exchange (HREM) simulations for the sampling of two macrocyclic hosts for guest delivery, characterized by induced fit (phenyl-based host) and conformation selection (naphthyl-based host) and of the ODR-BRD4(I) drug-receptor system where the ligand can assume two main poses. For the optimization of the HREM simulation, we have proposed and tested an on-the-fly iterative scheme for equalizing the acceptance ratio along the replica progression at a constant replica number resulting in a moderate impact of the sampling efficiency. Concerning standard MD, we have found that, while splitting the total allocated simulation time in short MD replicates can reproduce the sampling efficiency of HREM in the phenyl-based host and in the ODR-BRD4(I) complex, in the naphthyl-based macrocycle, characterized by long-lived metastable states, enhanced sampling techniques are the only viable alternative for a reliable canonical sampling of the rugged conformational landscape.

Optimal Replica-Exchange Molecular Simulations in Combination with Evolution Strategies

We have incorporated Evolution Strategies into the Replica-Exchange Monte Carlo simulation method to predict the phase behavior of several example fluids. The replica-exchange method allows one system to exchange temperatures with its neighbors to search for the most stable structure relatively efficiently in a single simulation. However, if the temperature intervals of the replicas are not positioned carefully, there is an issue that local exchange does not occur. Our results for a simple Lennard-Jones fluid and the liquid-crystal Yukawa model demonstrate the utility of the approach when compared to conventional methods. When Evolution Strategies were applied to the Replica-Exchange Monte Carlo simulation, the problem of a significant localized decrease in exchange probability near the phase transition was avoided. By obtaining the optimal temperature intervals, the system efficiently traverses a broader parameter space with a small number of replicas. This is equivalent to accelerating molecular simulations with limited computational resources and can be useful when attempting to predict the phase behavior of complex systems.

Computational methods for exploring protein conformations

Proteins are dynamic molecules that can transition between a potentially wide range of structures comprising their conformational ensemble. The nature of these conformations and their relative probabilities are described by a high-dimensional free energy landscape. While computer simulation techniques such as molecular dynamics simulations allow characterisation of the metastable conformational states and the transitions between them, and thus free energy landscapes, to be characterised, the barriers between states can be high, precluding efficient sampling without substantial computational resources. Over the past decades, a dizzying array of methods have emerged for enhancing conformational sampling, and for projecting the free energy landscape onto a reduced set of dimensions that allow conformational states to be distinguished, known as collective variables (CVs), along which sampling may be directed. Here, a brief description of what biomolecular simulation entails is followed by a more detailed exposition of the nature of CVs and methods for determining these, and, lastly, an overview of the myriad different approaches for enhancing conformational sampling, most of which rely upon CVs, including new advances in both CV determination and conformational sampling due to machine learning.

Pristine and Hydroxylated Fullerenes Prevent the Aggregation of Human Islet Amyloid Polypeptide and Display Different Inhibitory Mechanisms

Protein aggregation, involving the formation of dimers, oligomers, and fibrils, is associated with many human diseases. Type 2 diabetes is one of the common amyloidosis and linked with the aggregation of human islet amyloid polypeptide (hIAPP). A series of nanoparticles are reported to be able to interact with proteins and enhance/inhibit protein aggregation. However, the effects of C₆₀ (a model system of hydrophobic nanoparticle) and C₆₀(OH)₈ (a hydroxylated fullerene) on hIAPP aggregation remain unknown. In this study, we investigate the influences of pristine fullerene C₆₀ and hydroxylated C₆₀ on the dimerization of hIAPP using molecular dynamics (MD) simulations. Extensive replica exchange molecular dynamics (REMD) simulations show that isolated hIAPP dimers adopt β-sheet structure containing the amyloid-precursor (β-hairpin). Both C₆₀ and C₆₀(OH)₈ notably inhibit the β-sheet formation of hIAPP dimer and induce the formation of collapsed disordered coil-rich conformations. Protein—nanoparticle interaction analyses reveal that the inhibition of hIAPP aggregation by C₆₀ is mainly via hydrophobic and aromatic-stacking interactions, while the prevention of hIAPP aggregation by C₆₀(OH)₈ is mostly through collective hydrogen bonding and aromatic-stacking interactions. Conventional MD simulations indicate that both C₆₀ and C₆₀(OH)₈ weaken the interactions within hIAPP protofibril and disrupt the β-sheet structure. These results provide mechanistic insights into the possible inhibitory mechanism of C₆₀ and C₆₀(OH)₈ toward hIAPP aggregation, and they are of great reference value for the screening of potent amyloid inhibitors.

Enhanced Molecular Dynamics Simulations of Intrinsically Disordered Proteins

Molecular dynamics simulations represent a powerful tool to gain insights into structural and dynamical features of biomolecular systems. Nevertheless, their recognized limitation in terms of achievable timescales becomes particularly severe when dealing with slow processes. In such cases, the employment of enhanced sampling methods, which allow accelerating the characterization of rare events in a timeframe consistent with conventional computational resources, results as crucial. In particular, such advanced techniques have proven highly valuable in the context of protein folding and, specifically, to explore the conformational ensemble spanned by intrinsically disordered proteins (IDPs). Here, we describe how to set up molecular dynamics simulations with one of these enhanced sampling approaches (namely, Parallel Tempering Metadynamics in the Well-Tempered Ensemble) using the N_TAIL peptide as a test case.

Binding of Copper Ions with Octapeptide Region in Prion Protein: Simulations with Charge Transfer Model

Copper ions are important cofactors of many metalloproteins. The binding dynamics of proteins to the copper ion is important for biological functions but is less understood at the microscopic level. What are the key factors determining the recognition and the stabilization of the copper ion during the binding? Our work investigates the binding dynamics of the copper ion with a simple system (the N-terminus of PrP) using simulation methods. To precisely characterize the protein?ion interaction, we build up an effective copper?peptide force field based on quantum chemistry calculations. In our model, the effects of charge transfer, protonation/deprotonation, and induced polarization are considered. With this force field, we successfully characterize the local structures and the complex interactions of the octapeptide around the copper ion. Furthermore, using an enhanced sampling method, the binding/unbinding processes of the copper ion with the octapeptide are simulated. Free-energy landscapes are generated in consequence, and multiple binding pathways are characterized. It is observed that various native ligands contribute differently to the binding processes. Some residues are related to the capture of the ion (behaving like ?arm?s), and some others contribute to the stabilization of the coordination structure (acting like ?core?s). These different interactions induce various pathways. Besides, a nonnative binding ligand is determined, and it has essential contributions and modulations to the binding pathways. With all these results, the picture of copper?octapeptide binding is outlined. These features are believed to happen in many ion?peptide interactions, such as the cooperative stabilization between the coordinations with neighboring backbone nitrogens and an auxiliary intermediate coordination with the neighboring oxygen from the N-terminal direction. We believe that our studies are valuable to understand the complicated ion?peptide binding processes.

Introduction of a bounded penalty function in contact-assisted simulations of protein structures to omit false restraints

Contact-assisted simulations, the contacts being predicted or determined experimentally, have become very important in the determination of the structures of proteins and other biological macromolecules. In this work, the effect of contact-distance restraints on the simulated structures was investigated with the use of multiplexed replica exchange simulations with the coarse-grained UNRES force field. A modified bounded flat-bottom restraint function that does not generate a gradient when a restraint cannot be satisfied was implemented. Calculations were run with (i) a set of four small proteins, with contact restraints derived from experimental structures, and (ii) selected CASP11 and CASP12 targets, with restraints as used at prediction time. The bounded penalty function largely omitted false contacts, which were usually inconsistent. It was found that at least 20% of correct contacts must be present in the restraint set to improve model quality with respect to unrestrained simulations. © 2019 Wiley Periodicals, Inc.

Alchemical Free-Energy Calculations by Multiple-Replica λ-Dynamics: The Conveyor Belt Thermodynamic Integration Scheme

A new method is proposed to calculate alchemical free-energy differences based on molecular dynamics (MD) simulations, called the conveyor belt thermodynamic integration (CBTI) scheme. As in thermodynamic integration (TI), K replicas of the system are simulated at different values of the alchemical coupling parameter λ. The number K is taken to be even, and the replicas are equally spaced on a forward-turn-backward-turn path, akin to a conveyor belt (CB) between the two physical end-states; and as in λ-dynamics (λD), the λ-values associated with the individual systems evolve in time along the simulation. However, they do so in a concerted fashion, determined by the evolution of a single dynamical variable Λ of period 2π controlling the advance of the entire CB. Thus, a change of Λ is always associated with K/2 equispaced replicas moving forward and K/2 equispaced replicas moving backward along λ. As a result, the effective free-energy profile of the replica system along Λ is periodic of period 2 πK^-1, and the magnitude of its variations decreases rapidly upon increasing K, at least as K^-1 in the limit of large K. When a sufficient number of replicas is used, these variations become small, which enables a complete and quasi-homogeneous coverage of the λ-range by the replica system, without application of any biasing potential. If desired, a memory-based biasing potential can still be added to further homogenize the sampling, the preoptimization of which is computationally inexpensive. The final free-energy profile along λ is calculated similarly to TI, by binning of the Hamiltonian λ-derivative as a function of λ considering all replicas simultaneously, followed by quadrature integration. The associated quadrature error can be kept very low owing to the continuous and quasi-homogeneous λ-sampling. The CBTI scheme can be viewed as a continuous/deterministic/dynamical analog of the Hamiltonian replica-exchange/permutation (HRE/HRP) schemes or as a correlated multiple-replica analog of the λD or λ-local elevation umbrella sampling (λ-LEUS) schemes. Compared to TI, it shares the advantage of the latter schemes in terms of enhanced orthogonal sampling, i.e. the availability of variable-λ paths to circumvent conformational barriers present at specific λ-values. Compared to HRE/HRP, it permits a deterministic and continuous sampling of the λ-range, is expected to be less sensitive to possible artifacts of the thermo- and barostating schemes, and bypasses the need to carefully preselect a λ-ladder and a swapping-attempt frequency. Compared to λ-LEUS, it eliminates (or drastically reduces) the dead time associated with the preoptimization of a biasing potential. The goal of this article is to provide the mathematical/physical formulation of the proposed CBTI scheme, along with an initial application of the method to the calculation of the hydration free energy of methanol.

Selective enhanced sampling in dihedral energy facilitates overcoming the dihedral energy increase in protein folding and accelerates the searching for protein native structure

A dihedral-energy-based selective enhanced sampling method (D-SITSMD) is presented with improved capabilities for searching a protein's natively folded structure and for providing the underlying folding pathway.

Influence of fullerenol on hIAPP aggregation: amyloid inhibition and mechanistic aspects

C60(OH)₂₄inhibits hIAPP aggregation by suppressing the fibril-prone structure and destabilizes hIAPP protofibrils by binding to the amyloid core region.

The Inhibitory Effect of Hydroxylated Carbon Nanotubes on the Aggregation of Human Islet Amyloid Polypeptide Revealed by a Combined Computational and Experimental Study

Fibrillar deposits formed by the aggregation of the human islet amyloid polypeptide (hIAPP) are the major pathological hallmark of type 2 diabetes mellitus (T2DM). Inhibiting the aggregation of hIAPP is considered the primary therapeutic strategy for the treatment of T2DM. Hydroxylated carbon nanoparticles have received great attention in impeding amyloid protein fibrillation owing to their reduced cytotoxicity compared to the pristine ones. In this study, we investigated the influence of hydroxylated single-walled carbon nanotubes (SWCNT-OHs) on the first step of hIAPP aggregation: dimerization by performing explicit solvent replica exchange molecular dynamics (REMD) simulations. Extensive REMD simulations demonstrate that SWCNT-OHs can dramatically inhibit interpeptide β-sheet formation and completely suppress the previously reported β-hairpin amyloidogenic precursor of hIAPP. On the basis of our simulation results, we proposed that SWCNT-OH can hinder hIAPP fibrillation. This was further confirmed by our systematic turbidity measurements, thioflavin T fluorescence, circular dichroism (CD), transmission electron microscope (TEM), and atomic force microscopy (AFM) experiments. Detailed analyses of hIAPP-SWCNT-OH interactions reveal that hydrogen bonding, van der Waals, and π-stacking interactions between hIAPP and SWCNT-OH significantly weaken the inter- and intrapeptide interactions that are crucial for β-sheet formation. Our collective computational and experimental data reveal not only the inhibitory effect but also the inhibitory mechanism of SWCNT-OH against hIAPP aggregation, thus providing new clues for the development of future drug candidates against T2DM.

Temperature-induced collapse of a disordered peptide observed by three sampling methods in molecular dynamics simulations

The conformational ensembles of a disordered peptide, polyglutamine Q15, over a wide temperature range were sampled using multiple replicates of conventional molecular dynamics (cMD) simulations as well as two enhanced sampling methods, temperature replica exchange (TREMD) and replica exchange with solute tempering (REST). The radius of gyration, asphericity, secondary structure, and hydrogen bonding patterns were used for the comparison of the sampling methods. Overall, the three sampling methods generated similar conformational ensembles, with progressive collapse at higher temperatures. Although accumulating the longest simulation time (90 μs), cMD at room temperature missed a small subspace that was sampled by both TREMD and REST. This subspace was high in α-helical content and separated from the main conformational space by an energy barrier. REST used less simulation time than TREMD (36 μs versus 42 μs), and this gap is expected to widen significantly for larger disordered proteins. We conclude that REST is the method of choice for conformational sampling of intrinsically disordered proteins.

A generalized linear response framework for expanded ensemble and replica exchange simulations

Expanded ensemble simulation is a powerful technique for enhancing sampling over a range of thermodynamic parameters. However, although the premise is relatively simple, running successful simulations in practice still presents something of an ad hoc challenge. Three main difficulties exist: (1) the selection of the thermodynamic states, (2) the selection of the sampling weights, and (3) efficient sampling of the expanded parameter space. Here we consider these problems in the context of a pairwise linear response approach to the work fluctuation theorem. The approach offers comprehensive tactics for addressing the three difficulties and can be used as either an alternative or a complement to replica exchange simulations. Importantly, the results are trivially implemented for multi-dimensional parameter spaces and they recover results from the literature aimed at the special cases of simulated/parallel tempering and replica exchange umbrella sampling. Illustrative examples are shown using the NAMD simulation engine.

Online Optimization of Total Acceptance in Hamiltonian Replica Exchange Simulations

Replica exchange is a widely used sampling strategy in molecular simulation. While a variety of methods exist to optimize parameters for temperature replica exchange, less is known about how to optimize parameters for more general Hamiltonian replica exchange simulations. We present an algorithm for the online optimization of total acceptance for both temperature and Hamiltonian replica exchange simulations using stochastic gradient descent. We optimize the total acceptance, a heuristic objective function capturing the efficiency of replica exchange. Our approach is general and has several desirable properties, including: (1) it makes few assumptions about the system of interest, (2) optimization occurs online without the requirement of presimulation, and (3) most importantly, it readily generalizes to systems where there are multiple control parameters (e.g., temperatures, force constants, etc.) that determine the Hamiltonian of each replica. We explore some general properties of the algorithm on a simple harmonic oscillator system, and demonstrate its effectiveness on a more complex data-guided protein folding simulation.

Assemblies of amyloid-β30-36 hexamer and its G33V/L34T mutants by replica-exchange molecular dynamics simulation

The aggregation of amyloid-β peptides is associated with the pathogenesis of Alzheimer’s disease, in which the 30–36 fragments play an important part as a fiber-forming hydrophobic region. The fibrillar structure of Aβ_30–36 has been detected by means of X-ray diffraction, but its oligomeric structural determination, biophysical characterization, and pathological mechanism remain elusive. In this study, we have investigated the structures of Aβ_30–36 hexamer as well as its G33V and L34T mutants in explicit water environment using replica-exchange molecular dynamics (REMD) simulations. Our results show that the wild-type (WT) Aβ_30–36 hexamer has a preference to form β-barrel and bilayer β-sheet conformations, while the G33V or L34T mutation disrupts the β-barrel structures: the G33V mutant is homogenized to adopt β-sheet-rich bilayers, and the structures of L34T mutant on the contrary get more diverse. The hydrophobic interaction plays a critical role in the formation and stability of oligomeric assemblies among all the three systems. In addition, the substitution of G33 by V reduces the β-sheet content in the most populated conformations of Aβ_30–36 oligomers through a steric effect. The L34T mutation disturbs the interpeptide hydrogen bonding network, and results in the increased coil content and morphological diversity. Our REMD runs provide structural details of WT and G33V/L34T mutant Aβ_30–36 oligomers, and molecular insight into the aggregation mechanism, which will be helpful for designing novel inhibitors or amyloid-based materials.

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.

Comparing generalized ensemble methods for sampling of systems with many degrees of freedom

We compare two standard replica exchange methods using temperature and dielectric constant as the scaling variables for independent replicas against two new corresponding enhanced sampling methods based on non-equilibrium statistical cooling (temperature) or descreening (dielectric). We test the four methods on a rough 1D potential as well as for alanine dipeptide in water, for which their relatively small phase space allows for the ability to define quantitative convergence metrics. We show that both dielectric methods are inferior to the temperature enhanced sampling methods, and in turn show that temperature cool walking (TCW) systematically outperforms the standard temperature replica exchange (TREx) method. We extend our comparisons of the TCW and TREx methods to the 5 residue met-enkephalin peptide, in which we evaluate the Kullback-Leibler divergence metric to show that the rate of convergence between two independent trajectories is faster for TCW compared to TREx. Finally we apply the temperature methods to the 42 residue amyloid-β peptide in which we find non-negligible differences in the disordered ensemble using TCW compared to the standard TREx. All four methods have been made available as software through the OpenMM Omnia software consortium (http://www.omnia.md/).

Enhanced conformational sampling using replica exchange with concurrent solute scaling and hamiltonian biasing realized in one dimension

Replica exchange (REX) is a powerful computational tool for overcoming the quasi-ergodic sampling problem of complex molecular systems. Recently, several multidimensional extensions of this method have been developed to realize exchanges in both temperature and biasing potential space or the use of multiple biasing potentials to improve sampling efficiency. However, increased computational cost due to the multidimensionality of exchanges becomes challenging for use on complex systems under explicit solvent conditions. In this study, we develop a one-dimensional (1D) REX algorithm to concurrently combine the advantages of overall enhanced sampling from Hamiltonian solute scaling and the specific enhancement of collective variables using Hamiltonian biasing potentials. In the present Hamiltonian replica exchange method, termed HREST-BP, Hamiltonian solute scaling is applied to the solute subsystem, and its interactions with the environment to enhance overall conformational transitions and biasing potentials are added along selected collective variables associated with specific conformational transitions, thereby balancing the sampling of different hierarchical degrees of freedom. The two enhanced sampling approaches are implemented concurrently allowing for the use of a small number of replicas (e.g., 6 to 8) in 1D, thus greatly reducing the computational cost in complex system simulations. The present method is applied to conformational sampling of two nitrogen-linked glycans (N-glycans) found on the HIV gp120 envelope protein. Considering the general importance of the conformational sampling problem, HREST-BP represents an efficient procedure for the study of complex saccharides, and, more generally, the method is anticipated to be of general utility for the conformational sampling in a wide range of macromolecular systems.

Critical Nucleus Structure and Aggregation Mechanism of the C-terminal Fragment of Copper-Zinc Superoxide Dismutase Protein

The aggregation of the copper-zinc superoxide dismutase (SOD1) protein is linked to familial amyotrophic lateral sclerosis, a progressive neurodegenerative disease. A recent experimental study has shown that the (147)GVIGIAQ(153) SOD1 C-terminal segment not only forms amyloid fibrils in isolation but also accelerates the aggregation of full-length SOD1, while substitution of isoleucine at site 149 by proline blocks its fibril formation. Amyloid formation is a nucleation-polymerization process. In this study, we investigated the oligomerization and the nucleus structure of this heptapeptide. By performing extensive replica-exchange molecular dynamics (REMD) simulations and conventional MD simulations, we found that the GVIGIAQ hexamers can adopt highly ordered bilayer β-sheets and β-barrels. In contrast, substitution of I149 by proline significantly reduces the β-sheet probability and results in the disappearance of bilayer β-sheet structures and the increase of disordered hexamers. We identified mixed parallel-antiparallel bilayer β-sheets in both REMD and conventional MD simulations and provided the conformational transition from the experimentally observed parallel bilayer sheets to the mixed parallel-antiparallel bilayer β-sheets. Our simulations suggest that the critical nucleus consists of six peptide chains and two additional peptide chains strongly stabilize this critical nucleus. The stabilized octamer is able to recruit additional random peptides into the β-sheet. Therefore, our simulations provide insights into the critical nucleus formation and the smallest stable nucleus of the (147)GVIGIAQ(153) peptide.

The effects of organic solvents on the folding pathway and associated thermodynamics of proteins: a microscopic view

Protein folding is subject to the effects of solvation environment. A variety of organic solvents are used as additives for in vitro refolding of denatured proteins. Examination of the solvent effects on protein folding could be of fundamental importance to understand the molecular interactions in determining protein structure. This article investigated the folding of α-helix and β-hairpin structures in water and the solutions of two representative refolding additives (methanol (MeOH) and 1-Ethyl-3-methylimidazolium chloride (EMIM-Cl) ionic liquid) using REMD simulations. For both α-helix and β-hairpin in MeOH/water solution or α-helix in EMIM-Cl/water solution, the transient structures along the folding pathway are consistent with the counterparts in water but the relative statistical weights are changed, leading to the decrease in the overall folding free energy barrier. Accordingly, MeOH promotes the folding of both α-helix and β-hairpin but EMIM-Cl ionic liquid only promotes the folding of α-helix, consistent with experimental observations. The present study reveals for the first time the trivial effects on folding route but significant effects on folding thermodynamics from MeOH and EMIM-Cl, explaining the function of protein refolding additives and testifying the validity of the folding mechanism revealed by in vitro protein folding study using refolding additives.

Simulated Solute Tempering

For the enhanced conformational sampling in molecular dynamics (MD) simulations, we present "simulated solute tempering" (SST) which is an easy to implement variant of simulated tempering. SST extends conventional simulated tempering (CST) by key concepts of "replica exchange with solute tempering" (REST, Liu et al. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 13749). We have applied SST, CST, and REST to molecular dynamics (MD) simulations of an alanine octapeptide in explicit water. The weight parameters required for CST and SST are determined by two different formulas whose performance is compared. For SST only one of them yields a uniform sampling of the temperature space. Compared to CST and REST, SST provides the highest exchange probabilities between neighboring rungs in the temperature ladder. Concomitantly, SST leads to the fastest diffusion of the simulation system through the temperature space, in particular, if the "even-odd" exchange scheme is employed in SST. As a result, SST exhibits the highest sampling speed of the investigated tempering methods.

Amphiphilic Peptides A6K and V6K Display Distinct Oligomeric Structures and Self-Assembly Dynamics: A Combined All-Atom and Coarse-Grained Simulation Study

Amphiphilic peptides can self-assemble into ordered nanostructures with different morphologies. However, the assembly mechanism and the structures of the early assemblies prior to nanostructure formation remain elusive. In this study, we investigated the oligomeric structures of two amphiphilic heptapeptides A6K and V6K by all-atom explicit-solvent replica-exchange molecular dynamics (REMD) simulations, and then examined the assembly dynamics of large aggregates by coarse-grained (CG) MD simulations. Our 200 ns REMD simulations show that A6K peptides predominantly adopt loosely packed disordered coil aggregates, whereas V6K peptides mostly assemble into compact β-sheet-rich conformations, consistent with the signal measured experimentally in aqueous solution. Well-organized β-sheet-rich conformations, albeit with low population, are also populated for V6K octamers, including bilayer β-sheets and β-barrels. These ordered β-sheet-rich conformations are observed for the first time for amphiphilic peptides. Our 10-μs CG-MD simulations on 200 peptide chains demonstrate that A6K and V6K peptides follow two different self-assembly processes, and the former form monolayer lamellas while the latter assemble into plate-like assemblies. CG-MD simulations also show that V6K peptides display higher assembly capability than A6K, in support of our all-atom REMD simulation results. Interpeptide interaction analyses reveal that the marked differences in oligomeric structures and assembly dynamics between A6K and V6K result from the subtle interplay of competition among hydrophobic, hydrogen-bonding, and electrostatic interactions of the two peptides. Our study provides structural and mechanistic insights into the initial self-assembly process of A6K and V6K at the molecular level.

Exploring biomolecular dynamics and interactions using advanced sampling methods

Molecular dynamics (MD) and Monte Carlo (MC) simulations have emerged as a valuable tool to investigate statistical mechanics and kinetics of biomolecules and synthetic soft matter materials. However, major limitations for routine applications are due to the accuracy of the molecular mechanics force field and due to the maximum simulation time that can be achieved in current simulations studies. For improving the sampling a number of advanced sampling approaches have been designed in recent years. In particular, variants of the parallel tempering replica-exchange methodology are widely used in many simulation studies. Recent methodological advancements and a discussion of specific aims and advantages are given. This includes improved free energy simulation approaches and conformational search applications.

Superposition-Enhanced Estimation of Optimal Temperature Spacings for Parallel Tempering Simulations

Effective parallel tempering simulations rely crucially on a properly chosen sequence of temperatures. While it is desirable to achieve a uniform exchange acceptance rate across neighboring replicas, finding a set of temperatures that achieves this end is often a difficult task, in particular for systems undergoing phase transitions. Here we present a method for determination of optimal replica spacings, which is based upon knowledge of local minima in the potential energy landscape. Working within the harmonic superposition approximation, we derive an analytic expression for the parallel tempering acceptance rate as a function of the replica temperatures. For a particular system and a given database of minima, we show how this expression can be used to determine optimal temperatures that achieve a desired uniform acceptance rate. We test our strategy for two atomic clusters that exhibit broken ergodicity, demonstrating that our method achieves uniform acceptance as well as significant efficiency gains.

Aβ(16-22) peptides can assemble into ordered β-barrels and bilayer β-sheets, while substitution of phenylalanine 19 by tryptophan increases the population of disordered aggregates

A recent experimental study reported that termini-uncapped Aβ(16-22) (with sequence KLVFFAE) peptides self-assembled into nanofibrils at pH 2.0. The oligomerization of this uncapped peptide at atomic level in acidic pH condition remains to be determined, as computational studies mainly focus on the self-assembly of capped Aβ(16-22) peptides at neutral pH condition. In this study, using replica exchange molecular dynamics (REMD) simulations with explicit solvent, we investigated the octameric structures of the uncapped Aβ(16-22) and its F19W variant at acidic pH condition. Our simulations reveal that the Aβ(16-22) octamers adopt various conformations, including closed β-barrels, bilayer β-sheets, and disordered aggregates. The closed β-barrel conformation is particularly interesting, as the cylindrical β-barrel has been reported recently as a cytotoxic species. Interpeptide contact probability analyses between all pairs of residues reveal that the hydrophobic and aromatic stacking interactions between F19 residues play an essential role in the formation of β-barrels and bilayer β-sheets. The importance of F19 and the steric effect on the structures of Aβ(16-22) octamers are further examined by REMD simulation of F19W mutant. This REMD run shows that substitution of F19 by W with a more bulky aromatic side chain significantly reduces the β-sheet content and in turn enhances the population of disordered aggregates, indicating that the steric effect significantly affect the self-assembly of low molecular weight Aβ(16-22) oligomers.

Advanced replica-exchange sampling to study the flexibility and plasticity of peptides and proteins

Molecular dynamics (MD) simulations are ideally suited to investigate protein and peptide plasticity and flexibility simultaneously at high spatial (atomic) and high time resolution. However, the applicability is still limited by the force field accuracy and by the maximum simulation time that can be routinely achieved in current MD simulations. In order to improve the sampling the replica-exchange (REMD) methodology has become popular and is now the most widely applied advanced sampling approach. Many variants of the REMD method have been designed to reduce the computational demand or to enhance sampling along specific sets of conformational variables. An overview on recent methodological advances and discussion of specific aims and advantages of the approaches will be given. Applications in the area of free energy simulations and advanced sampling of intrinsically disordered peptides and proteins will also be discussed. This article is part of a Special Issue entitled: The emerging dynamic view of proteins: Protein plasticity in allostery, evolution and self-assembly.

Replica exchanging self-guided Langevin dynamics for efficient and accurate conformational sampling

This work presents a replica exchanging self-guided Langevin dynamics (RXSGLD) simulation method for efficient conformational searching and sampling. Unlike temperature-based replica exchanging simulations, which use high temperatures to accelerate conformational motion, this method uses self-guided Langevin dynamics (SGLD) to enhance conformational searching without the need to elevate temperatures. A RXSGLD simulation includes a series of SGLD simulations, with simulation conditions differing in the guiding effect and/or temperature. These simulation conditions are called stages and the base stage is one with no guiding effect. Replicas of a simulation system are simulated at the stages and are exchanged according to the replica exchanging probability derived from the SGLD partition function. Because SGLD causes less perturbation on conformational distribution than high temperatures, exchanges between SGLD stages have much higher probabilities than those between different temperatures. Therefore, RXSGLD simulations have higher conformational searching ability than temperature based replica exchange simulations. Through three example systems, we demonstrate that RXSGLD can generate target canonical ensemble distribution at the base stage and achieve accelerated conformational searching. Especially for large systems, RXSGLD has remarkable advantages in terms of replica exchange efficiency, conformational searching ability, and system size extensiveness.

Transmembrane helix assembly by window exchange umbrella sampling

A method of window exchange umbrella sampling molecular dynamics simulation is employed for transmembrane helix assembly. An analytical expression for the average acceptance probability between neighboring windows is derived and combined with the first passage time optimization method to predetermine a parameter set in an optimal range. With the parameter set, the method provides a substantially more efficient sampling of helix-helix interfaces together with the potential of mean force along the helix-helix distance of a transmembrane helix-dimer model, compared to the umbrella sampling method.

Continuous Constant pH Molecular Dynamics in Explicit Solvent with pH-Based Replica Exchange

A computational tool that offers accurate pKa values and atomically detailed knowledge of protonation-coupled conformational dynamics is valuable for elucidating mechanisms of energy transduction processes in biology, such as enzyme catalysis and electron transfer as well as proton and drug transport. Toward this goal we present a new technique of embedding continuous constant pH molecular dynamics within an explicit-solvent representation. In this technique we make use of the efficiency of the generalized-Born (GB) implicit-solvent model for estimating the free energy of protein solvation while propagating conformational dynamics using the more accurate explicit-solvent model. Also, we employ a pH-based replica exchange scheme to significantly enhance both protonation and conformational state sampling. Benchmark data of five proteins including HP36, NTL9, BBL, HEWL, and SNase yield an average absolute deviation of 0.53 and a root mean squared deviation of 0.74 from experimental data. This level of accuracy is obtained with 1 ns simulations per replica. Detailed analysis reveals that explicit-solvent sampling provides increased accuracy relative to the previous GB-based method by preserving the native structure, providing a more realistic description of conformational flexibility of the hydrophobic cluster, and correctly modeling solvent mediated ion-pair interactions. Thus, we anticipate that the new technique will emerge as a practical tool to capture ionization equilibria while enabling an intimate view of ionization coupled conformational dynamics that is difficult to delineate with experimental techniques alone.

A Replica Exchange Molecular Dynamics Simulation of a Single Polyethylene Chain: Temperature Dependence of Structural Properties and Chain Conformational Study at the Equilibrium Melting Temperature

The conformational properties of a finite length polyethylene chain were explored over a wide range of temperatures using a replica exchange molecular dynamics simulation providing high quality simulation data representative for the equilibrium behavior of the chain molecule. The radial distribution function (RDF) and the structure factor S(q) of the chain as a function of temperature are analyzed in detail. The different characteristic peaks in the RDF and S(q) were assigned to specific distances in the chain and structural changes occurring with the temperature. In S(q), a peak characteristic for the order in the solid state was found and used to determine the equilibrium melting temperature. A detailed scaling analysis of the structure factor covering the full q range was performed according to the work of Hammouda. In the Θ region, a quantitative analysis of the full structure factor was done using the equivalent Kuhn chain, which enabled us to assign the Θ region of our chain and to demonstrate, in our particular case, the failure of the Gaussian chain approach. The chain conformational properties at the equilibrium melting temperature are discussed using conformational distribution functions, using the largest principal component of the radius of gyration and shape parameters as order parameters. We demonstrate that for the system studied here, the Landau free energy expression based on this conformational distribution information leads to erroneous conclusions concerning the thermodynamic transition behavior. Finally, we focus on the instantaneous conformational properties at the equilibrium melting temperature and give a detailed analysis of the conformational shapes using different shape parameters and a simulation snapshot. We show that the chain does not only take the lamellar rod-like and globular conformational shapes, typical of the solid and liquid states, but can also explore many other conformational states, including the toroidal conformational state. It is the first demonstration that a flexible molecule like PE can also take a toroidal conformational state, which is normally linked to stiffer chains.

Computing Relative Free Energies of Solvation using Single Reference Thermodynamic Integration Augmented with Hamiltonian Replica Exchange

This paper introduces an efficient single-topology variant of Thermodynamic Integration (TI) for computing relative transformation free energies in a series of molecules with respect to a single reference state. The presented TI variant that we refer to as Single-Reference TI (SR-TI) combines well-established molecular simulation methodologies into a practical computational tool. Augmented with Hamiltonian Replica Exchange (HREX), the SR-TI variant can deliver enhanced sampling in select degrees of freedom. The utility of the SR-TI variant is demonstrated in calculations of relative solvation free energies for a series of benzene derivatives with increasing complexity. Noteworthy, the SR-TI variant with the HREX option provides converged results in a challenging case of an amide molecule with a high (13-15 kcal/mol) barrier for internal cis/trans interconversion using simulation times of only 1 to 4 ns.

Highly polar environments catalyze the unfolding of PrP C helix 1

The first alpha-helix (H1) likely plays an important role in the conversion of the cellular prion protein (PrP(C)) into its pathogenic isoform (PrP(Sc)). In this conversion, H1 may either have to unfold or may represent a site of intermolecular contact. A recent molecular dynamics simulation suggested that H1 can unfold if it is detached from the protein core (Hirschberger et al. in Biophys J 90:3908, 2006). It has been hypothesized that the high dielectric constant epsilon (S) of the bulk water environment facilitates the unfolding of H1. To check this hypothesis, we performed a number of replica exchange molecular dynamics simulations of an H1 peptide in solvents of different epsilon (S). We found that the equilibrium helix fraction in water is less than 40%, in agreement with previous experimental findings, and that the helix unfolds much faster in water than in less polar solvents. The kinetically stabilizing effect of the organic solvents is largely unspecific and correlates well with their dielectric constant epsilon (S).

Error and efficiency of simulated tempering simulations

We derive simple analytical expressions for the error and computational efficiency of simulated tempering (ST) simulations. The theory applies to the important case of systems whose dynamics at long times is dominated by the slow interconversion between two metastable states. An extension to the multistate case is described. We show that the relative gain in efficiency of ST simulations over regular molecular dynamics (MD) or Monte Carlo (MC) simulations is given by the ratio of their reactive fluxes, i.e., the number of transitions between the two states summed over all ST temperatures divided by the number of transitions at the single temperature of the MD or MC simulation. This relation for the efficiency is derived for the limit in which changes in the ST temperature are fast compared to the two-state transitions. In this limit, ST is most efficient. Our expression for the maximum efficiency gain of ST simulations is essentially identical to the corresponding expression derived by us for replica exchange MD and MC simulations [E. Rosta and G. Hummer, J. Chem. Phys. 131, 165102 (2009)] on a different route. We find quantitative agreement between predicted and observed efficiency gains in a test against ST and replica exchange MC simulations of a two-dimensional Ising model. Based on the efficiency formula, we provide recommendations for the optimal choice of ST simulation parameters, in particular, the range and number of temperatures, and the frequency of attempted temperature changes.

Strengths and weaknesses of parallel tempering

Parallel tempering, also known as replica exchange Monte Carlo, is studied in the context of two simple free-energy landscapes. The first is a double-well potential defined by two macrostates separated by a barrier. The second is a "golf course" potential defined by microstates having two possible energies with exponentially more high-energy states than low-energy states. The equilibration time for replica exchange is analyzed for both systems. For the double-well system, parallel tempering with a number of replicas that scales as the square root of the barrier height yields exponential speedup of the equilibration time. On the other hand, replica exchange yields only marginal speedup for the golf course system. For the double-well system, the free-energy difference between the two wells has a large effect on the equilibration time. Nearly degenerate wells equilibrate much more slowly than strongly asymmetric wells. It is proposed that this difference in equilibration time may lead to a bias in measuring overlaps in spin glasses. These examples illustrate the strengths and weaknesses of replica exchange and may serve as a guide for understanding and improving the method in various applications.

Error and efficiency of replica exchange molecular dynamics simulations

We derive simple analytical expressions for the error and computational efficiency of replica exchange molecular dynamics (REMD) simulations (and by analogy replica exchange Monte Carlo simulations). The theory applies to the important case of systems whose dynamics at long times is dominated by the slow interconversion between two metastable states. As a specific example, we consider the folding and unfolding of a protein. The efficiency is defined as the rate with which the error in an estimated equilibrium property, as measured by the variance of the estimator over repeated simulations, decreases with simulation time. For two-state systems, this rate is in general independent of the particular property. Our main result is that, with comparable computational resources used, the relative efficiency of REMD and molecular dynamics (MD) simulations is given by the ratio of the number of transitions between the two states averaged over all replicas at the different temperatures, and the number of transitions at the single temperature of the MD run. This formula applies if replica exchange is frequent, as compared to the transition times. High efficiency of REMD is thus achieved by including replica temperatures in which the frequency of transitions is higher than that at the temperature of interest. In tests of the expressions for the error in the estimator, computational efficiency, and the rate of equilibration we find quantitative agreement with the results both from kinetic models of REMD and from actual all-atom simulations of the folding of a peptide in water.