Equilibration of experimentally determined protein structures for molecular dynamics simulation

Preceding molecular dynamics simulations of biomolecular interactions, the molecule of interest is often equilibrated with respect to an initial configuration. This so-called equilibration stage is required because the input structure is typically not within the equilibrium phase space of the simulation conditions, particularly in systems as complex as proteins, which can lead to artifactual trajectories of protein dynamics. The time at which nonequilibrium effects from the initial configuration are minimized—what we will call the equilibration time—marks the beginning of equilibrium phase-space exploration. Note that the identification of this time does not imply exploration of the entire equilibrium phase space. We have found that current equilibration methodologies contain ambiguities that lead to uncertainty in determining the end of the equilibration stage of the trajectory. This results in equilibration times that are either too long, resulting in wasted computational resources, or too short, resulting in the simulation of molecular trajectories that do not accurately represent the physical system. We outline and demonstrate a protocol for identifying the equilibration time that is based on the physical model of Normal Mode Analysis. We attain the computational efficiency required of large-protein simulations via a stretched exponential approximation that enables an analytically tractable and physically meaningful form of the root-mean-square deviation of atoms comprising the protein. We find that the fitting parameters (which correspond to physical properties of the protein) fluctuate initially but then stabilize for increased simulation time, independently of the simulation duration or sampling frequency. We define the end of the equilibration stage—and thus the equilibration time—as the point in the simulation when these parameters attain constant values. Compared to existing methods, our approach provides the objective identification of the time at which the simulated biomolecule has entered an energetic basin. For the representative protein considered, bovine pancreatic trypsin inhibitor, existing methods indicate a range of $0.2–10\mathrm{ns}$ of simulation until a local minimum is attained. Our approach identifies a substantially narrower range of $4.5–5.5\mathrm{ns}$, which will lead to a much more objective choice of equilibration time.

Figure 1

Outline of the steps in a typical MD protein simulation, illustrated for a steered molecular dynamics simulation of biotin unbinding from streptavidin [10]. The first step is to import an initial structure from the PDB. Next, solvent is added and the structure is energy-minimized via a non-MD algorithm such as steepest-descents minimization to remove close contacts. Then the system is equilibrated: the minimized experimental structure undergoes MD simulation to bring it into equilibrium with the simulation conditions. Only after this equilibration stage can production simulations be run that will yield realistic results. We are proposing a new procedure for the equilibration stage, which is indicated by a shaded background.

Figure 2

Raw trajectory of solvated BPTI at $300\mathrm{K}$ as a function of time (black), showing the corresponding KWW model fit (gray). The correlation coefficient for this fit is 0.843, indicating good agreement between the data and the model.

Figure 3

Bovine pancreatic trypsin inhibitor was chosen because it is a small (58 amino acids) protein, which has been widely studied computationally. Thus, large simulated times could be considered for reasonable real-time durations. This image was produced using the UCSF Chimera package from the Resource for Biocomputing, Visualization, and Informatics at the University of California, San Francisco (supported by NIH P41 RR-01081) [36].

Figure 4

For simulated BPTI at $T=300\mathrm{K}$, parameters of the KWW model of the RMSD referenced to $t_{\mathrm{ref}}=0\mathrm{ns}$ (dashed line) and $t_{\mathrm{ref}}=11\mathrm{ns}$ (solid line) became independent of simulation duration when fitting Eq. (3) to only $5\mathrm{ns}$ of the trajectory ($t=0–5\mathrm{ns}$ or $t=11–16\mathrm{ns}$, respectively). However, the RMSD referenced to $t=2\mathrm{ns}$ (dotted line) required $20\mathrm{ns}$ of the trajectory ($t=2–22\mathrm{ns}$) before the fitting parameters $A_e$, ${\tau }_e$, and $\beta$ stabilized to unique values.

Figure 5

Consideration of BPTI fluctuations via the KWW model. (A) The complexity parameter $\beta$ decreases before increasing to a steady value of about 0.43 at $5\mathrm{ns}$. (B) The exponential prefactor $A_e$ (solid line) and the effective time constant ${\tau }_e$ (dashed line) fluctuated before attaining steady values at $5\mathrm{ns}$. These parameters were obtained at each reference time $t_{\mathrm{ref}}$ by fitting Eq. (3) to the $20\mathrm{ns}$ of the trajectory following each $t_{\mathrm{ref}}$, as this was the maximum sampling duration required for the KWW model of the RMSD to attain a stable, unique functional form.

Figure 6

Identification of an equilibrated BPTI structure via the RMSD plateau method. (A) Examining the trajectory between 0 and $1\mathrm{ns}$ (bottom) gives an equilibration time of $\sim 200\mathrm{ps}$, but examining the trajectory between 0 and $40\mathrm{ns}$ (top) gives an equilibration time of $\sim 10\mathrm{ns}$. (B) Ambiguity arises from examining the trajectory between 0 and $2\mathrm{ns}$ due to a double plateau, one at RMSD $\sim 0.2\mathrm{nm}$ and the other at RMSD $\sim 0.33\mathrm{nm}$ (dashed lines). It is not obvious which, if either, corresponds to the end of nonequilibrium fluctuations.

Figure 7

Identification of an equilibrated BPTI structure via the Stella method of sharp decreases in RMSD plateau values. (A) Here, two different reference structures ($t_{\mathrm{ref}}=1$ and $2\mathrm{ns})$ both give distinct drops in the plateau value, leading to ambiguity in determining the proper equilibration time. Black: RMSD ($t_{\mathrm{ref}}=0\mathrm{ns}$, $t$); dark gray: RMSD ($t_{\mathrm{ref}}=1\mathrm{ns}$, $t$); medium gray: RMSD ($t_{\mathrm{ref}}=2\mathrm{ns}$, $t$). (B) While the RMSD plateau values of all the reference structures are distinct at $t=5\mathrm{ns}$ (left arrow), differences decrease when considering fluctuations to $t=10\mathrm{ns}$ (right arrow). Examining only part of the data shown, e.g., from $0–5\mathrm{ns}$, would lead to choosing an equilibration time from this trajectory, but when longer simulations are considered it becomes unclear as to whether equilibration has occurred according to this method. Black: RMSD ($t_{\mathrm{ref}}=0\mathrm{ns}$, $t$); dark gray: RMSD ($t_{\mathrm{ref}}=1\mathrm{ns}$, $t$); medium gray: RMSD ($t_{\mathrm{ref}}=2\mathrm{ns}$, $t$); light gray: RMSD ($t_{\mathrm{ref}}=3\mathrm{ns}$, $t$).

Figure 8

The simple RMSD-based method and the Stella method both result in a considerable range in potential equilibration times and thus in potential structures for solvated BPTI, reflecting the ambiguity present in those methods. Our KWW-based method gives a comparatively unambiguous equilibration time of $4.5–5.5\mathrm{ns}$.

Figure 9

Consideration of fluctuations in streptavidin alpha-carbon atoms via the KWW model. (A) The complexity parameter $\beta$ fluctuates before attaining to a steady value of about 0.35 at $15\mathrm{ns}$. (B) The exponential prefactor $A_e$ (solid line) and the effective time constant ${\tau }_e$ (dashed line) fluctuated before attaining steady values at $15\mathrm{ns}$. These parameters were obtained at each reference time $t_{\mathrm{ref}}$ by fitting Eq. (3) to the trajectory following each $t_{\mathrm{ref}}$, checking at each $t_{\mathrm{ref}}$ that the parameters did not depend on the duration of trajectory used.

Figure 10

Determining the applicability of the KWW approximation. (A) The RMSD of our computationally modeled protein (BPTI) on a logarithmic time scale (B) RMSD of a single simple harmonic oscillator (solid): KWW model is applicable, with the number of harmonic oscillators equal to one. RMSD of a system with two discrete narrow peaks in the time-constant distribution (dashed): KWW model is not applicable. RMSD of a system with broad overlapping peaks in the time-constant distribution (dotted): KWW model is applicable. The RMSD of our computationally modeled protein BPTI most resembles the dotted line, indicating that the time-constant distribution likely consists of broad, overlapping peaks. Therefore, we conclude that the KWW approximation is applicable to our system.