Accelerating atomistic simulations of proteins using multiscale enhanced sampling with independent tempering

Introduction

Efficient sampling of relevant conformational space is crucial for reliable molecular simulations of protein structure, dynamics and interactions. Despite recent breakthroughs that greatly extend the reach of traditional physics-based atomistic simulations^1–4, enhanced sampling methods are often preferred, if not required, to bridge substantial gaps remained between simulation and biological timescales^5–9. This is particularly true for the studies of so-called intrinsically disordered proteins (IDPs)^10–12, where the molecular basis of function is largely determined by the disordered ensemble and how it may respond to various cellular signals such as cofactor binding and post-translational modifications^13–16. Reliable simulations of IDPs depend critically on the ability to sample the broad manifold of functionally relevant conformational space^17–21. Enhanced sampling aims to yield statistically meaningful structural ensembles with less computation, typically by accelerating the crossing of energy barriers. Temperature replica exchange (T-RE) in particular has emerged as a premier method for general enhanced sampling of biomolecular structures^22–24. It relies on (parallel) tempering and can accelerate sampling as long as the transition activation enthalpy is positive^25–31. However, it is also recognized that T-RE is severely limited by cooperative transitions such as protein folding^32–34. This limitation could not be completely overcome by T-RE variants that focus on improving replica exchange or diffusion in the temperature space³⁵. The fundamental reason is that the free energy barriers of cooperative conformational transitions are dominated by entropic components^{36, 37}, which renders tempering (i.e., raising the temperature) ineffective for driving faster transitions.

To overcome entropic sampling bottlenecks, a multiscale enhanced sampling (MSES) method has been developed, which takes advantage of the reduced conformational space as well as fast structural transition of coarse-grained (CG) models to accelerate the sampling of atomistic (AT) protein energy landscapes^38–40. Briefly, the simulation box contains both AT and CG representations of the system, which do not interact directly with each other. The AT and CG components of the hybrid system interact only through a MSES coupling potential (U_MSES) that is applied to selected “essential” degrees of freedom⁴⁰:

where U_hybrid is the potential energy of the hybrid system, and U_AT and U_CG are potentials for AT and CG models, respectively. The purpose of U_MSES is to couple the conformational fluctuations of AT and CG simulations, along the selected degrees of freedom. U_MSES usually takes the form of appropriate restrain potentials imposed on some measures of the structural difference between AT and CG representations. Such coupling scheme makes MSES more tolerant to inevitable artifacts of CG models and more scalable to large systems including explicit solvent all-atom simulations. The coupling parameter λ ranges from 0 to 1. When λ is large, the AT and CG models are strongly coupled, and faster structural transitions at the CG level are expected to accelerate AT conformational transitions. Hamiltonian replica exchange is utilized to remove the bias introduced by the coupling, such that proper thermodynamic ensembles could be recovered at both AT and CG levels at the limit condition of λ = 0. In addition, T-RE is incorporated into MSES to further enhance the sampling efficiency³⁸. The temperature/Hamiltonian replica exchange (H/T-RE) simulation protocol in MSES involves two parameters, temperature and λ, which are adjusted together when exchange between neighboring replicas occurs (i.e., 1D replica exchange).

Previous benchmark simulations using β-hairpins³⁸ or IDPs³⁹ have shown that MSES is able to significantly improve the efficiency of sampling atomistic conformations compared with T-RE. To further increase the efficiency of MSES, one may need to further accelerate large-scale conformational transitions at both AT and CG levels in strongly coupled conditions, and at the same time, allow such coupled, faster dynamics to be effectively communicated to the uncoupled limit. Here we describe a simple but effective strategy to achieve these objectives, by controlling the effective temperatures of AT and CG simulations independently. This method is hereafter referred to as MSES with independent tempering (MSES-IT). It employs the Hamiltonian scaling approach similar to the one used in replica exchange with solute tempering (REST)⁴¹ to set the effective temperature of CG model in all conditions near its melting temperature, where reversible structural transition rate at CG level is expected to be maximized. Applied to a nontrivial β-hairpin model peptide, MSES-IT proved to be more effective in driving and accelerating atomistic conformational rearrangements. Moreover, the rate of round trips between the low and high temperature extremes is greater in MSES-IT than in the original MSES, promoting more effective communication of accelerated dynamics between the coupled and uncoupled conditions. Together, these effects lead to faster convergence in the simulated structural ensembles by MSES-IT.

Method

In MSES-IT, the potential energy of the hybrid system in Eq. 1 is rewritten as,

where λ_AT and λ_CG are the energy scaling factors between 0 and 1 for the AT and CG simulations, such that their effective temperatures are increased by 1/λ_AT and 1/λ_CG, respectively. Both scaling factors can be independently adjusted in the H/T-RE protocol of MSES-IT to optimize the rate of conformational transitions and efficiency of dynamic coupling between AT and CG simulations. In this work, we test whether simply setting the CG simulation temperatures to near its melting temperature could provide improvement compared to the original MSES algorithm.

In the current MSES and MSES-IT simulations, the coupling potential was applied to Cα-Cα distances between nine residue pairs that form native contacts, which are believed to be the most essential degrees of freedom for GB1p reversible folding/unfolding transitions³⁸. To reduce energy penalty for large structural deviations between atomistic and CG models and ensure uniform exchange acceptance probability across all possible neighboring replicas, the coupling potential, U_MSES, was smoothly switched from a quadratic form for small structural deviations to the soft asymptote for large deviations³⁹:

where $d_i^{AT}$ and $d_i^{CG}$ are Cα-Cα distances for the i-th native contact in the atomistic and CG models, respectively, and Δd_i is the difference between them. At the distance threshold d_s, U_MSES begins to smoothly switch from a harmonic form to the soft asymptote (see Figure 2 of reference 39 for a plot of Eq. 3). The switching exponent s indicates how fast the limiting force, f_max, can be approached at large |Δd_i|. Parameters A and B are calculated by requiring both MSES coupling energy and force to be continuous when |Δd_i| is at the threshold distance d_s. Parameters used in the present work were k = 1.0 kcal/mol/Ų, s = 1, d_s = 2.0 Å, and f_max = 0.1 kcal/mol/Å. The coupling potential was further scaled by a coupling parameter λ (see Eq. 1) ranging from 0 to 1, and λ values of the eight replicas were 0, 0.10, 0.22, 0.35, 0.49, 0.64, 0.81 and 1.00, respectively.

Computational Details

A β-hairpin IDP GB1p (GEWTY DDATK TFTVT E, ~ 42% folded at 278 K⁴²), was used to evaluate the sampling efficiency of MSES-IT. Although this peptide is short, its folding time scale of 6–20 μs^{43, 44} is similar to that of many small proteins. Note that because of the highly flexible DATK loop in GB1p sequence, this peptide could sample a large conformational space, and its folding involves a large entropic cost (~ −40 cal/mol/K) as confirmed by previous experimental^44–46 and simulation⁴⁷ studies. Unlike helix/coil transitions, forming a beta-hairpin is highly cooperative and involves forming multiple long-range contacts. Therefore, it provides a very stringent system for benchmarking the efficiency of sampling algorithms, especially in overcoming the entropic barriers of protein folding. As a comparison, simulations using T-RE and the original MSES algorithm were also performed. All T-RE, MSES and MSES-IT simulations were carried out using CHARMM^{48, 49} together with a modified version of MMTSB⁵⁰. For each protocol, two independent simulations were performed, one starting from fully folded state (i.e., control run) and the other one starting from fully extended state (i.e., folding run). This allowed us to rigorously evaluate simulation convergence. Each simulation had eight replicas, with their temperatures spaced exponentially between 300 K and 450 K. λ_AT was set to 1.0 for all replicas in MSES-IT, while λ_CG was scaled such that the effective temperature of CG simulations in all replicas was 389 K, the melting temperature of the CG model used in this work. The atomistic model used in all simulations was implicit solvent GBSW force field⁵¹, while in MSES and MSES-IT simulations, a topology-based coarse-grained model, the Gō-like model calibrated in our previous work³⁸, was also used to drive the structural transitions of the atomistic model.

For all simulations, Langevin dynamics with a friction coefficient of 0.1 ps⁻¹ was performed. The equation of motion was integrated with a time step of 2 fs. SHAKE algorithm⁵² was used to constrain the length of all bonds involving hydrogen atoms. Exchange of replicas with their neighbors were attempted every 2 ps. Averaged exchange acceptance ratio was ~35% in all simulations. The production simulation time was 800 ns/replica in all case, making the current study one of the most extensive implicit solvent simulations of GB1p.

To quantify the number of folding/unfolding transitions of GB1p, we computed the number of native hydrogen bonds in the atomistic model during simulations. For this, we first computed the distances between seven pairs of backbone N and O atoms that are involved in hydrogen bonding in the fully folded state. In comparison to the fully folded state, distance violations were first calculated for the seven pairs of N-O atoms in reference to the distances in the native folded structure, which then can be converted to number of native hydrogen bonds using Eq. 4:

where Δd_i is the distance violation for the i-th pair of N-O atoms, and N is the number of native hydrogen bonds. The summand smoothly switches from 1 to 0 as the distance violation Δd_i increases, allowing the total number of hydrogen bond formed be continuously ranging from 0 to 7. This helps avoid observing spurious transitions between folded and unfolded states in cases where a single distance cutoff is directly used to identify contact formation. The atomistic structural transition rate was calculated as the total number of reversible folding-unfolding transitions divided by the total simulation time, where folded and unfolded states were defined as N ≥ 6.0 and N ≤ 1.0, respectively (see Figure 1). Results reported in this study are averaged values over all replicas.

To quantify convergence rate of each sampling algorithm, we computed simulation error as a function of trajectory length. Each trajectory was first divided into multiple segments of length t. The probability distribution of native hydrogen bond number P(N) was calculated for each segment. The associated mean absolute error was computed with respect to the reference P(N), which was calculated using the whole trajectory and averaged over all six (T-RE, MSES, MSES-IT control and folding) simulations. Average and standard deviation of mean absolute error were reported here by considering all trajectory segments of length t.

Results

We first validate that both MSES and MSES-IT are able to remove the bias of MSES coupling between AT and CG simulations and generate correct statistical ensembles at the limit condition of λ = 0 via Hamiltonian replica exchange. In Figure 1, we compare probability distributions of the number of native hydrogen bonds formed in the atomistic model obtained from T-RE, MSES and MSES-IT control and folding simulations. Because of the very long simulation time of 800 ns, T-RE is able to achieve apparent convergence in the simulated ensembles derived from control and folding simulations (Figure1, blue traces), thus providing a good reference for assessing the correctness of atomistic ensembles derived from MSES protocols ((Figure 1, green and red traces). The root mean squared differences (RMSDs) among all distributions is no greater than 0.02. That is, there is no systematic deviation in results derived from either MSES or MSES-IT simulations, indicating that no thermodynamic bias was introduced at the atomistic level due to the coupling to CG simulations in MSES or MSES-IT protocols.

Next, we examine the efficacy of MSES-IT in accelerating large-scale conformational transitions at the atomistic level. Figure 2 summarizes the reversible folding-unfolding transition rates of GB1p in the GBSW implicit solvent. Clearly, coupling the atomistic model to even the topology-based CG model in either MSES or MSES-IT simulations significantly accelerated structural transitions at the atomistic level, generating over 20-fold more reversible conformational transitions compared to T-RE (Figure 2A). Importantly, the averaged folding-unfolding transition rates over all replicas are the highest in MSES-IT simulations. Closer inspection shows that individual replicas consistently underwent faster reversible folding-unfolding transitions in MSES-IT (Figure 2B). Therefore, setting the effective temperature of CG model to near its melting temperature in MSES-IT did further accelerate atomistic structural transitions compared to the original MSES, where the temperatures of AT and CG simulations are always the same. Accelerated large scale transitions should allow faster convergence of the resulting atomistic conformational ensembles. This will be further evaluated later.

Another important indicator of replica exchange sampling efficiency is the replica mixing rate, which can be measured as the rate of replica round trips between the lowest and highest conditions. There have been many efforts towards maximizing this quantity to increase the sampling efficiency of T-RE, such as adaptively adjusting the distribution of simulation temperatures^{53, 54} or using biasing potentials^{55, 56}. With independent tempering, MSES-IT achieves the fastest replica mixing rate among three simulation protocols (see Figure 3). In MSES, structural divergence between atomistic and CG models may lead to replicas being trapped at low-temperature (and weakly coupled) conditions, thus lowering replica mixing rate³⁹. As illustrated in Figure 4, faster folding/unfolding transitions at both CG and AT levels in MSES-IT reduce the structural divergence between the two models, particularly under the conditions of low temperature/weak MSES coupling (small λ_MSES). Therefore, it is less likely for replicas being trapping in low-temperature and weakly coupled regime, thus improving replica mixing efficiency. Rapid replica mixing, in turn, allows the fast dynamics driven by CG models under the coupled conditions to be more efficiently exchanged to the uncoupled limit, further increasing sampling efficiency.

Faster replica mixing together with larger numbers of atomistic conformational transitions, driven by coupling to CG simulations near the melting temperature, further improve the ability of MSES-IT to generate converged structural ensembles. As shown in Figure 5, the mean absolute error in the calculated probability distribution of native hydrogen bond number of GB1p decays to 0.02 within 400 ns MSES-IT simulations, which is faster than both T-RE and MSES. Note that the current error estimation is noisy and associated with large uncertainties, especially when trajectory length is long and there are limited number of available data points. Nonetheless, fitting the curve using a single exponential function f(x) = Ae^−x/τ also confirms that the error in MSES-IT decays fastest with smallest half-time τ. This demonstrates the superior ability of MSES-IT in generating well converged protein conformational ensembles.

Conclusions

MSES provides an effective strategy to exploit faster dynamics in CG simulations that accelerate the sampling of atomistic conformational space and overcome entropic barriers of cooperative conformational transitions. Maximizing the MSES sampling efficiency requires the CG model to drive atomistic structural transitions as fast as possible. Furthermore, faster dynamics in the coupled conditions needs to be efficiently exchanged to the uncoupled/unbiased condition. To better achieve these goals, we have developed a simple but effective strategy by allowing independent tempering of the AT and CG simulations in MSES. Briefly, the effective temperatures of AT and CG simulations in all conditions are independently controlled by scaling the AT and CG potentials accordingly. In particular, we showed that setting the CG simulations to the melting temperature would maximize the structural transition rate at the CG level, providing additional acceleration of reversible folding and unfolding transitions of a non-trivial β-hairpin peptide GB1p. The benchmark showed that the replica mixing rate was also increased, allowing more efficient communication of dynamic transitions between the coupled and uncoupled conditions. These improvements did translate into faster convergence in the atomistic structural ensembles, supporting that MSES-IT provides superior sampling efficiency over the original MSES scheme.

Although structural transitions usually are the fastest near the melting temperature, it is possible that melting temperature may not be the optimal effective temperature for CG model to maximally accelerate atomistic structural transitions in the coupled conditions. The ability to enhance atomistic structural transition rate in MSES-IT also depends on the CG model used. Further improvement of CG modeling quality should be able to derive even faster conformational rearrangement⁵⁷. There is room for improving the MSES-IT sampling efficiency by further optimization of the H/T-RE protocol in conjunction with better CG models. In principal we could set the effective temperatures of CG simulations to different values under different conditions. Nevertheless, current work demonstrates that independent tempering of the AT and CG models is an effective strategy to more effectively drive atomistic structural transitions and improve MSES sampling efficiency.

Data Availability Statements

The data that support the findings of this study are available from the corresponding author upon reasonable request.