Stiff-spring approximation revisited: inertial effects in non-equilibrium trajectories

Abstract: Use of harmonic guiding potentials is the most common method for implementing steered molecular dynamics (SMD) simulations, performed to obtain potentials of mean force (PMFs) of molecular systems using non-equilibrium work (NEW) theorems. Harmonic guiding potentials are also the natural choice in single molecule force spectroscopy experiments. The stiff spring approximation (SSA) of Schulten and coworkers enables to use the work performed along SMD trajectories to obtain the PMF. We discuss and demonstrate how a high spring constant, k, required for the validity of the SSA can violate another requirement of this theory, i.e., the validity of Brownian dynamics of the system. Violation of the Brownian condition results in the introduction of kinetic energy contributions to the external work, performed during SMD simulations. These inertial effects result in skewed work distributions, rather than the Gaussian distributions predicted by SSA. The inertial effects also result in broader work distributions, which worsen the effect of the skewness when calculating work averages. Remarkably, our results strongly suggest that the skew and width of work distributions are independent of the average drift velocity and physical asymmetries. The skew and broadening of work distributions result in biased estimation of the PMF. The bias manifests itself in the form of a systematic error that increases with simulation time. We discuss the proper upper limit for k, such that the inertial effects are avoided. This limit, used together with the relation for the lower limit of k, enables to conduct accurate steering while satisfying the Brownian dynamics. Furthermore, we argue and demonstrate that using the peak-value (rather than the statistical mean) of the work distributions vastly reduces the bias in the calculated PMFs and improves the accuracy.

• The paper demonstrates that using overly stiff springs in SMD introduces inertial effects that skew work distributions, biasing PMF calculations.
• It shows that the skewness in work distributions arises primarily from high spring constants, regardless of drift speed and system asymmetries.
• It introduces a peak-finding method that replaces the traditional mean with the peak of work distributions, significantly reducing bias in PMF estimations.

Inertial Effects on PMF Calculations Using the Stiff-Spring Approximation

The stiff-spring approximation (SSA) is widely used in steered molecular dynamics (SMD) simulations to obtain potentials of mean force (PMFs) using Jarzynski's equality and other non-equilibrium work (NEW) theorems. This paper (1607.07430) revisits the SSA, discussing how using a high spring constant, $k$, can violate the Brownian dynamics requirement, leading to inertial effects that skew work distributions and bias PMF estimations. It demonstrates that the skewness and broadening of work distributions are primarily due to using springs that are too stiff, independent of average drift velocity and system asymmetries. Furthermore, the paper introduces a "peak-finding" method that uses the peak value of work distributions instead of the statistical mean to vastly reduce bias in calculated PMFs.

Theoretical Background

Jarzynski's equality relates the free energy difference $\Delta \Phi_{AB}$ between two macrostates $A$ and $B$ to the exponential average of the external work performed on the system, $W_{A\rightarrow B}$:

$$e^{-\beta \Delta \Phi_{AB}} = \langle e^{-\beta W_{A\rightarrow B}} \rangle$$,

where $\beta \equiv (k_B T)^{-1}$. SSA asserts that the distribution of external work will be Gaussian when using a sufficiently stiff spring. Kosztin's forward-reverse (FR) method provides a simpler mechanism for understanding work distributions and obtaining PMFs. The FR method exploits the equality of average dissipative work in forward and reverse directions:

$$\Delta \Phi_{AB} = \frac{\langle W \rangle_F - \langle W \rangle_R}{2}$$.

Inertial Effects and Skewed Work Distributions

The paper demonstrates that a high spring constant $k$ can lead to significant deviations from Brownian dynamics, introducing kinetic energy contributions to the external work. This results in skewed work distributions, rather than the Gaussian distributions predicted by SSA. The inertial effects are more pronounced in larger biomolecular simulations, where longer sampling is required for convergence, leading to a systematic error that increases with simulation time.

Figure 1: Reaction path for an HHC-36 peptide molecule (sequence KRWWKWWRR) as it is steered to go towards the surface of a patch of mixed (3:1) POPE/POPG membrane.

Figure 1 shows the reaction path for steering an HHC-36 peptide toward a POPE/POPG membrane. The work distributions obtained from SMD simulations with different spring constants vividly exhibit the skewness and broadening caused by using stiffer springs. The skewness values ($\alpha$) are calculated as:

$\alpha = \frac{\mu_3}{\sigma^3}$,

where $\mu_3$ is the third moment around the mean, and $\sigma$ is the standard deviation of work.

(Figure 2)

Figure 2: Work distributions obtained from SMD simulations with an HHC-36 peptide steered to go toward a POPE/POPG membrane patch, using three different steering spring constants.

Figure 2 shows that increasing the spring constant leads to right-skewed distributions, indicating that inertial effects are present. This figure shows work distributions obtained from simulations with different spring constants, showing the skewness and broadening of the distributions with stiffer springs. The authors also introduce the notion of scaled work and the zooming factor to better visualize these distributions.

Dependence of Skewness on Physical Parameters

The paper investigates the influence of physical asymmetries and drift speed on work distributions, using a system consisting of a $\text{C}_{60}$ buckminsterfullerene molecule approaching a silicon dioxide ($\text{SiO}_2$) slab in an aqueous environment.

(Figure 3)

Figure 3: Depiction of the reaction path for a $\text{C}_{60}$ buckminsterfullerene molecule, as it is restrained to go toward a silicon dioxide ($\text{SiO}_2$) slab and back.

Figure 3 depicts the reaction path for the $\text{C}_{60}$-silica system, showing the setup for studying the effects of steering speed and spring constant on work distributions.

(Figure 4)

Figure 4: Forward work distributions ($\text{C}_{60}$ going toward the silica slab) in the bin with $z\in[20,25]$, for the system shown in Fig. 3.

Figure 4 shows that variations in drift speed have minimal impact on the width or skewness of the distributions, whereas changes in the spring constant dramatically alter the width. The mass of the steered objects also plays a role, as demonstrated by simulations with modified ion masses, where increasing the mass significantly reduces skewness. The skewness appears only when both steered objects are restrained, and is less pronounced when one object is constrained.

Biased PMFs due to Inertial Effects

The paper demonstrates that inertial effects can lead to biased PMF estimations. Using longer simulation times with stiffer springs exacerbates the bias, indicating a systematic error in the free energy profile.

(Figure 5)

Figure 5: K$^+$-Na$^+$ work distributions for the reaction path bin with $z\in[4.0, 5.0]$ from seven series of simulations with different steering spring constants and masses.

Figure 5 shows K$^+$-Na$^+$ work distributions with varying spring constants and masses, highlighting the impact of these parameters on the skewness of the distributions. The impact of skewness of the work distributions from a SMD simulation can be understood by examining the PMFs in Figure 6.

(Figure 6)

Figure 6: PMFs obtained from the same set of SMD simulations of a peptide-membrane system, as described under Fig. 2, at the far region ($z\in[45,55]$) where there is essentially no interaction and thus a flat PMF is expected.

Figure 6 shows that using larger spring constants yields sloped (biased) PMFs, while softer springs align with the equilibrium PMF. These biased PMFs can be attributed to the skewness of the work distributions, where deviations from the correct average forward and reverse works are unequal for each bin.

(Figure 7)

Figure 7: PMFs obtained from the same set of SMD simulations of a peptide-membrane system, as described under Fig. 2.

Figure 7 further illustrates the biased PMFs obtained from simulations with larger spring constants in the peptide-membrane system.

Fine-Tuning the Stiff Spring Approximation

The paper provides criteria for selecting the appropriate spring constant, balancing the need for precise steering with the requirement for Brownian dynamics. The upper limit for $k$ is determined by:

$k < \frac{9 \pi^2 \eta^2 r^2}{m}$,

where $\eta$ is the shear viscosity of the fluid, $r$ is the radius of the steered object, and $m$ is its mass. The lower limit for $k$ is:

$k > \frac{k_B T}{(\delta x)^2}$,

where $\delta x$ is the maximum tolerable deviation from the prescribed trajectory.

The Peak-Finding Method

To mitigate the bias introduced by inertial effects, the paper introduces a "peak-finding" method that uses the peak value of work distributions instead of the statistical mean in the FR formula. The peak location serves as a better estimate of the unbiased $\langle W \rangle$, reducing the impact of skewed distributions. This method involves establishing work distributions for each bin and accurately determining the peak location using a linear curve-fitting procedure.

(Figure 8)

Figure 8: PMFs for the interaction of the HHC-36 peptide with a POPE/POPG membrane patch, obtained using the peak-finding method.

Figure 8 demonstrates the improvement achieved by using the peak-finding method, resulting in less biased PMFs, even with high $k$ values.

(Figure 9)

Figure 9: K$^+$-Na$^+$ PMFs calculated using bin-passing (a and c) and peak-finding (b and d) methods, for systems where the ion masses are left at their natural values (a and b) and modified (c and d).

Figure 9 shows K$^+$-Na$^+$ PMFs calculated using both bin-passing and peak-finding methods, demonstrating the effectiveness of the peak-finding approach in reducing bias and improving convergence. The software for building work distributions and finding peak locations is publicly available.

Conclusion

This paper (1607.07430) highlights the importance of carefully selecting the spring constant in SMD simulations to avoid inertial effects that can skew work distributions and bias PMF estimations. It demonstrates that the skewness of work distributions is primarily due to using springs that are too stiff. The introduced "peak-finding" method offers a robust approach for reducing bias and improving the accuracy of PMF calculations, particularly in large biomolecular systems.