Computational Estimation of Microsecond to Second Atomistic Folding Times

Abstract

Despite the development of massively parallel computing hardware including inexpensive graphics processing units (GPUs), it has remained infeasible to simulate the folding of atomistic proteins at room temperature using conventional molecular dynamics (MD) beyond the microsecond scale. Here, we report the folding of atomistic, implicitly solvated protein systems with folding times τ ranging from ∼10 μs to ∼100 ms using the weighted ensemble (WE) strategy in combination with GPU computing. Starting from an initial structure or set of structures, WE organizes an ensemble of GPU-accelerated MD trajectory segments via intermittent pruning and replication events to generate statistically unbiased estimates of rate constants for rare events such as folding; no biasing forces are used. Although the variance among atomistic WE folding runs is significant, multiple independent runs are used to reduce and quantify statistical uncertainty. Folding times are estimated directly from WE probability flux and from history-augmented Markov analysis of the WE data. Three systems were examined: NTL9 at low solvent viscosity (yielding τ_f = 0.8-9 μs), NTL9 at water-like viscosity (τ_f = 0.2-2 ms), and Protein G at low viscosity (τ_f = 3-200 ms). In all cases, the folding time, uncertainty, and ensemble properties could be estimated from WE simulation; for Protein G, this characterization required significantly less overall computing than would be required to observe a single folding event with conventional MD simulations. Our results suggest that the use and calibration of force fields and solvent models for precise estimation of kinetic quantities is becoming feasible.

Figure 1:

The WE procedure and comparison to regular MD simulation. (A) A schematic of the WE simulation procedure is shown with two-dimensional binning for protein folding. Three iterations (red, then blue, then green) are shown based on a target number of 4 trajectories per bin, illustrating the “statistical ratcheting” effect which is possible without applying biasing forces. Note that a set of new trajectories is shown only for those parent trajectories that reached a new bin. (B) A brute-force MD simulation of NTL9 leads to the “MD-partly folded” structure (black structure) with a Cα-RMSD of 6.5 Å with respect to the folded crystal structure (blue structures at right) after 7 μs of simulation time. By contrast, a WE simulation starting from the same initial structure (blue structure at left) samples the “WE-folded” NTL9 structure with Cα-RMSD < 1 Å (red structure on the right panel). The WE simulation time is the aggregate time including all trajectory segments, representing a fair comparison using roughly the same amount of computational resources.

Figure 2:

Rate constant estimations for NTL9 folding using 2D WE method with solvent viscosity (γ) set to 5 ps⁻¹. The red lines show the nominal 95% Credibility Region (CR) as a function of molecular time from Bayesian bootstrapping based on direct WE rate constant estimates, which were windowed averages of the previous 1 ns of molecular time for each of the 10 independent simulations (see Figure S1). The green lines show the 95% CR for rate constants obtained by the haMSM method. The experimental rate constant is shown in gold, but note that the low viscosity used in these simulations is expected to yield overly fast kinetics.

Figure 3:

Rate constant estimations for NTL9 folding using 1D WE method with solvent viscosity (γ) set to 80 ps⁻¹. The red lines show the nominal 95% Credibility Region (CR) as a function of molecular time from Bayesian bootstrapping based on direct WE rate constant estimates, which were windowed averages of the previous 1 ns of molecular time for each of the 30 independent simulations (see Figure S2). The green lines show the 95% CR for rate constants obtained by the haMSM method. The experimental rate constant is shown in gold.

Figure 4:

Rate constant estimations for Protein G folding using 2D WE method with solvent viscosity (γ) set to 5 ps⁻¹. The red lines show the nominal 95% Credibility Region (CR) as a function of molecular time from Bayesian bootstrapping based on direct WE rate constant estimates, which were windowed averages of the previous 1 ns of molecular time for each of the 15 independent simulations (see Figure S3). The green lines show the 95% CR for rate constants obtained by the haMSM method. The experimental rate constant is shown in gold, but note that the low viscosity used in these simulations is expected to yield overly fast kinetics.

Figure 5:

A set of example NTL9 (A) and Protein G (B) structures with decreasing Cα-RMSDs from left to right obtained from a continuous trajectory along with the folded crystal structure. Residues are colored based on their native secondary structures in violet (α-helix), green (β-sheet), and cyan (loops). Native backbone hydrogen bonds are indicated as dashed lines, if they emerge in the structure shown.