Generalized ensemble and tempering simulations: a unified view

A density of states-based approach to determine temperature-dependent aggregation rates

Here, we establish an approach to determine temperature-dependent aggregation rates in terms of thermostatistical quantities, which can be obtained directly from flat-histogram and statistical temperature algorithms considering the density of states of the system. Our approach is validated through simulations of an Ising-like model with anisotropically interacting particles at temperatures close to its first-order phase transition. Quantitative comparisons between the numerically obtained forward and reverse rates to approximate analytical expressions corroborate its use as a model-independent approach.

Exploring Protocols to Build Reservoirs to Accelerate Temperature Replica Exchange MD Simulations

Temperature replica exchange molecular dynamics (REMD) is a widely used enhanced sampling method for accelerating biomolecular simulations. During the past 2 decades, several variants of REMD have been developed to further improve the rate of conformational sampling of REMD. One such variant, reservoir REMD (RREMD), was shown to improve the rate of conformational sampling by around 5-20×. Despite the significant increase in the sampling speed, RREMD methods have not been widely used because of the difficulties in building the reservoir and also because of the code not being available on the graphics processing units (GPUs). In this work, we ported the Amber RREMD code onto GPUs making it 20× faster than the central processing unit code. Then, we explored protocols for building Boltzmann-weighted reservoirs as well as non-Boltzmann reservoirs and tested how each choice affects the accuracy of the resulting RREMD simulations. We show that, using the recommended protocols outlined here, RREMD simulations can accurately reproduce Boltzmann-weighted ensembles obtained by much more expensive conventional temperature-based REMD simulations, with at least 15× faster convergence rates even for larger proteins (>50 amino acids) compared to conventional REMD.

Flat-histogram method comparison on the two-dimensional Ising model

We compare the convergence of several flat-histogram methods applied to the two-dimensional Ising model, including the recently introduced stochastic approximation with a dynamic update factor (SAD) method. We compare this method to the Wang-Landau (WL) method, the 1/t variant of the WL method, and standard stochastic approximation Monte Carlo (SAMC). In addition, we consider a procedure WL followed by a "production run" with fixed weights that refines the estimation of the entropy. We find that WL followed by a production run does converge to the true density of states, in contrast to pure WL. Three of the methods converge robustly: SAD, 1/t-WL, and WL followed by a production run. Of these, SAD does not require a priori knowledge of the energy range. This work also shows that WL followed by a production run performs superior to other forms of WL while ensuring both ergodicity and detailed balance.

Role of substrate-product frustration on enzyme functional dynamics

Natural enzymes often have enormous catalytic power developed by evolution. Revealing the underlying physical strategy used by enzymes to achieve high catalysis efficiency is one of the central focuses in the field of biological physics. Our recent work demonstrated that multisubstrate enzymes can utilize steric frustration encountered in the substrate-product cobound complex to overcome the bottleneck of the enzymatic cycle [W. Li et al., Phys. Rev. Lett. 122, 238102 (2019)10.1103/PhysRevLett.122.238102]. However, the key atomic-level interactions by which the steric frustration contributes to the enzymatic cycle remain elusive. In this work we study the microscopic mechanism for the role of the substrate-product frustration on the key physical steps in the enzymatic cycle of adenylate kinase (AdK), a multisubstrate enzyme catalyzing the reversible phosphoryl transfer reaction ATP+AMP⇋ADP+ADP. By using atomistic molecular dynamics simulations with enhanced sampling, we showed that the competitive interactions from the phosphate groups of the substrate ATP and product ADP in the ATP-ADP cobound complex of the AdK lead to local frustration in the binding pockets. Such local frustration disrupts the hydrogen bond network around the binding pockets, which causes lowered barrier height for the opening of the enzyme conformations and expedited release of the bottleneck product ADP. Our results directly demonstrated from the atomistic level that the local frustration in the active sites of the enzyme can be utilized to facilitate the key physical steps of the enzymatic cycle, providing numerical evidence to the predictions of the previous theoretical work.

Replica sub-permutation method for molecular dynamics and monte carlo simulations

We propose an improvement of the replica-exchange and replica-permutation methods, which we call the replica sub-permutation method (RSPM). Instead of considering all permutations, this method uses a new algorithm referred to as sub-permutation to perform parameter transition. The RSPM succeeds in reducing the number of combinations between replicas and parameters without the loss of sampling efficiency. For comparison, we applied the replica sub-permutation, replica-permutation, and replica-exchange methods to a β-hairpin mini protein, chignolin, in explicit water. We calculated the transition ratio and number of tunneling events in the parameter space, the number of folding-unfolding events, the autocorrelation function, and the autocorrelation time as measures of sampling efficiency. The results indicate that among the three methods, the proposed RSPM is the most efficient in both parameter and conformational spaces. © 2019 Wiley Periodicals, Inc.

Riemann metric approach to optimal sampling of multidimensional free-energy landscapes

Exploring the free-energy landscape along reaction coordinates or system parameters λ is central to many studies of high-dimensional model systems in physics, e.g., large molecules or spin glasses. In simulations this usually requires sampling conformational transitions or phase transitions, but efficient sampling is often difficult to attain due to the roughness of the energy landscape. For Boltzmann distributions, crossing rates decrease exponentially with free-energy barrier heights. Thus, exponential acceleration can be achieved in simulations by applying an artificial bias along λ tuned such that a flat target distribution is obtained. A flat distribution is, however, an ambiguous concept unless a proper metric is used and is generally suboptimal. Here we propose a multidimensional Riemann metric, which takes the local diffusion into account, and redefine uniform sampling such that it is invariant under nonlinear coordinate transformations. We use the metric in combination with the accelerated weight histogram method, a free-energy calculation and sampling method, to adaptively optimize sampling toward the target distribution prescribed by the metric. We demonstrate that for complex problems, such as molecular dynamics simulations of DNA base-pair opening, sampling uniformly according to the metric, which can be calculated without significant computational overhead, improves sampling efficiency by 50%-70%.

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.

Multifunnel Landscape of the Fold-Switching Protein RfaH-CTD

Proteins such as the transcription factor RfaH can change biological function by switching between distinct three-dimensional folds. RfaH regulates transcription if the C-terminal domain folds into a double helix bundle and promotes translation when this domain assumes a β-barrel form. This fold-switch has been also observed for the isolated C-terminal domain, dubbed by us as RfaH-C-terminal domain (RfaH-CTD), and is studied here with a variant of the replica-exchange-with-tunneling approach recently introduced by us. We use the enhanced sampling properties of this technique to map the free-energy landscape of RfaH-CTD and to propose a mechanism for the conversion process.

Efficient Round-Trip Time Optimization for Replica-Exchange Enveloping Distribution Sampling (RE-EDS)

Replica-exchange enveloping distribution sampling (RE-EDS) allows the efficient estimation of free-energy differences between multiple end-states from a single molecular dynamics (MD) simulation. In EDS, a reference state is sampled, which can be tuned by two types of parameters, i.e., smoothness parameters(s) and energy offsets, such that all end-states are sufficiently sampled. However, the choice of these parameters is not trivial. Replica exchange (RE) or parallel tempering is a widely applied technique to enhance sampling. By combining EDS with the RE technique, the parameter choice problem could be simplified and the challenge shifted toward an optimal distribution of the replicas in the smoothness-parameter space. The choice of a certain replica distribution can alter the sampling efficiency significantly. In this work, global round-trip time optimization (GRTO) algorithms are tested for the use in RE-EDS simulations. In addition, a local round-trip time optimization (LRTO) algorithm is proposed for systems with slowly adapting environments, where a reliable estimate for the round-trip time is challenging to obtain. The optimization algorithms were applied to RE-EDS simulations of a system of nine small-molecule inhibitors of phenylethanolamine N-methyltransferase (PNMT). The energy offsets were determined using our recently proposed parallel energy-offset (PEOE) estimation scheme. While the multistate GRTO algorithm yielded the best replica distribution for the ligands in water, the multistate LRTO algorithm was found to be the method of choice for the ligands in complex with PNMT. With this, the 36 alchemical free-energy differences between the nine ligands were calculated successfully from a single RE-EDS simulation 10 ns in length. Thus, RE-EDS presents an efficient method for the estimation of relative binding free energies.

Simulating Protein Fold Switching by Replica Exchange with Tunneling

Recent experiments suggest that an amino acid sequence encodes not only the native fold of a protein but also other forms that are essential for its function or are important during folding or association. These various forms populate a multifunnel folding and association landscape where mutations, changes in environment, or interaction with other molecules switch between the encoded folds. We introduce replica exchange with tunneling as a way to efficiently simulate switching between distinct folds of proteins and protein aggregates. The correctness and efficiency of our approach are demonstrated in a series of simulations covering a wide range of proteins, from a small 11-residue large designed peptide to two 56-residue large mutants of the A and B domains of protein G.

Replica-exchange-with-tunneling for fast exploration of protein landscapes

While the use of replica-exchange molecular dynamics in protein simulations has become ubiquitous, its utility is limited in many practical applications. We propose to overcome some shortcomings that hold back its use in settings such as multi-scale or explicit solvent simulations by integrating ideas of hybrid MC/MD into the replica-exchange protocol. This Replica-Exchange-with-Tunneling method is tested by simulating the Trp-cage protein, a system often used in molecular biophysics for testing sampling techniques.

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.

Communication: Multiple atomistic force fields in a single enhanced sampling simulation

The main concerns of biomolecular dynamics simulations are the convergence of the conformational sampling and the dependence of the results on the force fields. While the first issue can be addressed by employing enhanced sampling techniques such as simulated tempering or replica exchange molecular dynamics, repeating these simulations with different force fields is very time consuming. Here, we propose an automatic method that includes different force fields into a single advanced sampling simulation. Conformational sampling using three all-atom force fields is enhanced by simulated tempering and by formulating the weight parameters of the simulated tempering method in terms of the energy fluctuations, the system is able to perform random walk in both temperature and force field spaces. The method is first demonstrated on a 1D system and then validated by the folding of the 10-residue chignolin peptide in explicit water.

Simulating Replica Exchange: Markov State Models, Proposal Schemes, and the Infinite Swapping Limit

Replica exchange molecular dynamics is a multicanonical simulation technique commonly used to enhance the sampling of solvated biomolecules on rugged free energy landscapes. While replica exchange is relatively easy to implement, there are many unanswered questions about how to use this technique most efficiently, especially because it is frequently the case in practice that replica exchange simulations are not fully converged. A replica exchange cycle consists of a series of molecular dynamics steps of a set of replicas moving under different Hamiltonians or at different thermodynamic states followed by one or more replica exchange attempts to swap replicas among the different states. How the replica exchange cycle is constructed affects how rapidly the system equilibrates. We have constructed a Markov state model of replica exchange (MSMRE) using long molecular dynamics simulations of a host-guest binding system as an example, in order to study how different implementations of the replica exchange cycle can affect the sampling efficiency. We analyze how the number of replica exchange attempts per cycle, the number of MD steps per cycle, and the interaction between the two parameters affects the largest implied time scale of the MSMRE simulation. The infinite swapping limit is an important concept in replica exchange. We show how to estimate the infinite swapping limit from the diagonal elements of the exchange transition matrix constructed from MSMRE "simulations of simulations" as well as from relatively short runs of the actual replica exchange simulations.

Enhancing pairwise state-transition weights: A new weighting scheme in simulated tempering that can minimize transition time between a pair of conformational states

Simulated tempering (ST) is a widely used enhancing sampling method for Molecular Dynamics simulations. As one expanded ensemble method, ST is a combination of canonical ensembles at different temperatures and the acceptance probability of cross-temperature transitions is determined by both the temperature difference and the weights of each temperature. One popular way to obtain the weights is to adopt the free energy of each canonical ensemble, which achieves uniform sampling among temperature space. However, this uniform distribution in temperature space may not be optimal since high temperatures do not always speed up the conformational transitions of interest, as anti-Arrhenius kinetics are prevalent in protein and RNA folding. Here, we propose a new method: Enhancing Pairwise State-transition Weights (EPSW), to obtain the optimal weights by minimizing the round-trip time for transitions among different metastable states at the temperature of interest in ST. The novelty of the EPSW algorithm lies in explicitly considering the kinetics of conformation transitions when optimizing the weights of different temperatures. We further demonstrate the power of EPSW in three different systems: a simple two-temperature model, a two-dimensional model for protein folding with anti-Arrhenius kinetics, and the alanine dipeptide. The results from these three systems showed that the new algorithm can substantially accelerate the transitions between conformational states of interest in the ST expanded ensemble and further facilitate the convergence of thermodynamics compared to the widely used free energy weights. We anticipate that this algorithm is particularly useful for studying functional conformational changes of biological systems where the initial and final states are often known from structural biology experiments.

Coupled Nosé-Hoover equations of motion to implement a fluctuating heat-bath temperature

The Nosé-Hoover (NH) equation provides a universal and powerful computer simulation protocol to realize an equilibrium canonical temperature for a target physical system. Here we demonstrate a general formalism to couple such NH equations. We provide a coupled NH equation that is constructed by coupling the NH equation of a target physical system and the NH equation of a temperature system. Thus, in contrast to the conventional single NH equation, the heat-bath temperature is a dynamical variable. The temperature fluctuations are not ad hoc, but instead are generated by the newly defined temperature system, and the statistical distribution of the temperature is completely described with an arbitrarily given probability function. The current equations of motion thus describe the physical system that develops with a predistributed fluctuating temperature, which allows enhanced sampling of the physical system. Since the total system is governed by a prescribed distribution, the equilibrium of the physical system is also reconstructed by reweighting. We have formulated a scheme for specifically setting the distribution of the dynamical inverse temperature and demonstrate the statistical relationship between the dynamical and physical temperatures. The statistical features, dynamical properties, and sampling abilities of the current method are demonstrated via the distributions, trajectories, dynamical correlations, and free energy landscapes for both a model system and a biomolecular system. These results indicated that the current coupled NH scheme works well.

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.

Towards an Optimal Flow: Density-of-States-Informed Replica-Exchange Simulations

Replica exchange (RE) is one of the most popular enhanced-sampling simulations technique in use today. Despite widespread successes, RE simulations can sometimes fail to converge in practical amounts of time, e.g., when sampling around phase transitions, or when a few hard-to-find configurations dominate the statistical averages. We introduce a generalized RE scheme, density-of-states-informed RE, that addresses some of these challenges. The key feature of our approach is to inform the simulation with readily available, but commonly unused, information on the density of states of the system as the RE simulation proceeds. This enables two improvements, namely, the introduction of resampling moves that actively move the system towards equilibrium and the continual adaptation of the optimal temperature set. As a consequence of these two innovations, we show that the configuration flow in temperature space is optimized and that the overall convergence of RE simulations can be dramatically accelerated.

Simulations in generalized ensembles through noninstantaneous switches

Generalized-ensemble simulations, such as replica exchange and serial generalized-ensemble methods, are powerful simulation tools to enhance sampling of free energy landscapes in systems with high energy barriers. In these methods, sampling is enhanced through instantaneous transitions of replicas, i.e., copies of the system, between different ensembles characterized by some control parameter associated with thermodynamical variables (e.g., temperature or pressure) or collective mechanical variables (e.g., interatomic distances or torsional angles). An interesting evolution of these methodologies has been proposed by replacing the conventional instantaneous (trial) switches of replicas with noninstantaneous switches, realized by varying the control parameter in a finite time and accepting the final replica configuration with a Metropolis-like criterion based on the Crooks nonequilibrium work (CNW) theorem. Here we revise these techniques focusing on their correlation with the CNW theorem in the framework of Markovian processes. An outcome of this report is the derivation of the acceptance probability for noninstantaneous switches in serial generalized-ensemble simulations, where we show that explicit knowledge of the time dependence of the weight factors entering such simulations is not necessary. A generalized relationship of the CNW theorem is also provided in terms of the underlying equilibrium probability distribution at a fixed control parameter. Illustrative calculations on a toy model are performed with serial generalized-ensemble simulations, especially focusing on the different behavior of instantaneous and noninstantaneous replica transition schemes.

Computational Calorimetry: High-Precision Calculation of Host-Guest Binding Thermodynamics

We present a strategy for carrying out high-precision calculations of binding free energy and binding enthalpy values from molecular dynamics simulations with explicit solvent. The approach is used to calculate the thermodynamic profiles for binding of nine small molecule guests to either the cucurbit[7]uril (CB7) or β-cyclodextrin (βCD) host. For these systems, calculations using commodity hardware can yield binding free energy and binding enthalpy values with a precision of ∼0.5 kcal/mol (95% CI) in a matter of days. Crucially, the self-consistency of the approach is established by calculating the binding enthalpy directly, via end point potential energy calculations, and indirectly, via the temperature dependence of the binding free energy, i.e., by the van't Hoff equation. Excellent agreement between the direct and van't Hoff methods is demonstrated for both host-guest systems and an ion-pair model system for which particularly well-converged results are attainable. Additionally, we find that hydrogen mass repartitioning allows marked acceleration of the calculations with no discernible cost in precision or accuracy. Finally, we provide guidance for accurately assessing numerical uncertainty of the results in settings where complex correlations in the time series can pose challenges to statistical analysis. The routine nature and high precision of these binding calculations opens the possibility of including measured binding thermodynamics as target data in force field optimization so that simulations may be used to reliably interpret experimental data and guide molecular design.

Toward Improved Force-Field Accuracy through Sensitivity Analysis of Host-Guest Binding Thermodynamics

Improving the capability of atomistic computer models to predict the thermodynamics of noncovalent binding is critical for successful structure-based drug design, and the accuracy of such calculations remains limited by nonoptimal force field parameters. Ideally, one would incorporate protein-ligand affinity data into force field parametrization, but this would be inefficient and costly. We now demonstrate that sensitivity analysis can be used to efficiently tune Lennard-Jones parameters of aqueous host-guest systems for increasingly accurate calculations of binding enthalpy. These results highlight the promise of a comprehensive use of calorimetric host-guest binding data, along with existing validation data sets, to improve force field parameters for the simulation of noncovalent binding, with the ultimate goal of making protein-ligand modeling more accurate and hence speeding drug discovery.

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.

Theory of Adaptive Optimization for Umbrella Sampling

We present a theory of adaptive optimization for umbrella sampling. With the analytical bias force constant obtained from the constrained thermodynamic length along the reaction coordinate, the windows are distributed to optimize the overlap between neighbors. Combining with the replica exchange method, we propose a method of adaptive window exchange umbrella sampling. The efficiency gain in sampling by the present method originates from the optimal window distribution, in which windows are concentrated to the region where the free energy is steep, as well as consequently improved random walk.

Replica exchange with guided annealing for accelerated sampling of disordered protein conformations

We critically examine a recently proposed convective replica exchange (cRE) method for enhanced sampling of protein conformation based on theoretical and numerical analysis. The results demonstrate that cRE and related replica exchange with guided annealing (RE-GA) schemes lead to unbalanced exchange attempt probabilities and break detailed balance whenever the system undergoes slow conformational transitions (relative to the temperature diffusion timescale). Nonetheless, numerical simulations suggest that approximate canonical ensembles can be generated for systems with small conformational transition barriers. This suggests that RE-GA maybe suitable for simulating intrinsically disordered proteins, an important class of newly recognized functional proteins. The efficacy of RE-GA is demonstrated by calculating the conformational ensembles of intrinsically disordered kinase inducible domain protein. The results show that RE-GA helps the protein to escape nonspecific compact states more efficiently and provides several fold speedups in generating converged and largely correct ensembles compared to the standard temperature RE.

A Hamiltonian replica exchange molecular dynamics (MD) method for the study of folding, based on the analysis of the stabilization determinants of proteins

Herein, we present a novel Hamiltonian replica exchange protocol for classical molecular dynamics simulations of protein folding/unfolding. The scheme starts from the analysis of the energy-networks responsible for the stabilization of the folded conformation, by means of the energy-decomposition approach. In this framework, the compact energetic map of the native state is generated by a preliminary short molecular dynamics (MD) simulation of the protein in explicit solvent. This map is simplified by means of an eigenvalue decomposition. The highest components of the eigenvector associated with the lowest eigenvalue indicate which sites, named "hot spots", are likely to be responsible for the stability and correct folding of the protein. In the Hamiltonian replica exchange protocol, we use modified force-field parameters to treat the interparticle non-bonded potentials of the hot spots within the protein and between protein and solvent atoms, leaving unperturbed those relative to all other residues, as well as solvent-solvent interactions. We show that it is possible to reversibly simulate the folding/unfolding behavior of two test proteins, namely Villin HeadPiece HP35 (35 residues) and Protein A (62 residues), using a limited number of replicas. We next discuss possible implications for the study of folding mechanisms via all atom simulations.

Charge-leveling and proper treatment of long-range electrostatics in all-atom molecular dynamics at constant pH

Recent development of constant pH molecular dynamics (CpHMD) methods has offered promise for adding pH-stat in molecular dynamics simulations. However, until now the working pH molecular dynamics (pHMD) implementations are dependent in part or whole on implicit-solvent models. Here we show that proper treatment of long-range electrostatics and maintaining charge neutrality of the system are critical for extending the continuous pHMD framework to the all-atom representation. The former is achieved here by adding forces to titration coordinates due to long-range electrostatics based on the generalized reaction field method, while the latter is made possible by a charge-leveling technique that couples proton titration with simultaneous ionization or neutralization of a co-ion in solution. We test the new method using the pH-replica-exchange CpHMD simulations of a series of aliphatic dicarboxylic acids with varying carbon chain length. The average absolute deviation from the experimental pK(a) values is merely 0.18 units. The results show that accounting for the forces due to extended electrostatics removes the large random noise in propagating titration coordinates, while maintaining charge neutrality of the system improves the accuracy in the calculated electrostatic interaction between ionizable sites. Thus, we believe that the way is paved for realizing pH-controlled all-atom molecular dynamics in the near future.

Two Dimensional Window Exchange Umbrella Sampling for Transmembrane Helix Assembly

The method of window exchange umbrella sampling molecular dynamics (WEUSMD) with a pre-optimized parameter set was recently used to obtain the most probable conformations and the energetics of transmembrane (TM) helix assembly of a generic TM sequence. When applied to glycophorin A TM domain (GpA-TM) using the restraint potentials along the helix-helix distance, however, tight interfacial packing of GpA-TM resulted in insufficient conformational sampling at short helix-helix separation. This sampling issue is addressed by extending the WEUSMD into two dimensions with the restraint potentials along the helix-helix distance and crossing angle. The two-dimensional WEUSMD results demonstrate that the incomplete sampling in the one-dimensional WEUSMD arises from high barriers along the crossing angle between the GpA-TM helices. Together with the faster convergence in both the assembled conformations and the potential of mean force, the 2D-WEUSMD can be a general and efficient approach in computational studies of TM helix assembly.

A convective replica-exchange method for sampling new energy basins

Replica-exchange is a powerful simulation method for sampling the basins of a rugged energy landscape. The replica-exchange method's sampling is efficient because it allows replicas to perform round trips in temperature space, thereby visiting both low and high temperatures in the same simulation. However, replicas have a diffusive walk in temperature space, and the round trip rate decreases significantly with the system size. These drawbacks make convergence of the simulation even more difficult than it already is when bigger systems are tackled. Here, we present a simple modification of the exchange method. In this method, one of the replicas steadily raises or lowers its temperature. We tested the convective replica-exchange method on three systems of varying complexity: the alanine dipeptide in implicit solvent, the GB1 β-hairpin in explicit solvent and the Aβ(25-35) homotrimer in a coarse grained representation. For the highly frustrated Aβ(25-35) homotrimer, the proposed "convective" replica-exchange method is twice as fast as the standard method. It discovered 24 out of 27 free-energy basins in less than 500 ns. It also prevented the formation of groups of replicas that usually form on either side of an exchange bottleneck, leading to a more efficient sampling of new energy basins than in the standard method.

Modeling Structural Flexibility of Proteins with Go-Models

Structure-based models are an efficient tool for folding studies of proteins since by construction their energy landscape is only minimal frustrated. However, their intrinsic drawback is a lack of structural flexibility as usually only one target structure is employed to construct the potentials. Hence, a Go-model may not capture differences in mutation-induced protein dynamics, if - as in the case of the disease-related A629P mutant of the Menkes protein ATP7A - the structural differences between mutant and wild type are small. In this work, we introduced three implementations of Go-models that take into account the flexibility of proteins in the NMR ensemble. Comparing the wild type and the mutant A629P of the 75-residue large 6th domain Menkes protein, we find that these new Go-potentials lead to broader distributions than Go-models relying on a single member of the NMR-ensemble. This allows us to detect the transient unfolding of a loosely formed β1β4-sheet in the mutant protein. Our results are consistent with previous simulations using physical force field and an explicit solvent, and suggests a mechanism by which these mutations cause Menkes disease. In addition, the improved Go-models suggest differences in the folding pathway between wild type and mutant, an observation that was not accessible to simulations of this 75-residue protein with a physical all-atom force field and explicit solvent.

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.

Folding simulations of the A and B domains of protein G

We study wild type and mutants of the A and B domain of protein G using all-atom Go-models. Our data substantiate the usefulness of such simulation for probing the folding mechanism of proteins and demonstrate that multifunnel versions of such models also allow probing of more complicated funnel landscapes. In our case, such models reproduce the experimentally observed distributions of the GA98 and GB98 mutants which differ only by one residue but fold into different structures. They also reveal details on the folding mechanism in these two proteins.

Velocity-scaling optimized replica exchange molecular dynamics of proteins in a hybrid explicit/implicit solvent

We propose a scheme for replica exchange molecular dynamics of proteins in explicit solvent that minimizes the number of required replicas using velocity rescaling. Our approach relies on a hybrid method where the protein evolves at each temperature in an explicit solvent, but replica exchange moves utilize an implicit solvent term. The two terms are coupled through the velocity rescaling. We test the efficiency of this approach for a common test case, the trp-cage protein.

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.

Serial tempering without exchange

Serial tempering is a computational method that turns the temperature T (or more generally any independent λ parameter) into a dynamical variable. It is shown that, under conditions for which this variable is fast, serial tempering is equivalent to the umbrella sampling method with a single effective potential. This equivalence is demonstrated using both a small one-dimensional system and a small solvated peptide. The suggestion is then made to replace the serial tempering protocol with the equivalent umbrella sampling calculation. This approach, serial tempering without exchange (STeWiE), has the same performance as serial tempering in the limit that exchanges are frequent, is simpler to implement, and has fewer adjustable parameters than conventional serial tempering. The equivalence of serial tempering and STeWiE also provides a convenient route for estimating and optimizing the performance of serial tempering simulations and other generalized-ensemble methods.

Balancing simulation accuracy and efficiency with the Amber united atom force field

We have analyzed the quality of a recently proposed Amber united-atom model and its overall efficiency in ab initio folding and thermodynamic sampling of two stable beta-hairpins. It is found that the mean backbone structures are quite consistent between the simulations in the united-atom and its corresponding all-atom models in Amber. More importantly, the simulated beta turns are also consistent between the two models. Finally, the chemical shifts on H alpha are highly consistent between simulations in the two models, although the simulated chemical shifts are lower than experiment, indicating less structured peptides, probably due to the omission of the hydrophobic term in the simulations. More interestingly, the stabilities of both beta-hairpins at room temperature are similar to those derived from the NMR measurement, whether the united-atom or the all-atom model is used. Detailed analysis shows high percentages of backbone torsion angles within the beta region and high percentages of native contacts. Given the reasonable quality of the united-atom model with respect to experimental data, we have further studied the simulation efficiency of the united-atom model over the all-atom model. Our data shows that the united-atom model is a factor of 6-8 faster than the all-atom model as measured with the ab initio first pass folding time for the two tested beta-hairpins. Detailed structural analysis shows that all ab initio folded trajectories enter the native basin, whether the united-atom model or the all-atom model is used. Finally, we have also studied the simulation efficiency of the united-atom model as measured in terms of how fast thermodynamic convergence can be achieved. It is apparent that the united-atom simulations reach convergence faster than the all-atom simulations with respect to both mean potential energies and mean native contacts. These findings show that the efficiency of the united-atom model is clearly beyond the per-step dynamics simulation of about 2 over the all-atom model. Thus, reasonable reduction of a protein model can be achieved with improved sampling efficiency while still preserving a high level of accuracy for applications in both ab initio folding and thermodynamic sampling. This study motivates us to develop more simplified protein models with sufficient consistency with the all-atom models for enhanced conformational sampling.

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.

Microcanonical replica exchange molecular dynamics simulation of proteins

We present microcanonical replica exchange molecular dynamics simulations as an alternative to canonical ones. Its advantage is the easily tunable high acceptance rate for replica exchange. We present the theory, comment on its actual implementation, and demonstrate its application for a common test case, the trp-cage protein.

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.