Path probability ratios for Langevin dynamics—Exact and approximate

Kieninger, S.; Keller, B. G.

doi:10.1063/5.0038408

Path reweighting is a principally exact method to estimate dynamic properties from biased simulations—provided that the path probability ratio matches the stochastic integrator used in the simulation. Previously reported path probability ratios match the Euler–Maruyama scheme for overdamped Langevin dynamics. Since molecular dynamics simulations use Langevin dynamics rather than overdamped Langevin dynamics, this severely impedes the application of path reweighting methods. Here, we derive the path probability ratio M_L for Langevin dynamics propagated by a variant of the Langevin Leapfrog integrator. This new path probability ratio allows for exact reweighting of Langevin dynamics propagated by this integrator. We also show that a previously derived approximate path probability ratio M_approx differs from the exact M_L only by $O (ξ^{4} Δ t^{4})$ and thus yields highly accurate dynamic reweighting results. (Δt is the integration time step, and ξ is the collision rate.) The results are tested, and the efficiency of path reweighting is explored using butane as an example.

I. INTRODUCTION

Molecular dynamics are astonishingly complex and occur in a wide range of length and timescales.^1–3 To elucidate the mechanisms by which different parts of a molecular system interact and how macroscopic properties arise from these interactions, molecular dynamics (MD) simulations have become an indispensable tool.^4–9 Because the timescales covered by MD simulations are often orders of magnitude lower than the slowest timescale of the system, a wide variety of enhanced sampling techniques have been developed, which distort the dynamics of the simulation such that rare molecular transitions occur more frequently. This can be achieved by raising the temperature or by adding a bias to the potential energy function.^10,11 How to extract the correct values of dynamical properties (mean-first passage times, residence times, binding rates, or transition probabilities) from these accelerated dynamics is an open question and a very active field of research.

The goal of dynamical reweighting methods is to estimate dynamical properties of the system at a target state $\tilde{S}$ from a trajectory generated at simulation state S. S could correspond to a higher temperature or to a biased potential. Starting points for the derivation of dynamical reweighting methods are Kramers rate theory,^12–15 the likelihood function for estimating the transition probabilities from MD trajectories,^16–19 or a discretization of the Fokker–Planck equation.^7,20–22 The methods differ in the ease of use and the severity of the assumptions they make.²³

A principally exact formalism to reweight dynamic properties is path reweighting methods, which have been reported already early in Refs. 24–28. In path reweighting methods, the trajectory generated at state S is split into short paths ω. Then, the path probability ${\tilde{P}}_{L} (ω; Δ t | (x_{0}, v_{0}))$ of a given ω at the target state $\tilde{S}$ is calculated by reweighting the path probability P_L(ω; Δt|(x₀, $v_{0}$ ⁠)) of ω at the simulation state S,

{\tilde{P}}_{L} (ω; Δ t | (x_{0}, v_{0})) \approx M \cdot P_{L} (ω; Δ t | (x_{0}, v_{0})) .

(1)

(x₀, $v_{0}$ ⁠) is the initial state of the path ω, and Δt is the integration time step. M(ω) is the path probability ratio or reweighting factor. Equation (1) is exact if the path probability ratio $M = {\tilde{P}}_{L} (ω; Δ t | (x_{0}, v_{0})) / P_{L} (ω; Δ t | (x_{0}, v_{0}))$ is derived from the numerical integration scheme used to generate ω. The mathematical basis for path reweighting methods is the Girsanov theorem,^29,30 or else, they can be derived from the Onsager–Machlup action.^24–27,31 A prerequisite for path reweighting is that a stochastic integrator is used in the MD simulation, e.g., a Langevin thermostat.

However, it has been challenging to apply path reweighting to simulations of large molecular systems. For example, the variance of the reweighting estimators increases rapidly with increasing path length such that for long paths, reweighting becomes inefficient compared to direct simulation of the target state. Combining path reweighting techniques with Markov state models (MSMs) alleviates this problem.^32–35 In MSMs,^36–42 the dynamics of the system is represented by transitions between discrete states in the conformational space of the molecular system, where the lag time τ of the transition is much shorter than the slow timescales of the system. Thus, only short paths of length τ are needed to estimate and reweight the transition probabilities.

Second, a number of technical difficulties arise. The path probability ratio M decreases exponentially with the path length τ such that the standard numerical accuracy is quickly exceeded. This problem can be solved by using high precision arithmetic libraries.³⁵ To calculate the path probability ratio M, one needs to know the trajectory and the random numbers of the stochastic integrator at every integration time step. Writing this information to disk at every integration time step is not a workable option. We, therefore, proposed to calculate the path reweighting factor “on-the-fly” during the simulation and to write out intermediate results at regular intervals, e.g., whenever the positions are written to disk. The additional storage requirements and computational costs for the “on-the-fly”-calculations are negligible compared to the overall cost of the simulation.^34,35 Having solved the technical challenges, we tested the path reweighting method on several peptides using path lengths of up to τ = 600 ps.^34,35 Applications to larger systems and longer path lengths are likely within reach.

Yet, the equation for the path probability ratio M poses a barrier to a more widespread use of path reweighting methods. Because M is derived from the stochastic integration scheme used to simulate the system, one cannot readily apply a path probability ratio derived for one integration scheme to a simulation generated by another integration scheme.

In temperature reweighting, i.e., when simulation and target state differ in the temperature, only the random term of the stochastic integrator is affected by the change in temperature. Path probability ratios for temperature reweighting have been constructed by rescaling the normal distributions of the random or noise terms of the stochastic integration scheme.^32,43

In potential reweighting, i.e., when simulation and target state differ in the potential energy function, one needs to account for changes in the drift terms of the stochastic integration scheme. The path probability ratio M_o for the Euler–Maruyama scheme for overdamped Langevin dynamics has been reported multiple times.^24–26,33 However, the dynamics of large molecular systems is better reproduced by Langevin dynamics, and MD programs implement a wide variety of Langevin integration schemes.^44–53 The time-continuous Onsager–Machlup action for Langevin dynamics has been reported,²⁷ but to the best of our knowledge, path probability ratios for Langevin integration schemes M_L have not yet been reported. Thus, exact path reweighting for Langevin dynamics has not been possible, so far.

In Refs. 34 and 35, we demonstrated that path reweighting can be applied to biased simulations of large molecular systems, nonetheless. We used an approximate path probability ratio M_approx that is based on the path probability ratio for the Euler–Maruyama scheme but uses the random numbers that are generated during the Langevin MD simulation. We tested M_approx extensively, and for low-dimensional model systems and for molecular systems, this approximate path probability ratio yielded very accurate results. In these two publications, we used a variant of the Langevin Leapfrog integration scheme developed by Izaguirre, Sweet, and Pande⁴⁹ to propagate the system. Both the Langevin Leapfrog integration scheme and its variant are implemented in OpenMM⁵⁴ (see Appendix A). We will abbreviate the variant by the “ISP scheme.”

In this contribution, we derive the path probability ratio M_L for Langevin dynamics propagated by a variant of the Langevin Leapfrog integrator.⁴⁹ M_L allows for exact reweighting of Langevin dynamics (Sec. IV). We analyze why M_approx is an excellent approximation to M_L (Sec. VI), and we discuss whether there are scenarios in which M_o is a viable approximation to M_L (Sec. V). The general framework of the path reweighting equations and the corresponding equations for the Euler–Maruyama scheme are summarized in Secs. II and III. Section VIII reports the computational details.

II. PATH REWEIGHTING

The path probability P(ω; Δt|(x₀, $v_{0}$ ⁠)) is the probability to generate a time-discretized path ω = (x₀, x₁, …, x_n) starting in a pre-defined initial state (x₀, $v_{0}$ ⁠) at the simulation potential V(x). The notation emphasizes that the probability is conditioned on an initial state (x₀, $v_{0}$ ⁠) and that the path has been generated with a fixed time step Δt, whereas ω is the argument of the function. In short, P(ω; Δt|(x₀, $v_{0}$ ⁠)) maps a path in position space to a probability. Its functional form depends on the integration scheme used to generate ω and the potential energy function.

The path probability ratio is the ratio between the probability $\tilde{P} (ω; Δ t | (x_{0}, v_{0}))$ to generate a path ω at a target potential,

\tilde{V} (x) = V (x) + U (x),

(2)

and the probability P(ω; Δt|(x₀, $v_{0}$ ⁠)) to generate the same path ω at the simulation potential V(x),

M (ω; Δ t | (x_{0}, v_{0})) = \frac{\tilde{P} (ω; Δ t | (x_{0}, v_{0}))}{P (ω; Δ t | (x_{0}, v_{0}))} .

(3)

The potential energy function U(x) is usually called perturbation or bias.

In integration schemes for stochastic dynamics, random numbers are used to propagate the system. If a single random number is drawn per integration step, then the probability to generate ω is equal to the probability P(η) to generate the corresponding random number sequence η = (η₀, η₁, …, η_n−1),

P (ω; Δ t | (x_{0}, v_{0})) = P (η),

(4)

where ω and η are linked by the equations for the integration scheme. Since the random numbers η_k are drawn from a Gaussian distribution with zero mean and unit variance, the functional form of P(η) is

P (η) = N \exp (- \frac{1}{2} \sum_{k = 0}^{n - 1} η_{k}^{2}), N = {(\frac{1}{2 π})}^{\frac{n}{2}} .

(5)

P(η) is a function that maps a random number sequence to a probability. One can interpret Eq. (4) as a change in variables from ω to η, where the change is defined by the equations for the integration scheme.

Suppose that η is the random number sequence needed to generate ω at a simulation potential V(x). To generate the same path at a target potential $\tilde{V} (x)$ ⁠, one would need a different random number sequence $\tilde{η} = ({\tilde{η}}_{0}, {\tilde{η}}_{1}, \dots, {\tilde{η}}_{n - 1})$ ⁠, with

{\tilde{η}}_{k} = η_{k} + Δ η_{k} .

(6)

Δη_k is the random number difference, and it depends on the integration scheme and the difference between the two potentials. The random number probability ratio is the ratio between the probability of drawing η and the probability of drawing ${\tilde{η}}_{k}$ ⁠,

\begin{matrix} \frac{P (\tilde{η})}{P (η)} & = & \frac{N \exp (- \frac{1}{2} \sum_{k = 0}^{n - 1} {(η_{k} + Δ η_{k})}^{2})}{N \exp (- \frac{1}{2} \sum_{k = 0}^{n - 1} η_{k}^{2})} \\ = & \exp (- \sum_{k = 0}^{n - 1} η_{k} \cdot Δ η_{k}) \cdot \exp (- \frac{1}{2} \sum_{k = 0}^{n - 1} {(Δ η_{k})}^{2}) . \end{matrix}

(7)

Mathematically, the following has happened in the previous paragraph. The path ω remained unchanged. The functional form of the path probability has changed as $\tilde{P} (ω; Δ t | (x_{0}, v_{0}))$ because the potential energy enters the equations for the integration scheme. Likewise, the change in variables from ω to $\tilde{η}$ has changed. The functional form of the random number probability remains the same [Eq. (5)]. The analogon to Eq. (4) for the target potential is

\tilde{P} (ω; Δ t | (x_{0}, v_{0})) = P (\tilde{η}),

(8)

where ω and $\tilde{η}$ are linked by the equations for the integration scheme using $\tilde{V} (x)$ ⁠. Given the two changes in variables for the simulation and the target potential, the path probability ratio [Eq. (3)] and the random number probability ratio [Eq. (7)] are equal. Note that Eq. (3) is a ratio of two different functions that have the same argument ω, whereas Eq. (7) is the ratio of the same function with different arguments η and $\tilde{η}$ ⁠.

Equation (7) is of little practical use because $\tilde{η}$ is not available from a simulation at the simulation state. However, the random number difference Δη_k can be expressed as a function of ω, and the random number probability ratio can thus be expressed as a function of ω and η,

M (ω, η; Δ t | (x_{0}, v_{0})) = \frac{P (\tilde{η})}{P (η)} .

(9)

For a path ω and the corresponding random number sequence η that was used to generate this path, we will use the following equality:

M (ω, η; Δ t | (x_{0}, v_{0})) = M (ω; Δ t | (x_{0}, v_{0})) .

(10)

The functional form and the value of the properties introduced in this section depend strongly on the integration scheme. In Sec. III, we summarize the equations for the Euler–Maruyama scheme for overdamped Langevin dynamics. In Sec. IV, we derive the corresponding equations for the ISP integration scheme for Langevin dynamics (see Table I). Throughout this manuscript, properties associated with Langevin dynamics are subscripted with L, and properties associated with overdamped Langevin dynamics are subscripted with o.

TABLE I.

References to the equations for the properties introduced in Sec. II.

		Overdamped Langevin	Langevin
Equation of motion		Eq. (11)	Eq. (19)
Integration scheme		Eq. (12)	Eqs. (20) and (21)
Path probability	P(ω; Δt\|(x₀, $v_{0}$ ⁠))	Eq. (13)	Eq. (22)
Path probability ratio	M(ω; Δt\|(x₀, $v_{0}$ ⁠))	Eq. (14)	Eq. (23)
Random number	η_k	Eq. (15)	Eq. (24)
Random number difference	Δη_k	Eq. (17)	Eq. (26)
Random number probability ratio	M(ω, η; Δt\|(x₀, $v_{0}$ ⁠))	Eq. (18)	Eq. (27)

		Overdamped Langevin	Langevin
Equation of motion		Eq. (11)	Eq. (19)
Integration scheme		Eq. (12)	Eqs. (20) and (21)
Path probability	P(ω; Δt\|(x₀, $v_{0}$ ⁠))	Eq. (13)	Eq. (22)
Path probability ratio	M(ω; Δt\|(x₀, $v_{0}$ ⁠))	Eq. (14)	Eq. (23)
Random number	η_k	Eq. (15)	Eq. (24)
Random number difference	Δη_k	Eq. (17)	Eq. (26)
Random number probability ratio	M(ω, η; Δt\|(x₀, $v_{0}$ ⁠))	Eq. (18)	Eq. (27)

III. OVERDAMPED LANGEVIN DYNAMICS

A. Equation of motion and integration scheme

Consider a one particle system that moves in a one-dimensional position space with temperature T and potential energy function V. The overdamped Langevin equation of motion is

\dot{x} (t) = - \frac{\nabla V (x (t))}{ξ m} + \sqrt{\frac{2 k_{B} T}{ξ m}} η (t),

(11)

with particle mass m, position x, velocity $v = \dot{x}$ ⁠, and Boltzmann constant k_B. x(t) ∈ Ω_o is the state of the system at time t, where $Ω_{o} \subset R$ is the state space of the system. The collision rate ξ (in units of s⁻¹) models the interaction with the thermal bath. $η (t) \in R$ describes an uncorrelated Gaussian white noise with unit variance centered at zero, which is scaled by the volatility $\sqrt{\frac{2 k_{B} T}{ξ m}}$ ⁠.

A numerical algorithm to calculate an approximate solution to Eq. (11) is the Euler–Maruyama integration scheme,^30,55

x_{k + 1} = x_{k} - \frac{\nabla V (x_{k})}{ξ m} Δ t + \sqrt{\frac{2 k_{B} T}{ξ m}} \sqrt{Δ t} η_{o, k},

(12)

where Δt is the time step, x_k is the position, and η_o,k is the random number at iteration k. The random numbers are drawn from a Gaussian distribution with zero mean and unit variance. For k = 0, …, n − 1, Eq. (12) yields a time-discretized overdamped Langevin path ω_o = (x₀, x₁, …, x_n), which starts at the pre-defined initial position x₀. Note that while the state of the system at iteration k is defined by the position x_k, the progress to x_k+1 depends on x_k and on the value of the random number η_o,k. The random number sequence that was used to generate a specific ω_o is denoted by η_o = (η_o,0, …, η_o,n−1).

B. Path probability and path probability ratio

The probability to observe a path ω_o generated by the Euler–Maruyama scheme [Eq. (12)] is^28,34,56,57

\begin{align} P_{o} (ω_{o}; Δ t | x_{0}) \\ = {[\sqrt{\frac{ξ m}{4 π k_{B} T Δ t}}]}^{n} \\ \cdot \exp (- \frac{ξ m}{4 k_{B} T Δ t} \sum_{k = 0}^{n - 1} {(x_{k + 1} - x_{k} + \frac{Δ t}{ξ m} \nabla V (x_{k}))}^{2}) . \end{align}

(13)

For the Euler–Maruyama scheme, the path probability P_o(ω_o; Δt|x₀) does not depend on the initial velocity; hence, we dropped $v_{0}$ in the notation. However, it does depend on the potential energy function V(x) that has been used in Eq. (12) to generate the path ω_o.

The path probability that the same path ω_o has been generated at a target potential $\tilde{V} (x)$ [Eq. (2)] is ${\tilde{P}}_{o} (ω_{o}; Δ t | x_{0})$ ⁠, which is obtained by replacing the potential V(x) with $\tilde{V} (x)$ in Eq. (13). The ratio between the two path probabilities is

\begin{matrix} M_{o} (ω_{o}; Δ t | x_{0}) \\ = \frac{{\tilde{P}}_{o} (ω_{o}; Δ t | x_{0})}{P_{o} (ω_{o}; Δ t | x_{0})} \\ = \exp (- \frac{\sum_{k = 0}^{n - 1} (x_{k + 1} - x_{k}) (\nabla \tilde{V} (x_{k}) - \nabla V (x_{k}))}{2 k_{B} T}) \\ \times \exp (- \frac{\sum_{k = 0}^{n - 1} (\nabla {\tilde{V}}^{2} (x_{k}) - \nabla V^{2} (x_{k})) Δ t}{4 k_{B} T ξ m}) . \end{matrix}

(14)

Equation (14) is a function of the path ω_o and does not depend on the random number sequence η_o explicitly. It is equivalent to Eq. (B4) in Ref. 34.

C. Random numbers and random number probability ratio

Given ω_o, the sequence of random numbers η_o that was used to generate ω_o at the simulation potential V(x) can be back-calculated by rearranging Eq. (12) for η_o,k,

η_{o, k} = \sqrt{\frac{ξ m}{2 k_{B} T Δ t}} (x_{k + 1} - x_{k} + \frac{\nabla V (x_{k})}{ξ m} Δ t) .

(15)

We remark that the path probability [Eq. (13)] can formally be derived by inserting Eq. (15) into Eq. (5). Since Eq. (15) defines a coordinate transformation from x_k to η_o,k, one needs to normalize with respect to the new coordinates in order to obtain the correct normalization constant. The random number sequence ${\tilde{η}}_{o}$ needed to generate ω_o at a target potential $\tilde{V} (x)$ is calculated by inserting Eq. (2) into Eq. (15),

\begin{matrix} {\tilde{η}}_{o, k} & = & \sqrt{\frac{ξ m}{2 k_{B} T Δ t}} (x_{k + 1} - x_{k} + \frac{\nabla V (x_{k})}{ξ m} Δ t) + \sqrt{\frac{Δ t}{2 k_{B} T ξ m}} \nabla U (x_{k}) \\ = & η_{o, k} + Δ η_{o, k} . \end{matrix}

(16)

Equation (15) defines the change in variables from ω to η_o for the Euler–Maruyama scheme at the simulation potential. Likewise, Eq. (16) defines the change in variables from ω to ${\tilde{η}}_{o}$ at the target potential. The random number difference is

Δ η_{o, k} = \sqrt{\frac{Δ t}{2 k_{B} T ξ m}} \nabla U (x_{k}) .

(17)

It depends on the perturbation U(x), but not on the simulation potential V(x). Inserting Δη_o,k [Eq. (17)] into Eq. (7) yields the random number probability ratio,

\begin{matrix} M_{o} (ω_{o}, η_{o}; Δ t | x_{0}) \\ = \exp (- \sum_{k = 0}^{n - 1} \sqrt{\frac{Δ t}{2 k_{B} T ξ m}} \nabla U (x_{k}) \cdot η_{o, k}) \\ \cdot \exp (- \frac{1}{2} \sum_{k = 0}^{n - 1} \frac{Δ t}{2 k_{B} T ξ m} {(\nabla U (x_{k}))}^{2}) . \end{matrix}

(18)

Because of Eq. (10), Eqs. (14) and (18) are equal. However, the two probability ratios use different time-series and different information on the system to evaluate the path probability ratio. To evaluate Eq. (14), one needs the path ω_o, the simulation potential V(x), and the target potential $\tilde{V} (x)$ ⁠. To evaluate Eq. (18), one needs the path ω_o, the random number sequence for the simulation potential η_o, and the perturbation U(x). Because U(x) often only affects a few coordinates of the systems, i.e., it is low-dimensional, Eq. (18) is computationally more efficient. Besides the force calculation −∇V(x) needed to generate the path ω_o, it requires an additional force calculation −∇U(x) only along the coordinates that are affected by the perturbation. By contrast, Eq. (14) requires an additional force calculation on the entire system $- \nabla \tilde{V} (x)$ ⁠.

IV. LANGEVIN DYNAMICS

A. Equation of motion and integration scheme

Consider a one particle system that moves in a one-dimensional position space with temperature T and potential energy function V. The Langevin equation of motion is

m \ddot{x} (t) = - \nabla V (x (t)) - ξ m \dot{x} (t) + \sqrt{2 k_{B} T ξ m} η (t),

(19)

with particle mass m, position x, velocity $v = \dot{x}$ ⁠, acceleration $a = \ddot{x}$ ⁠, and Boltzmann constant k_B. The state of the system at time t is determined by the position and the velocity $(x (t), \dot{x} (t)) \in Ω_{L}$ ⁠, where $Ω_{L} \subset R^{2}$ is the state space of the system. The collision rate ξ (in units of s⁻¹) models the interaction with the thermal bath. $η \in R$ describes an uncorrelated Gaussian white noise with unit variance centered at 0, which is scaled by the volatility $\sqrt{2 k_{B} T ξ m}$ ⁠.

A numerical algorithm to calculate an approximate solution to Eq. (19) is the ISP scheme,⁴⁹

\begin{matrix} x_{k + 1} & = & x_{k} + \exp (- ξ Δ t) v_{k} Δ t - [1 - \exp (- ξ Δ t)] \frac{\nabla V (x_{k})}{ξ m} Δ t \\ + \sqrt{\frac{k_{B} T}{m} [1 - \exp (- 2 ξ Δ t)]} η_{L, k} Δ t, \end{matrix}

(20)

v_{k + 1} = \frac{x_{k + 1} - x_{k}}{Δ t},

(21)

where Δt is the time step, x_k is the position, $v$ _k is the velocity, and η_L,k is the random number at iteration k (see Appendix A). The random numbers are drawn from a Gaussian distribution with zero mean and unit variance. For k = 0, …, n − 1, Eqs. (20) and (21) yield a time-discretized Langevin path ω_L = ((x₀, $v_{0}$ ⁠), (x₁, $v_{1}$ ⁠), …, (x_n, $v_{n}$ ⁠)), which starts at the pre-defined initial state (x₀, $v_{0}$ ⁠). Note that while the state of the system at iteration k is defined by the tuple (x_k, $v$ _k) ∈ Ω_L, the progress to (x_k+1, $v$ _k+1) depends on (x_k, $v$ _k) and on the value of the random number η_L,k. The random number sequence that was used to generate a specific ω_L is denoted by η_L = (η_L,0, …, η_L,n−1).

The position x_k+1 is treated as a random variable because it directly depends on a random number [Eq. (20)], while the velocity $v$ _k+1 is calculated from the new position x_k+1 and the preceding position x_k. Because the velocity $v$ _k in Eq. (20) is determined by the positions x_k and x_k−1 [Eq. (21)], it carries a small memory effect into the time-evolution of x.

B. Path probability and path probability ratio

The probability to generate a path ω_L by the ISP scheme [Eqs. (20) and (21)] at the simulation potential V(x) is

\begin{matrix} P_{L} (ω_{L}; Δ t | (x_{0}, v_{0})) & = & [\prod_{k = 0}^{n - 1} δ (v_{k + 1} - \frac{x_{k + 1} - x_{k}}{Δ t})] \cdot {[\sqrt{\frac{m}{2 π k_{B} T Δ t^{2} (1 - \exp (- 2 ξ Δ t))}}]}^{n} \\ \times \exp (- \sum_{k = 0}^{n - 1} \frac{m {(x_{k + 1} - x_{k} - \exp (- ξ Δ t) v_{k} Δ t + (1 - \exp (- ξ Δ t)) \frac{\nabla V (x_{k})}{ξ m} Δ t)}^{2}}{2 k_{B} T (1 - \exp (- 2 ξ Δ t)) Δ t^{2}}) . \end{matrix}

(22)

The derivation of Eq. (22) is shown in Appendixes B and C. Appendix B explains the strategy for the derivation, and Appendix C shows how to solve the integrals that appear in the derivation.

The path probability ${\tilde{P}}_{L} (ω_{L}; Δ t | (x_{0}, v_{0}))$ to generate a path ω_L by the ISP scheme at the target potential is obtained by inserting $\tilde{V} (x)$ [Eq. (2)] into Eq. (22). The path probability ratio for overdamped Langevin dynamics is

\begin{matrix} M_{L} (ω_{L}; Δ t | (x_{0}, v_{0})) \\ = \frac{{\tilde{P}}_{L} (ω_{L}; Δ t | (x_{0}, v_{0}))}{P_{L} (ω_{L}; Δ t | (x_{0}, v_{0}))} \\ = \exp (- \frac{\sum_{k = 0}^{n - 1} (x_{k + 1} - x_{k}) (\nabla \tilde{V} (x_{k}) - \nabla V (x_{k}))}{k_{B} T ξ (1 + \exp (- ξ Δ t)) Δ t}) \\ \cdot \exp (\frac{\sum_{k = 0}^{n - 1} v_{k} (\nabla \tilde{V} (x_{k}) - \nabla V (x_{k}))}{k_{B} T ξ (1 + \exp (ξ Δ t))}) \\ \cdot \exp (- \frac{\exp (ξ Δ t) - 1}{\exp (ξ Δ t) + 1} \cdot \frac{\sum_{k = 0}^{n - 1} (\nabla {\tilde{V}}^{2} (x_{k}) - \nabla V^{2} (x_{k}))}{2 k_{B} T ξ^{2} m}) . \end{matrix}

(23)

Analogous to Eq. (14), Eq. (23) is a function of the path ω_L and does not depend on the random number sequence η_L.

C. Random numbers and random number probability ratio

Given ω_L, the sequence of random numbers η_L, which was used to generate ω_L at the simulation potential V(x), can be back-calculated by rearranging Eq. (20) for η_L,k,

\begin{matrix} η_{L, k} & = & \sqrt{\frac{m}{k_{B} T (1 - \exp (- 2 ξ Δ t)) Δ t^{2}}} \\ \times (x_{k + 1} - x_{k} - \exp (- ξ Δ t) v_{k} Δ t) \\ (+ (1 - \exp (- ξ Δ t)) \frac{\nabla V (x_{k})}{ξ m} Δ t) . \end{matrix}

(24)

The random number sequence ${\tilde{η}}_{L}$ needed to generate ω_L at a target potential $\tilde{V} (x)$ is calculated by inserting Eq. (2) into Eq. (24),

\begin{matrix} {\tilde{η}}_{L, k} & = & \sqrt{\frac{m}{k_{B} T (1 - \exp (- 2 ξ Δ t)) Δ t^{2}}} \\ \times (x_{k + 1} - x_{k} - \exp (- ξ Δ t) (x_{k} - x_{k - 1})) \\ (+ (1 - \exp (- ξ Δ t)) \frac{\nabla V (x_{k})}{ξ m} Δ t) \\ + \sqrt{\frac{1}{k_{B} T ξ^{2} m}} \cdot \frac{1 - \exp (- ξ Δ t)}{\sqrt{1 - \exp (- 2 ξ Δ t)}} \nabla U (x_{k}) \\ = & η_{L, k} + Δ η_{L, k} . \end{matrix}

(25)

Equation (24) defines the change in variables from ω to η_L for the ISP scheme at the simulation potential. Likewise, Eq. (25) defines the change in variables from ω to ${\tilde{η}}_{L}$ at the target potential. The random number difference is

Δ η_{L, k} = \sqrt{\frac{1}{k_{B} T ξ^{2} m}} \cdot \frac{1 - \exp (- ξ Δ t)}{\sqrt{1 - \exp (- 2 ξ Δ t)}} \nabla U (x_{k}) .

(26)

Again, the random number difference depends on the perturbation potential U(x), but not on the simulation potential V(x). Inserting Δη_L,k [Eq. (26)] into Eq. (7) yields the random number probability ratio,

\begin{matrix} M_{L} (ω_{L}, η_{L}; Δ t | (x_{0}, v_{0})) \\ = \exp (- \frac{1 - \exp (- ξ Δ t)}{\sqrt{1 - \exp (- 2 ξ Δ t)}} \cdot \frac{\sum_{k = 0}^{n - 1} \nabla U (x_{k}) η_{L, k}}{\sqrt{k_{B} T ξ^{2} m}}) \\ \times \exp (- \frac{{(1 - \exp (- ξ Δ t))}^{2}}{1 - \exp (- 2 ξ Δ t)} \cdot \frac{\sum_{k = 0}^{n - 1} \nabla U^{2} (x_{k})}{2 k_{B} T ξ^{2} m}) . \end{matrix}

(27)

Analogous to the path probability ratio for overdamped Langevin dynamics, M_L(ω_L; Δt|(x₀, $v$ _o)) [Eq. (23)] and M_L(ω_L, η_L; Δt|(x₀, $v_{0}$ ⁠)) [Eq. (27)] yield the same path probability ratio for a given path ω_L that has been generated using the random number sequence η_L, but they use different arguments. Again, the path probability from random numbers M_L(ω_L, η_L; Δt|(x₀, $v_{0}$ ⁠)) requires an additional force calculation −∇U(x) only along the coordinates that are affected by the perturbation, making it computationally more efficient than M_L(ω_L; Δt|(x₀, $v_{0}$ ⁠)) in most cases.

V. COMPARING LANGEVIN AND OVERDAMPED LANGEVIN DYNAMICS

A. Test system

Our test system is a one-dimensional one particle system at the simulation potential V(x) (Fig. 1, orange line) and at the target potential $\tilde{V} (x)$ (Fig. 1, black line). The trajectories generated at V(x) will be reweighted to the target potential $\tilde{V} (x)$ ⁠. The black lines in Fig. 4(b) represent the first three dominant MSM eigenfunctions⁴⁰ associated with the target potential. The implied timescales³⁷ are t₀ = ∞, t₁ = 20.5 s, and t₂ = 6.0 s, which are shown as black lines in Fig. 4(c). Computational details are reported in Sec. VIII.

FIG. 1.

View large Download slide

Simulation potential V(x) (orange) and target potential $\tilde{V} (x)$ (black).

B. From random numbers η to paths ω_o and ω_L

Given a random number sequence η = (η₀, …, η_n−1) and a starting state (x₀, $v_{0}$ ⁠), one can use the Euler–Maruyama scheme to generate an overdamped Langevin path ω_o, or else, one can use the ISP scheme to generate a Langevin path ω_L. We discuss briefly how the difference between ω_o and ω_L depends on the combined parameter ξΔt, which can be interpreted as the number of collisions per time step.

In the limit of high friction $ξ m \dot{x} ≫ m \ddot{x}$ ⁠, the Langevin dynamics [Eq. (19)] approaches the overdamped Langevin dynamics [Eq. (11)]. More specifically, in Eq. (19), set $m \ddot{x} = 0$ ⁠, and rearranging yields Eq. (11). However, even though the equation of motion for Langevin dynamics converges to the equation of motion for overdamped Langevin dynamics, the ISP scheme [Eqs. (20) and (21)] does not converge to the Euler–Maruyama scheme [Eq. (12)] in the limit of high friction. By “high friction,” we denote the range of collision rates ξ for which e^−ξΔt ≈ 0 in Eq. (20), but $|) \frac{\nabla V}{ξ m} |) > 0$ ⁠. (As reference, e^−0.1 = 0.904, e⁻¹ = 0.368, and e⁻⁵ = 0.007.) If e^−ξΔt ≈ 0, then also e^−2ξΔt ≈ 0, and Eq. (20) becomes

x_{k + 1} \approx x_{k} - \frac{\nabla V (x_{k})}{ξ m} Δ t + \sqrt{\frac{k_{B} T}{m}} η_{L, k} Δ t .

(28)

The first two terms on the right-hand side are identical to the Euler–Maruyama scheme [Eq. (12)], but the random number term differs from the Euler–Maruyama scheme. Thus, even in the limit of high friction, the two algorithms yield different paths for a given random number sequence η. The difference between a Langevin path ω_L and an overdamped Langevin path ω_o can be scaled by the combined parameter ξΔt. For some value ξΔt > 1, the difference between the two paths becomes minimal before increasing again, but for no value of ξΔt, the two paths fully coincide.

When Langevin integration schemes are used as a thermostat in MD simulations, the optimal friction coefficient should reproduce the expected temperature fluctuations and therefore depends on the system and the simulation box.⁵⁸ Reported collision rates^49,50,59 (while keeping the time step at Δt = 0.002 ps) range from 0.1 ps⁻¹ to ∼100 ps⁻¹, corresponding to ξΔt = 0.0002 to ξΔt = 0.2. However, even for a large collision rate of 100 ps⁻¹, e^−ξΔt = e^−0.2 = 0.819 ≉ 0. For these two reasons—MD simulations are not conducted in the high-friction regime, and even in the high-friction regime, ω_o differs from ω_L—a simulation with the ISP scheme yields a materially different path ensemble than a simulation with the Euler–Maruyama scheme.

C. From a path ω to random numbers η_o and η_L

In Sec. V B, we showed that given a random number sequence η, the path generated by the Euler–Maruyama integration scheme for overdamped Langevin dynamics differs from the path generated by the ISP integration scheme for Langevin dynamics. More relevant for path reweighting is the reverse situation: Given a sample path ω = (x₀, …, x_n) in position space and the parameters of the dynamics (m, V, T, ξ, k_B, and Δt), how does the random number sequence η_o needed to generate ω with the Euler–Maruyama scheme [Eq. (12)] differ from the random number sequence η_L needed to generate the same ω with the ISP scheme [Eqs. (20) and (21)]? An equivalent question is as follows: How does the path probability that ω has been generated by the Euler–Maruyama scheme differ from the path probability that ω has been generated by the ISP scheme? How does this difference affect the path probability ratios between the simulation and a target potential? Figure 2 gives an overview of the quantities we will compare. Note that we dropped the index o or L from the path ω because ω is a given dataset, which will be analyzed using various approaches to calculate the path probabilities.

FIG. 2.

View large Download slide

Overview of path probabilities and path probability ratios for a sample path ω = (x₀, … x_n).

First, we need to discuss whether such a comparison between the ISP scheme and the Euler–Maruyama scheme is even possible. From an algorithmic view point, this is clearly possible because both integrators [Eqs. (12) and (20)] use a single random number per integration time step. The path probabilities are then equal to the probability of the different random number sequences η_L and η_o needed to generate ω. From a physical view point, the answer is not as clear because overdamped Langevin dynamics evolves in position space (x_k), whereas Langevin dynamics evolves in phase space (x_k, $v$ _k). The velocity $v$ _k enters the integration scheme [Eq. (20)] and the path probability [Eq. (22)]. However, $v$ _k is fully determined by the current position x_k and the previous position x_k−1 [Eq. (21)]. Thus, if the initial velocity $v_{0}$ is known, the position trajectory is enough to evaluate the path probability [Eq. (22)], and the comparison to overdamped Langevin dynamics is possible.

We consider the test system described in Sec. V A at the simulation potential V(x) (double-well potential) simulated by the ISP scheme for Langevin dynamics. With ξ = 50 s⁻¹ and Δt = 0.01 s, we have e^−ξΔt = e^−0.5 = 0.607 ≉ 0, meaning that the system is not in the high-friction limit. Figure 3(a) additionally shows that with these parameters $O (ξ m \dot{x}) \approx O (m \ddot{x})$ and also according to the criterion for the stochastic differential equation, the system is not in the high-friction limit.

FIG. 3.

View large Download slide

(a) The acceleration term $m \ddot{x}$ and the friction $ξ m \dot{x}$ for the test system at V(x). (b) Example path ω of length n = 10. (c) Random number sequences η_L (green solid), η_o (blue solid), ${\tilde{η}}_{L}$ (green dashed), and ${\tilde{η}}_{o}$ (blue dashed) that correspond to ω. (d) Path probabilities P_L(ω; Δt|(x₀, $v_{0}$ ⁠)) (green solid), P(ω; Δt|x₀) (blue solid), ${\tilde{P}}_{L} (ω; Δ t | (x_{0}, v_{0}) ()$ (green dashed), and ${\tilde{P}}_{o} (ω; Δ t | x_{0})$ (blue dashed). (e) Path probability ratios: M_L(ω, Δt|(x₀, $v_{0}$ ⁠)) (green) and M_o(ω; Δt|x₀) (blue).

Figure 3(b) shows a sample path ω = (x₀, x₁, …, x₁₀). Figure 3(c) shows the random numbers η_o needed to generate ω with the Euler–Maruyama scheme [blue solid line, calculated using Eq. (15)] and the random numbers η_L needed to generate ω with the ISP scheme [green solid line, calculated using Eq. (24)]. As expected for the low-friction regime, these two random number sequences differ markedly.

Consequently, the path probabilities differ. Figure 3(d) shows the unnormalized path probability for generating ω with the Euler–Maruyama scheme (blue solid line),

\begin{align} P_{o} (ω; Δ t | x_{0}) \\ \sim \exp (- \frac{ξ m}{4 k_{B} T Δ t} \sum_{k = 0}^{n - 1} {(x_{k + 1} - x_{k} + \frac{Δ t}{ξ m} \nabla V (x_{k}))}^{2}), \end{align}

(29)

and for generating ω with the ISP scheme (green solid line),

P_{L} (ω; Δ t | (x_{0}, v_{0})) \sim \exp (- \sum_{k = 0}^{n - 1} \frac{m {(x_{k + 1} - x_{k} - \exp (- ξ Δ t) v_{k} Δ t + (1 - \exp (- ξ Δ t)) \frac{\nabla V (x_{k})}{ξ m} Δ t)}^{2}}{2 k_{B} T (1 - \exp (- 2 ξ Δ t)) Δ t^{2}}),

(30)

where we omitted those factors from Eqs. (13) and (22) that cancel in the path probability ratio. We checked that the path probabilities are consistent with P(η_o) and P(η_L). The two path probabilities diverge from the first simulation step on. After ten integration time steps, they differ by two orders of magnitude. Clearly, P_L(ω; Δt|(x₀, $v_{0}$ ⁠)) cannot be used as an approximation for P_o(ω; Δt|x₀).

However, an interesting observation arises when we consider reweighting ω to the target potential $\tilde{V} (x)$ (triple-well potential). Figure 3(c) shows the random numbers ${\tilde{η}}_{o}$ needed to generate ω with the Euler–Maruyama scheme at $\tilde{V} (x)$ [blue dashed line, calculated using Eq. (16)] and the random numbers ${\tilde{η}}_{L}$ needed to generate ω with the ISP scheme at $\tilde{V} (x)$ [green dashed line, calculated using Eq. (25)]. The corresponding unnormalized path probabilities $\sim {\tilde{P}}_{o} (ω; Δ t | x_{0})$ and $\sim {\tilde{P}}_{L} (ω; Δ t | (x_{0}, v_{0}))$ are shown as dashed lines in Fig. 3(d). Strikingly, a change in the integration scheme from the Euler–Maruyama scheme to ISP has a much stronger influence on the random numbers and the path probability than the modification of the potential energy function. Figure 3(e) shows the path probability ratios, i.e., the ratio between the dashed and the solid lines in Fig. 3(d), for the Euler–Maruyama scheme M_o = M_o(ω; Δt|x₀) = M_o(ω, η_o; Δt|x₀) (blue line) and the ISP scheme M_L = M_L(ω; Δt|(x₀, $v_{0}$ ⁠)) = M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)) (green line). Because, within an integration scheme, the path probability does not change drastically when going from the simulation potential V(x) to the target potential $\tilde{V} (x)$ ⁠, both path probability ratios remain at ≈1 throughout the path and follow similar curves, that is, the path probability ratios for Langevin and overdamped Langevin dynamics are much more similar than the underlying path probabilities.

D. Path reweighting

We return to the scenario described in the Introduction and ask the following: are the two path probability ratios similar enough that we can use M_o as an approximation to M_L in Eq. (1)? Figure 4(a) compares different ways to calculate the path probability ${\tilde{P}}_{L} (ω; Δ t | (x_{0}, v_{0}))$ ⁠, i.e., the probability with which an example path ω would have been generated at the target potential $\tilde{V} (x)$ ⁠. The black line is the reference solution calculated by inserting $\tilde{V} (x)$ into Eq. (22). It is identical to the green dashed line in Fig. 3(d). The green line in Fig. 4(a) shows the reweighted path probability, where we used the exact path probability ratio for the ISP scheme, M_L(ω; Δt|(x₀, $v_{0}$ ⁠)) [Eq. (23)], in Eq. (1). As expected, this reweighted path probability coincides with the directly calculated path probability. The blue line shows the reweighted path probability, where we used the path probability ratio for the Euler–Maruyama scheme, M_o(ω; Δt|x₀) [Eq. (14)], as an approximation to M_L in Eq. (1). The path probability deviates from the reference solution, but overall follows a similar curve.

FIG. 4.

View large Download slide

(a) Reference and reweighted path probabilities for ω for Langevin dynamics. (b) Reference and reweighted first three dominant MSM left eigenfunctions l₁, l₂, and l₃ associated with $\tilde{V} (x)$ for Langevin dynamics. (c) Reference and reweighted implied timescales corresponding to l₂ and l₃.

Figure 4(a) merely serves to illustrate the concepts. With only ten steps, the example path ω is far too short to judge the accuracy of the two path probability ratios for reweighting dynamic properties. We, therefore, constructed MSMs for the target potential $\tilde{V} (x)$ ⁠. The reference solution has been generated from simulations at the target potential $\tilde{V} (x)$ using the ISP scheme. The dominant MSM eigenfunctions and associated implied timescales are shown as black lines in Figs. 4(b) and 4(c). Next, we ran simulations at the simulation potential V(x) using the ISP scheme and constructed a reweighted MSM using the exact reweighting factor M_L(ω; Δt|(x₀, $v_{0}$ ⁠)) [Eq. (23)]. The dominant MSM eigenfunctions are shown as green lines in Fig. 4(b). They exactly match the reference solution. The reweighted implied timescales are shown as green lines in Fig. 4(c) and are in good agreement with the reference solution. Finally, we used the simulation at V(x) to construct a reweighted MSM using the reweighting factor for the Euler–Maruyama scheme M_o(ω; Δt|x₀) [Eq. (14)]. The dominant MSM eigenfunctions are shown as blue lines in Fig. 4(b). The eigenfunctions differ considerably from the reference solution. Most notably, the stationary distribution is not reproduced correctly [blue line in the upper panel of Fig. 4(b)]. The left peak is reduced to a shoulder of the central peak, and the relative heights of central peak and the right peak do not match those of the reference solution. Likewise, the implied timescales [blue line in Fig. 4(c)] are severely underestimated. This indicates that using the path probability ratio for overdamped Langevin dynamics, M_o(ω; Δt|x₀), to reweight Langevin trajectories does not yield acceptable results.

VI. APPROXIMATE PATH PROBABILITY RATIO

A. Derivation and numerical results

With the results from Sec. IV, the exact random number probability ratio M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)) [Eq. (7)] for the ISP scheme is straightforward to evaluate from a simulation at V(x): the random number sequence η = η_L can be recorded during the simulation, and the random number difference Δη = Δη_L is given by Eq. (26). Inserting η_L and Δη_L into Eq. (7) yields M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)). However, Δη_L,k in Eq. (26) is specific to the ISP scheme. If one uses a different Langevin integration scheme to simulate the dynamics at V(x), one needs to adapt Eq. (26) via the strategy outlined in Sec. IV.

Fortunately, the random number difference for overdamped Langevin dynamics Δη_o,k [Eq. (17)] is approximately equal to Δη_L,k for any given perturbation U(x). Figure 3(c) already suggests that. In Appendix D, we show that the difference between $Δ η_{L, k}^{2}$ and $Δ η_{o, k}^{2}$ is, in fact, only of $O (ξ^{4} Δ t^{4})$ so that for ξΔt < 1, we can assume with high accuracy that

\begin{matrix} \begin{gathered} Δ η_{L, k} \approx Δ η_{o, k}, \\ \sqrt{\frac{1}{k_{B} T ξ^{2} m}} \frac{1 - \exp (- ξ Δ t)}{\sqrt{1 - \exp (- 2 ξ Δ t)}} \cdot \nabla U (x_{k}) \approx \sqrt{\frac{Δ t}{2 k_{B} T ξ m}} \cdot \nabla U (x_{k}) . \end{gathered} \end{matrix}

(31)

The difference between Δη_L,k and Δη_o,k is determined by the prefactors in front of ∇U(x_k) in Eq. (31), which are shown as a function of ξΔt in Fig. 5(b). For ξΔt < 1, the two curves are virtually identical.

FIG. 5.

View large Download slide

(a) Sketch of a step x_k → x_k+1 and the quantities of influence for Langevin and overdamped Langevin dynamics. (b) Prefactors of Δη_L,k and Δη_o,k as a function of ξΔt. (c) Absolute difference (absolute error) between the random numbers ⟨|η_o,k − η_L,k|⟩ and the random number differences ⟨|Δη_o,k − Δη_L,k|⟩ as a function of ξΔt. (d) Relative difference (relative error) between the random numbers ⟨|(η_o,k − η_L,k)/η_L,k|⟩ and the random number differences ⟨|(Δη_o,k − Δη_L,k)/Δη_L,k|⟩ as a function of ξΔt.

With the approximation in Eq. (31), we can derive an approximate random number probability ratio, by using the recorded η_L, but substituting Δη_L,k [Eq. (26)] by Δη_o,k [Eq. (17)] in Eq. (7), we obtain

\begin{matrix} M_{L} (ω, η_{L}; Δ t | (x_{0}, v_{0})) & \approx & M_{a p p r o x} (ω, η_{L}; Δ t | x_{0}) \\ = & \exp (- \sum_{k = 0}^{n - 1} \sqrt{\frac{Δ t}{2 k_{B} T ξ m}} \nabla U (x_{k}) \cdot η_{L, k}) \\ \cdot \exp (- \frac{1}{2} \sum_{k = 0}^{n - 1} \frac{Δ t}{2 k_{B} T ξ m} {(\nabla U (x_{k}))}^{2}) . \end{matrix}

(32)

Equation (32) has the same functional form as the random number probability ratio for the Euler–Maruyama scheme M_o(ω, η_o; Δt|x₀) [Eq. (18)], but it uses η_L, the random numbers generated during the ISP simulation, instead of η_o. Equation (32) is the approximation that we used in Refs. 34 and 35 because we had not yet derived M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)) [Eqs. (23) and (27)].

Figure 4 demonstrates the accuracy of the approximate random number probability ratio M_approx(ω, η_L; Δt|x₀) [Eq. (32)] for our test system. The orange dashed line in Fig. 4(a) shows the reweighted path probability for the short example path, where we used M_approx(ω, η_L; Δt|x₀) [Eq. (32)], in Eq. (1). It exactly matches the reference solution (black line).

Next, we constructed a reweighted MSM for the target potential $\tilde{V} (x)$ based on our simulations at the simulation potential V(x) using M_approx(ω, η_L; Δt|x₀) [Eq. (32)] to reweight the transition counts. The dominant MSM eigenfunctions of the reweighted MSM are shown as orange dashed lines in Fig. 4(b). They exactly match the reference solution. The reweighted implied timescales are shown as orange dashed lines in Fig. 4(c) and seem to match the reference solution even better than the ones calculated using the exact path probability ratio [green line in Fig. 4(c)]. However, the difference between the orange dashed line and the green line is likely within statistical uncertainty. In summary, M_approx(ω, η_L; Δt|x₀) is a highly accurate approximation to M_L(ω, η_L; Δt|x₀) for ξΔt < 1. Using M_approx(ω, η_L; Δt|x₀) instead of M_L(ω, η_L; Δt|x₀) could even have the following advantages: (i) the implementation is less error-prone because the functional form of M_approx is simpler than that of M_L and (ii) M_approx might be numerically more stable because the calculation of exponential function on the left-hand side of Eq. (31) is avoided.

B. Intuition

We discuss why M_approx(ω, η_L; Δt|x₀) is a better approximation to M_L(ω, η_L; Δt|x₀) than M_o(ω; Δt|x₀) = M_o(ω, η_o; Δt|x₀). Figure 5(a) shows one integration time step of a stochastic integration scheme from x_k to x_k+1 (black line). From k to k + 1, the system has progressed by Δx = x_k+1 − x_k. In the ISP scheme, this progress is composed of a progress

Δ x_{d r i f t, L} = \exp (- ξ Δ t) v_{k} Δ t - [1 - \exp (- ξ Δ t)] \frac{\nabla V (x_{k})}{ξ m} Δ t

(33)

due to the drift force and the velocity of the system [second and third terms on the right-hand side of Eq. (20)] and a progress

Δ x_{r a n d o m, L} = \sqrt{\frac{k_{B} T}{m} [1 - \exp (- 2 ξ Δ t)]} η_{L, k} Δ t

(34)

due to the random force [fourth term on the right-hand side of Eq. (20)] such that Δx = Δx_drift,L + Δx_random,L. Δx_drift,L and Δx_random,L are illustrated as green solid lines in Fig. 5(a). The probability of generating the step x_k → x_k+1 is determined by Δx_random,L, which is proportional to the random number η_L,k (green solid arrow).

With a different potential energy function $\tilde{V} (x)$ at x_k, the displacement due to the drift force differs from the original Δx_drift,L. To achieve the same overall displacement Δx, Δx_random,L needs to be adjusted (green dotted line). The corresponding random number ${\tilde{η}}_{L, k}$ is shown as a green dotted arrow, and the difference between the two random numbers Δη_L,k is shown as a red line. In path reweighting, one constructs ${\tilde{η}}_{L, k}$ by adding Δη_L,k to η_L,k,

{\tilde{η}}_{L, k} = η_{L, k} + Δ η_{L, k}

(35)

[analogous to Eq. (6)], which then yields the general form of the random number probability ratio in Eq. (7).

An analogous analysis applies to the Euler-Maruyama scheme, where the progress due to the drift force is

Δ x_{d r i f t, o} = - \frac{\nabla V (x_{k})}{ξ m} Δ t

(36)

[second term on the right-hand side of Eq. (12)], and the progress due to the random force is

Δ x_{r a n d o m, o} = \sqrt{\frac{2 k_{B} T}{ξ m}} \sqrt{Δ t} η_{o, k}

(37)

[third term on the right-hand side of Eq. (12)]. In Fig. 5(a), Δx_drift,o and Δx_random,o are illustrated as blue solid lines, and the random number is represented as a blue solid arrow. With a different potential energy function $\tilde{V} (x)$ at x_k, the progress due to the drift force differs from the original Δx_drift,o. To achieve the same overall progress Δx, Δx_random,o needs to be adjusted (blue dotted line). The corresponding random number ${\tilde{η}}_{o, k}$ is shown as a blue dotted arrow, and the difference between the two random numbers Δη_o,k is shown as an orange line.

In Sec. VI A, we have shown that Δη_L,k ≈ Δη_o,k (for ξΔt < 1). Thus, approximating Δη_L,k by Δη_o,k in Eq. (35), or, visually, approximating the red line by the orange line in Fig. 5(a), is valid. However, the displacement due to the drift Δx_drift,o in the Euler–Maruyama scheme can differ strongly from Δx_drift,L in the ISP scheme, and consequently, the random numbers needed to generate the same overall progress Δx differ

η_{L, k} ≉ η_{o, k}

(38)

[blue solid and green solid arrow in Fig. 5(a)]. Consequently, approximating η_L,k by η_o,k in Eq. (35), or visually approximating the green solid arrow by the blue solid arrow in Fig. 5(a), is not valid.

The exact random number probability ratio M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)) [Eq. (27)] uses the exact η_L recorded during the simulation and the exact Δη_L [Eq. (26)]. It therefore yields results that exactly match the reference solutions (green lines in Fig. 4). M_approx(ω, η_L; Δt|x₀) uses the exact η_L recorded during the simulation but approximates Δη_L,k by Δη_o,k. This introduces only a small error but still yields excellent reweighting results in our test system (orange dashed lines in Fig. 4). However, in M_o(ω; Δt|x₀) = M_o(ω, η_o; Δt|x₀), one additionally approximates η_L by η_o. The difference between η_L and η_o is much larger than the difference between Δη_L and Δη_o, and this additional approximation leads to the distorted reweighting results we observed as the blue lines in Fig. 4.

The proportions in Fig. 5(a) are not exaggerated. The black line in Fig. 5(c) shows the average absolute difference between the random numbers ⟨|η_o,k − η_L,k|⟩ as a function of ξΔt. Visually, this is the difference between the green solid arrow and the blue solid arrow in Fig. 5(a). The orange line in Fig. 5(c) shows the average absolute difference between the random number differences ⟨|Δη_o,k − Δη_L,k|⟩, i.e., the difference between the orange and the red line in Fig. 5(a). The graph has been calculated by averaging over a path with 10⁶ time steps. The standard deviations are shown as vertical bars. ⟨|Δη_o,k − Δη_L,k|⟩ is close to 0 for all values of ξΔt, whereas there is a substantial difference between η_L and η_o. ⟨|η_o,k − η_L,k|⟩ has a minimum at ξΔt ≈ 2 because the difference between the Euler–Maruyama scheme and the ISP scheme is minimal for ξΔt ≈ 2 (see Sec. V B). Figure 5(d) shows the corresponding average relative errors. For ξΔt > 1, ⟨|(η_o − η_L)/η_L|⟩ (black line) decreases in accordance with the decrease in the absolute difference ⟨|(η_o − η_L)|⟩ and ⟨|(Δη_o − Δη_L)/Δη_L|⟩ (orange line) increases, reflecting the fact that the approximation [Eq. (31)] does not hold for ξΔt > 1. However, for ξΔt < 1, the region in which MD simulations are conducted, the relative error for the random numbers is much larger than the relative error for the random number difference. This reinforces that the random numbers η_L,k should not be approximated in the path probability ratio but, instead, should be recorded from the simulation at V(x). By contrast, the random number difference Δη_L,k can reliably be approximated by Eq. (31).

VII. MOLECULAR EXAMPLE: BUTANE

The slowest degree of freedom in butane is the torsion around the C₂–C₃ bond, which exhibits three metastable states: the trans-conformation at ϕ = π and the two gauche-conformations at $ϕ = \pm \frac{1}{3} π$ ⁠. Consequently, butane has three dominant MSM eigenvectors, where l₁ corresponds to the stationary density and l₂ and l₃ represent slow transitions along ϕ [Fig. 6(a)]. Because the two gauche-conformations are equally populated, l₂ and l₃ are degenerate [Fig. 6(b)]. We simulated butane in implicit water at three different temperatures, T = 300 K, T = 200 K, and T = 150 K, using direct and biased simulations. As we lower the temperature, we expect that the relative population of the trans-conformation increases, but that otherwise, the overall shape and sign-structure of the eigenvectors remain unchanged.

FIG. 6.

View large Download slide

Dynamics of the torsion angle in butane at T = 300 K, 200 K, and 150 K. (a) Dominant left eigenfunctions l₁, l₂, and l₃ of the MSM along the torsion angle ϕ, obtained by evaluating direct simulations at the target potential, as well as by reweighting biased simulations. (b) Implied timescales corresponding to l₂ and l₃ in (a). Solid lines: mean and shaded area: standard deviation. Standard deviations for the eigenvectors are too small to be shown.

At T = 300 K and T = 200 K, the reweighting results using M_approx(ω, η_L; Δt|x₀) [Eq. (32), orange dashed line] or M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)) [Eq. (23), green solid line] match the MSM obtained by direct simulation. In particular, the eigenvectors are reproduced with very high precision. By contrast, the reweighted results using M_o(ω; Δt|x₀) [Eq. (14), blue line] deviate considerably from the reference MSMs obtained by direct simulations. The stationary distribution l₁ is not reproduced correctly, which then leads to further errors in the dominant eigenvectors l₂ and l₃. The associated implied timescales are underestimated. Moreover, for T = 200 K and T = 150 K, the use of M_o(ω; Δt|x₀) yielded numerically instable transition matrices for lag times of τ > 100 ps. This demonstrates that path reweighting with an appropriate path probability ratio, such as M_approx(ω, η_L; Δt|x₀) or M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)), yields accurate results. However, M_o(ω; Δt|x₀) should not be used as an approximation for the exact path probability ratio M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)).

Note that reweighting results using the approximate probability ratio M_approx(ω, η_L; Δt|x₀) are virtually indistinguishable from the results using the exact probability ratio M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)) for all three temperatures. This confirms our analysis that M_approx(ω, η_L; Δt|x₀) can be used as highly accurate approximation to M_L(ω, η_L; Δt|(x₀, $v_{0}$ ⁠)).

The variation of the temperature from 300 K to 200 K and 150 K illustrates under which circumstances path reweighting is an efficient method. At T = 300 K, many transitions across the torsion angle barriers are observed in the direct simulation. Path reweighting and direct simulation yield identical results. However, path reweighting has a larger statistical uncertainty. At T = 200 K, fewer transitions are observed in the direct simulations, which results in an increased statistical uncertainty in the direct MSM. Finally, at T = 150 K, the transitions in the direct simulation are insufficient to correctly sample the stationary density. The MSM of the direct simulation predicts a higher population for the gauche-conformation at $ϕ = + \frac{1}{3} π$ than for the gauche-conformation at $ϕ = - \frac{1}{3} π$ ⁠, which is clearly a sampling error. This error in the stationary density then leads to vastly incorrect estimates for l₂ and l₃. Additionally, the direct MSM predicts that the degeneracy is lifted. By contrast, the results of the reweighted MSM are in line with what we expect: the gauche-conformations are equally populated, the overall shapes of the dominant eigenvectors correspond to those of the eigenvectors at higher temperatures, and l₂ and l₃ are degenerate. In conclusion, path reweighting in combination with enhanced sampling techniques is a promising tool in situations, where the stationary density cannot be sampled accurately by direct simulation.

VIII. METHODS

A. Simulations of the test system

The test system is a one-dimensional one particle system with mass m = 1 kg and k_BT = 2.494 J (corresponding to k_B = 0.008 314 J/K and T = 300 K). The simulation potential (orange line in Fig. 1) is

V (x) = {(x^{2} - 1)}^{2},

(39)

and the target potential (black line in Fig. 1) is

\tilde{V} (x) = 4 {(x^{3} - \frac{3}{2} x)}^{2} - x^{3} + x .

(40)

For the results in Figs. 3–5, we simulated the system using the ISP scheme [Eqs. (20) and (21)] with a time step of Δt = 0.01 s. The initial conditions were x₀ = 1.50 m, $v_{0}$ = 0 m/s. The number of time steps N_t, the collision rate ξ, and the potential energy function used are summarized in Table II.

TABLE II.

Simulation parameters.

Figures	N_t	ξ	Potential
3(a)	10⁵	50 s⁻¹	V(x)
4(b) and 4(c)	10⁷	50 s⁻¹	V(x)
4(b) and 4(c)	10⁷	50 s⁻¹	$\tilde{V} (x)$
5(c) and 5(d)	10⁷	0.1 s⁻¹–1000 s⁻¹	V(x)

Figures	N_t	ξ	Potential
3(a)	10⁵	50 s⁻¹	V(x)
4(b) and 4(c)	10⁷	50 s⁻¹	V(x)
4(b) and 4(c)	10⁷	50 s⁻¹	$\tilde{V} (x)$
5(c) and 5(d)	10⁷	0.1 s⁻¹–1000 s⁻¹	V(x)

In Fig. 3(a), we computed the acceleration $\ddot{x} = a$ as $a_{k + 1} = \frac{v_{k + 1} - v_{k}}{Δ t}$ ⁠. Figure 3(b) displays the first ten steps of the simulation as an example path ω, and all quantities displayed in Figs. 3(c)–3(e) are calculated from this short path.

The absolute and relative differences of the random numbers in Fig. 5 were calculated as

⟨ | η_{o, k} - η_{L, k} | ⟩ = \frac{1}{N_{t} - 1} \sum_{k = 0}^{N_{t} - 1} | η_{o, k} - η_{L, k} |

(41)

and

⟨|\frac{η_{o, k} - η_{L, k}}{η_{L, k}}|⟩ = \frac{1}{N_{t} - 1} \sum_{k = 0}^{N_{t} - 1} |\frac{η_{o, k} - η_{L, k}}{η_{L, k}}| .

(42)

Analogous equations were used for ⟨|Δη_o,k − Δη_L,k|⟩ and ⟨|(Δη_o,k − Δη_L,k)/Δη_L,k|⟩. η_L,k was recorded during the simulation. We used Eq. (26) to calculate Δη_L,k, Eq. (15) to calculate η_o,k, and Eq. (17) to calculate Δη_o,k.

The reference MSM in Figs. 4(b) and 4(c) has been constructed from the simulation at the target potential $\tilde{V} (x)$ ⁠. The state space has been discretized using a regular grid of 100 microstates (S₁, …, S₁₀₀) in the range −1.7 ≤ x ≤ 1.6. Transition counts between microstates were calculated as

c_{i j} (τ) = \frac{1}{N_{t} - τ} \sum_{k = 0}^{N_{t} - τ} χ_{i} (x_{k}) χ_{j} (x_{k + τ}),

(43)

with

χ_{i} (x) = \{\begin{matrix} 1 if x \in S_{i} \\ 0 else, \end{matrix})

(44)

where x_k is the trajectory and lag time τ = 200 steps. The resulting count matrix C(τ) was symmetrized as C(τ) + C^⊤(τ) to enforce detailed balance and row-normalized to obtain the MSM transition matrix T(τ). The dominant MSM eigenvectors l_i and associated eigenvalues λ_i(τ) were calculated from T(τ) using a standard eigenvalue solver, and the implied timescales were calculated as t_i = −τ/ln(λ_i(τ)).

The reweighted MSMs in Figs. 4(b) and 4(c) have been constructed from the simulation at the simulation potential V(x) using the same grid and lag time as for the reference MSM. Transition counts between microstates were counted and reweighted as^34,35

\begin{align} {\tilde{c}}_{i j} (τ) = \frac{1}{N_{t} - τ} \sum_{k = 0}^{N_{t} - τ} W ((x_{k}, x_{k + 1}, \dots, x_{k + τ}); Δ t | (x_{k}, v_{k})) \\ χ_{i} (x_{k}) χ_{j} (x_{k + τ}) . \end{align}

(45)

The weight W is defined as

\begin{align} W ((x_{k}, x_{k + 1}, \dots, x_{k + τ}); Δ t | (x_{k}, v_{k})) \\ = g (x_{k}) \cdot M ((x_{k}, x_{k + 1}, \dots, x_{k + τ}); Δ t | (x_{k}, v_{k})), \end{align}

(46)

with M being the path probability ratio [Eq. (3)] and g being

g (x_{k}) = \exp (- \frac{U (x_{k})}{k_{B} T}),

(47)

where the perturbation U is defined in Eq. (2). The remaining procedure was analogous to the reference MSM.

B. Butane—Direct simulations

We performed all-atom MD simulations of n-butane in implicit water using the OpenMM 7.4.1⁵⁴ simulation package. The GAFF (Generalized Amber Force Field)⁶⁰ was used to model butane, and the GBSA (Generalized Born Surface Area) model⁶¹ was used to model implicit water. Interactions beyond 1 nm were truncated. The trajectory was propagated according to the ISP integration scheme for a 3N-dimensional system,

\begin{matrix} x_{k + 1}^{i} & = & x_{k}^{i} + \exp (- ξ Δ t) v_{k}^{i} Δ t - [1 - \exp (- ξ Δ t)] \frac{\nabla_{i} V (x_{k})}{ξ m_{i}} Δ t \\ + \sqrt{\frac{k_{B} T}{m_{i}} [1 - \exp (- 2 ξ Δ t)]} η_{L, k}^{i} Δ t, \end{matrix}

(48)

v_{k + 1}^{i} = \frac{x_{k + 1}^{i} - x_{k}^{i}}{Δ t},

(49)

with i = 1, 2, …, 3N and N being the number of atoms. $x_{k}^{i}$ ⁠, $v_{k}^{i}$ ⁠, and $η_{k}^{i}$ are the position, velocity, and random number along dimension i at iteration step k, m_i is the mass of dimension i, and ∇_iV(x_k) denotes the gradient of V(x_k) along dimension i measured at the position x_k, with $x \in R^{3 N}$ ⁠. We implemented the ISP integration scheme using the simtk.openmm.openmm.CustomIntegrator⁶² class of OpenMM. The collision rate was ξ = 10 ps⁻¹. The simulation time step was Δt = 0.002 ps. Positions were written to disk every txout = 50 steps = 0.1 ps. We generated three trajectories with 500 ns each at T = 300 K, T = 200 K, and T = 150 K. These direct simulations correspond to simulations at the target potential $\tilde{V} (x)$ ⁠.

For the analysis, we cut each trajectory into five pieces of length 100 ns. For each 100-ns-trajectory, we constructed a MSM following the procedure outlined in Sec. VIII A. As state space we chose the C₂–C₃ dihedral angle ϕ, which we discretized using a regular grid of 100 microstates in the range 0 ≤ ϕ ≤ 2π. This resulted in five MSMs for each temperature. Figure 6 shows the mean and the standard deviation of the first three left MSM eigenvectors (evaluated at lag time τ = 1 ps) and the mean and the standard deviations of the associated implied timescale.

C. Butane—Path reweighting

We biased the simulations along the C₂–C₃ dihedral angle ϕ. To generate the bias potential U(ϕ), we constructed a histogram of the free-energy function $\tilde{F} (ϕ)$ ⁠,

\tilde{F} (ϕ) = - k_{B} T \ln (\tilde{p} (ϕ)),

(50)

where $\tilde{p} (ϕ)$ is the stationary density along ϕ as measured from the 500 ns direct simulations at T = 300 K. Fitting the histogram with a third order Fourier series yielded

\begin{matrix} {\tilde{F}}_{300 K} (ϕ) & = & 8.985 + 3.122 \cos (ω ϕ) + 0.959 \cos (2 ω ϕ) \\ + 7.742 \cos (3 ω ϕ) + 0.095 \sin (ω ϕ) \\ + 0.047 \sin (2 ω ϕ) + 0.002 \sin (3 ω ϕ), \end{matrix}

(51)

with ω = 0.989. The same procedure for the simulation at T = 200 K yielded

\begin{matrix} {\tilde{F}}_{200 K} (ϕ) & = & 8.311 + 2.847 \cos (ω ϕ) + 0.841 \cos (2 ω ϕ) \\ + 7.697 \cos (3 ω ϕ) + 0.046 \sin (ω ϕ) \\ + 0.026 \sin (2 ω ϕ) + 0.004 \sin (3 ω ϕ), \end{matrix}

(52)

with ω = 0.989. ${\tilde{F}}_{300 K} (ϕ)$ and ${\tilde{F}}_{200 K} (ϕ)$ are almost identical. The simulation at T = 150 K did not yield a converged stationary density, and thus, no free-energy function was constructed for this temperature, and instead, ${\tilde{F}}_{300 K} (ϕ)$ was used.

The biased simulations were carried out with the potential

V_{α} (x) = \tilde{V} (x) - α \cdot \tilde{F} (ϕ (x)),

(53)

where $\tilde{V} (x)$ is the target potential and α ∈ [0, 1] specifies the bias strength. V_α(x) corresponds to the “simulation potential” within the terminology of this paper; thus,

U (ϕ (x)) = α \cdot \tilde{F} (ϕ (x)) .

(54)

α was set to 0.1 in all biased simulations, corresponding to “10% of the full metadynamics potential.” We carried out biased simulations at three temperatures T = 300 K, T = 200 K, and T = 150 K, with bias potentials $U_{300 K} (ϕ) = 0.1 \cdot {\tilde{F}}_{300 K} (ϕ)$ ⁠, $U_{200 K} (ϕ) = 0.1 \cdot {\tilde{F}}_{200 K} (ϕ)$ ⁠, and $U_{150 K} (ϕ) = 0.1 \cdot {\tilde{F}}_{300 K} (ϕ)$ ⁠. All other simulation parameters were as described in Sec. VIII B.

The path probability ratios for the biased simulations were calculated on-the-fly^34,35 and were written to disk at the same frequency txout as the positions. For the approximate path probability ratio M_approx, we calculated

\begin{align} M_{a p p r o x} (b) = & \sum_{i = 1}^{3 N} \sum_{k = (b - 1) \cdot t x o u t}^{b \cdot t x o u t - 1} (- \sqrt{\frac{Δ t}{2 k_{B} T ξ m_{i}}} \nabla_{i} U (x_{k}) η_{L, k}^{i}) \\ (- \frac{Δ t}{4 k_{B} T ξ m_{i}} {(\nabla_{i} U (x_{k}))}^{2}) \end{align}

(55)

and constructed the complete path probability ratio as

M_{a p p r o x} (ω, η_{L}; Δ t | x_{0}) = \exp (\sum_{b = 1}^{A} M_{a p p r o x} (b))

(56)

during the construction of the MSM, where $A \in N$ such that τ = A · txout·Δt.

For the Langevin path probability ratio M_L, we calculated the terms

M_{L, 1} (b) = \sum_{i = 1}^{3 N} \sum_{k = (b - 1) \cdot t x o u t}^{b \cdot t x o u t - 1} (x_{k + 1}^{i} - x_{k}^{i}) (\nabla_{i} \tilde{V} (x_{k}) - \nabla_{i} V (x_{k})),

(57)

M_{L, 2} (b) = \sum_{i = 1}^{3 N} \sum_{k = (b - 1) \cdot t x o u t}^{b \cdot t x o u t - 1} v_{k}^{i} (\nabla_{i} \tilde{V} (x_{k}) - \nabla_{i} V (x_{k})),

(58)

M_{L, 3} (b) = \sum_{i = 1}^{3 N} \sum_{k = (b - 1) \cdot t x o u t}^{b \cdot t x o u t - 1} \frac{({(\nabla_{i} \tilde{V} (x_{k}))}^{2} - {(\nabla_{i} V (x_{k}))}^{2})}{m_{i}}

(59)

and constructed the complete path probability ratio as

\begin{matrix} M_{L} (ω, Δ t | x_{0}) \\ = \exp [\sum_{b = 1}^{A} (- \frac{M_{L, 1} (b)}{k_{B} T ξ (1 + \exp (- ξ Δ t)) Δ t})) \\ ((+ \frac{M_{L, 2} (b)}{k_{B} T ξ (1 + \exp (ξ Δ t))} - \frac{\exp (ξ Δ t) - 1}{\exp (ξ Δ t) + 1} \frac{M_{L, 3} (b)}{2 k_{B} T ξ^{2}})] \end{matrix}

(60)

during the construction of the MSM, where $A \in N$ such that τ = A · txout·Δt.

For the overdamped Langevin path probability ratio M_o, we calculated the terms

M_{o, 1} (b) = \sum_{i = 1}^{3 N} \sum_{k = (b - 1) \cdot t x o u t}^{b \cdot t x o u t - 1} (x_{k + 1}^{i} - x_{k}^{i}) (\nabla_{i} \tilde{V} (x_{k}) - \nabla_{i} V (x_{k})),

(61)

M_{o, 2} (b) = \sum_{i = 1}^{3 N} \sum_{k = (b - 1) \cdot t x o u t}^{b \cdot t x o u t - 1} \frac{({(\nabla_{i} \tilde{V} (x_{k}))}^{2} - {(\nabla_{i} V (x_{k}))}^{2})}{m_{i}}

(62)

and constructed the complete path probability ratio as

M_{o} (ω, Δ t | x_{0}) = \exp [\sum_{b = 1}^{A} (- \frac{M_{o, 1} (b)}{2 k_{B} T} - \frac{M_{o, 2} (b) Δ t}{4 k_{B} T ξ})]

(63)

during the construction of the MSM, where $A \in N$ such that τ = A · txout·Δt.

For the analysis, we cut each trajectory into five pieces of length 100 ns. For each 100-ns-trajectory, we constructed a MSM following the procedure outlined in Sec. VIII A. As state space we chose the C₂–C₃ dihedral angle ϕ, which we discretized using a regular grid of 100 microstates in the range 0 ≤ ϕ ≤ 2π. Transition counts between microstates were counted and reweighted as described in Eq. (45) with x_k = ϕ_k and

g (ϕ_{k}) = \exp (- \frac{U (ϕ_{k})}{k_{B} T}) = \exp (- \frac{0.1 \cdot \tilde{F} (ϕ_{k})}{k_{B} T}),

(64)

where ϕ_k is the first entry in the path of length τ. This resulted in five reweighted MSMs for each temperature. Figure 6 shows the mean and the standard deviation of the first three left MSM eigenvectors (evaluated at lag time τ = 1 ps) and the mean and the standard deviations of the associated implied timescale.

Example scripts for simulation and analysis are included as the supplementary material.

IX. CONCLUSION AND OUTLOOK

We have presented two strategies to derive the path probability ratio M_L for the ISP scheme. In the first strategy, the correctly normalized path probability is derived by integrating out the random number η_k from the one-step transition probability. In the second strategy, the equations for the ISP scheme are solved for η_k, and the resulting transformation is used as a change in variables on the Gaussian probability density of the random numbers. This yields an unnormalized path probability. The path probability ratio M_L is then calculated as the ratio between the path probability at the target potential ${\tilde{P}}_{L} (ω_{L}; Δ t | (x_{0}, v_{0}))$ and the path probability at the simulation potential P_L(ω_L; Δt|(x₀, $v_{0}$ ⁠)).

With M_L, we are now able to perform exact path reweighting for trajectories generated by the ISP integration scheme. Moreover, the two strategies serve as a blueprint for deriving path probability ratios for other Langevin integration schemes, which use Gaussian white noise.^{44–47,49–53} Thus, path reweighting can now readily be applied to MD simulation conducted at the NVT ensemble thermostatted with a stochastic thermostat.

We compared the approximate path probability ratio M_approx that we used in earlier publications^34,35 to the exact path probability ratio M_L, both analytically and numerically. We showed that the two expressions only differ by $O (ξ^{4} Δ t^{4})$ ⁠. Thus, M_approx is an excellent approximation to M_L for Langevin MD simulations. To understand why the approximation is so good, we showed that the random number η_k needed to generate a given step x_k → x_k+1 is highly dependent on the integration scheme. However, Δη_k, the difference between the random number ${\tilde{η}}_{k}$ at $\tilde{V} (x)$ and the random number η_k at V(x), has about the same value in the ISP scheme and in the Euler–Maruyama scheme.

In M_approx, one uses the random numbers directly recorded during the simulation at V(x), which does not introduce any error and approximates Δη_k by the expression from the Euler–Maruyama scheme Δη_o,k to construct ${\tilde{η}}_{k}$ ⁠.

We have chosen the ISP algorithm for the present analysis in order to be consistent with our previous work.^34,35 However, the same strategy can be used to derive the path probability ratio for other Langevin integrators.^{44–47,49–53} Specifically, solve the integrator equations for the random number η_k; from there, derive an expression for Δη_k, record η_k during the simulation at V(x) and calculate Δη_k on the fly, and insert η_k and Δη_k into Eq. (7). For a large application of path reweighting, using a modern Langevin integrator is likely worthwhile, such as the the BAOAB method⁵¹ (or alternatively the VRORV method⁵³). This method is exceptionally efficient at sampling the configurational stationary distribution, which allows for increasing the time step.^51,53

It is tempting to speculate that Δη_k for other Langevin integration schemes could also have about the same value as Δη_o,k for the Euler–Maruyama scheme. This would open up a route to a general approximate path probability ratio M_≈ and would eliminate the problem that the path probability needs to be adapted for each integration scheme. On the other hand, the structure of the ISP scheme is closer to that of the Euler–Maruyama scheme than most other Langevin integrators. Whether the approximate path probability can indeed be generalized to these integrators is, therefore, not yet obvious and needs to be checked carefully.

Our one-dimensional test system and our molecular system showed that the accuracy of the reweighting sensitively depends on an accurate representation of η_k in the path probability ratio. For example, reweighting a Langevin path by the path probability ratio for the Euler–Maruyama scheme yielded very distorted results. Neither the MSM eigenvectors nor the implied timescales were reproduced correctly. It is, however, possible that the distortion is less severe in the limit of infinite sampling of the combined space of molecular states and random numbers (probably less relevant to actual applications) or if the dynamics is projected onto a reaction coordinate before the reweighted dynamical properties are evaluated (probably very relevant to actual applications).

We used path reweighting to reweight MSMs. The dynamical property that is reweighted to estimate a transition probability is a correlation function. It is important to point out that correlation functions are a combination of path ensemble averages, where the path is conditioned on a particular initial state (x₀, $v_{0}$ ⁠) and a phase-space ensemble average for the initial states. Thus, the total reweighting factor for MSMs is combined of the path probability ratio M for the path ensemble average and the Boltzmann probability ratio for the phase-space ensemble average g(x) [Eq. (47)].27^,32–34 Even though the reweighting of the path ensemble average can be made exact, by averaging over the initial states within a microstate, one assumes local equilibrium within this microstate.²³ Beyond local equilibrium, the formalism has been extended to reweighting transition probabilities from non-equilibrium steady-state simulations.⁶³

When is the combination of enhanced sampling and path reweighting more efficient than a direct simulation? This depends on the uncertainty of the transition counts estimated from a direct simulation [Eq. (43)] compared to the uncertainty of the reweighted transition counts [Eq. (45)]. The molecular example demonstrated that path reweighting is particularly useful if the stationary density cannot be sampled accurately by direct simulation with the available computer resources. Furthermore, the efficiency of path reweighting increases if the number of transitions at the enhanced sampling simulation is large compared to the direct simulation and if the weights W = g · M are not too small. The path probability ratio M decreases with the path length τ and with the dimensionality of the bias potential U. The path length is kept short by combining path reweighting with MSMs and can be further limited by using advanced MSM discretization techniques.^64–66 The bottleneck for the dimensionality U already occurs at the stage of sampling because most enhanced sampling techniques¹⁰ are limited to very low-dimensional biases in practice. Note that increasing the dimensionality of the overall system does not lower the efficiency of the path reweighting. The question of how strong the bias should be is more difficult to answer. Strong biases increase the transitions in the biased simulation but reduce both g and M. In Ref. 35, we empirically found that a bias of ca. 10% of the full metadynamics biasing potential yielded optimal results, but this will likely depend on the system. Here, we have restricted ourselves to systems with low barriers in the order of k_BT so that we could generate reference solutions by direct simulation. However, we believe that path reweighting is most useful for systems with large barriers that cannot be sampled by direct simulation. An example is the β-hairpin folding equilibrium in Ref. 35.

Path reweighting is closely related to path sampling techniques, in particular path sampling techniques that aim at optimizing the path action.^67–70 The combination of enhanced sampling, path sampling, and path reweighting might change the way we explore the molecular state space and investigate rare events.

SUPPLEMENTARY MATERIAL

See the supplementary material for an example OpenMM script and the corresponding Python3 scripts to construct a reweighted MSM.

DEDICATION

This paper is dedicated to Dr. Irina V. Gopich, a master of stochastic processes. Her work has influenced the way scientists in the field think about the dynamics of molecules—in simulation and in experiment.

ACKNOWLEDGMENTS

The authors would like to thank Luca Donati and Marcus Weber for helpful comments on this manuscript. This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy—EXC 2008—390540038—UniSysCat and through grant CRC 1114 “Scaling Cascades in Complex Systems,” Project Number 235221301, Project B05” Origin of scaling cascades in protein dynamics.”

DATA AVAILABILITY

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

APPENDIX A: LANGEVIN LEAPFROG AND THE ISP SCHEME

Izaguirre, Sweet, and Pande developed the following Langevin Leapfrog algorithm:

\begin{matrix} v_{k + \frac{1}{2}} & = & \exp (- ξ \frac{Δ t}{2}) v_{k} - [1 - \exp (- ξ \frac{Δ t}{2})] \frac{\nabla V (x_{k})}{ξ m} \\ + \sqrt{\frac{k_{B} T}{m} [1 - \exp (- ξ Δ t)]} η_{k}, \end{matrix}

(A1)

x_{k + 1} = x_{k} + v_{k + \frac{1}{2}} Δ t,

(A2)

\begin{matrix} v_{k + 1} & = & \exp (- ξ \frac{Δ t}{2}) v_{k + \frac{1}{2}} - [1 - \exp (- ξ \frac{Δ t}{2})] \frac{\nabla V (x_{k + 1})}{ξ m} \\ + \sqrt{\frac{k_{B} T}{m} [1 - \exp (- ξ Δ t)]} η_{k + 1} \end{matrix}

(A3)

[Eqs. (14)–(16) in Ref. 49]. First, the velocity $v_{k + \frac{1}{2}}$ is updated by a half step using $v$ _k, x_k, and a random number η_k [Eq. (A1)]. Then, the position update to x_k+1 is computed from x_k, assuming a constant velocity $v_{k + \frac{1}{2}}$ in the interval [k, k + 1] [Eq. (A2)]. Finally, the remaining half step of the velocities to $v$ _k+1 is computed using x_k+1, $v_{k + \frac{1}{2}}$ ⁠, and a new random number η_k+1 [Eq. (A3)].

This Langevin Leapfrog algorithm has been converted to the following full-step scheme in the C⁺⁺ CpuLangevinDynamics class of OpenMM:⁷¹

\begin{matrix} v_{k + 1} & = & \exp (- ξ Δ t) v_{k} - [1 - \exp (- ξ Δ t)] \frac{\nabla V (x_{k})}{ξ m} \\ + \sqrt{\frac{k_{B} T}{m} [1 - \exp (- 2 ξ Δ t)]} η_{k}, \end{matrix}

(A4)

x_{k + 1} = x_{k} + v_{k + 1} Δ t,

(A5)

where the velocities are propagated by a full step [i.e., Δt/2 in Eq. (A1) is replaced by Δt, and Δt in Eq. (A1) is replaced by 2Δt] and the position update is based on $v$ _k rather than on $v_{k + \frac{1}{2}}$ ⁠. The second half step for the velocities [Eq. (A3)] is omitted. This integration scheme only uses a single random number per iteration. Equations (A4) and (A5) are the integration scheme we used in Refs. 34 and 35. To distinguish it from the original Langevin Leapfrog scheme [Eqs. (A1)–(A3)], we will refer to Eqs. (A4) and (A5) as the “ISP scheme.”

To be able to analyze the path probability as a function of the positions, we rearrange Eqs. (A4) and (A5) such that we first update the positions using a stochastic step [replace $v$ _k+1 in Eq. (A5) by Eq. (A4)] and then update the velocity as finite difference [rearrange Eq. (A5) with respect to $v$ _k+1]. This yields Eqs. (20) and (21).

APPENDIX B: PATH PROBABILITY FOR LANGEVIN DYNAMICS

We derive the closed-form expression for P_L(ω_L; Δt|(x₀, $v_{0}$ ⁠)) in Eq. (22) from the integration scheme [Eqs. (20) and (21)] by following the approach in Ref. 57. As a first step, we derive a closed-form expression for the one-step probability P_L(x_k+1, $v$ _k+1; Δt|(x_k, $v$ _k)) of observing a step (x_k, $v$ _k) → (x_k+1, $v$ _k+1). According to Eqs. (20) and (21), the tuple (x_k+1, $v$ _k+1) at iteration step k + 1 is entirely determined by the tuple (x_k, $v$ _k) at iteration step k if additionally the random number η_k is known. Thus, P_L(x_k+1, $v$ _k+1; Δt|(x_k, $v$ _k, η_k)), i.e., the one-step probability with fixed random number η_k, is a Dirac delta function centered at (x_k+1, $v$ _k+1). Our strategy is to derive a closed-form expression for this Dirac delta function using Eqs. (20) and (21) and to integrate out the dependency on η_k. In this Appendix, we omit the index L in η_L,k to simplify the notation.

We reformulate the two-dimensional probability P_L(x_k+1, $v$ _k+1; Δt|(x_k, $v$ _k, η_k)) as a product of two one-dimensional probabilities,

\begin{matrix} P_{L} (x_{k + 1}, v_{k + 1}; Δ t | (x_{k}, v_{k}, η_{k})) & = & P_{L} (v_{k + 1}; Δ t | (x_{k + 1}, x_{k}, v_{k}, η_{k})) \\ \cdot P_{L} (x_{k + 1}; Δ t | (x_{k}, v_{k}, η_{k})) \end{matrix}

(B1)

using the rule P(A, B|C) = P(A|B, C) · P(B|C), with A = $v$ _k+1, B = x_k+1, and C = (x_k, $v$ _k, η_k). This rule is the extension of the conditional probability P(A, B) = P(A|B) · P(B) to an additional condition C. The first factor is a Dirac delta function constrained to Eq. (21),

\begin{align} P_{L} (v_{k + 1}; Δ t | (x_{k + 1}, x_{k}, v_{k}, η_{k})) & = P_{L} (v_{k + 1}; Δ t | (x_{k + 1}, x_{k})) \\ = δ (v_{k + 1} - \frac{x_{k + 1} - x_{k}}{Δ t}), \end{align}

(B2)

where the first equality emphasizes that $v$ _k+1 does not depend on η_k or $v$ _k in Eq. (21). Note that the probability of the velocity $v$ _k+1 [Eq. (B2)] does not depend on a random number, which mirrors our previous observation that $v$ _k+1 is not treated as a random variable in Eq. (21). The second factor in Eq. (B1) is a Dirac delta function constrained to Eq. (20),

\begin{matrix} P_{L} (x_{k + 1}; Δ t | (x_{k}, v_{k}, η_{k})) & = & δ (x_{k + 1} - x_{k} - \exp (- ξ Δ t) v_{k} Δ t) \\ (+ [1 - \exp (- ξ Δ t)] \frac{\nabla V (x_{k})}{ξ m} Δ t) \\ (- \sqrt{\frac{k_{B} T}{m} [1 - \exp (- 2 ξ Δ t)]} η_{k} Δ t) . \end{matrix}

(B3)

Reinserting the two factors into Eq. (B1) yields the desired closed-form expression for P_L(x_k+1, $v$ _k+1; Δt|(x_k, $v$ _k; η_k)). Since we know that the random numbers η_k are drawn from a Gaussian distribution P(η_k) with zero mean and unit variance

P (η_{k}) = N^{- 1} \exp (- \frac{η_{k}^{2}}{2}), N = \sqrt{2 π},

(B4)

we can average out the random number dependency in Eq. (B1) to obtain the one-step probability,

\begin{matrix} P_{L} (x_{k + 1}, v_{k + 1}; Δ t | (x_{k}, v_{k})) \\ = \int_{- \infty}^{\infty} d η_{k} P (η_{k}) P_{L} (x_{k + 1}, v_{k + 1}; Δ t | (x_{k}, v_{k}, η_{k})) \\ = δ (v_{k + 1} - \frac{x_{k + 1} - x_{k}}{Δ t}) \\ \cdot \int_{- \infty}^{\infty} d η_{k} P_{η} (η_{k}) P_{L} (x_{k + 1}; Δ t | (x_{k}, v_{k}, η_{k})) . \end{matrix}

(B5)

The challenge lies in solving the integral in this equation. The solution, which is detailed in Appendix C, yields the closed-form expression for the one-step probability,

\begin{matrix} P_{L} (x_{k + 1}, v_{k + 1}; Δ t | (x_{k}, v_{k})) & = & δ (v_{k + 1} - \frac{x_{k + 1} - x_{k}}{Δ t}) \cdot \sqrt{\frac{m}{2 π k_{B} T Δ t^{2} (1 - \exp (- 2 ξ Δ t))}} \\ \times \exp (- \frac{m {(x_{k + 1} - x_{k} - \exp (- ξ Δ t) v_{k} Δ t + (1 - \exp (- ξ Δ t)) \frac{\nabla V (x_{k})}{ξ m} Δ t)}^{2}}{2 k_{B} T (1 - \exp (- 2 ξ Δ t)) Δ t^{2}}) . \end{matrix}

(B6)

Applying the Chapman–Kolmogorov equation⁷² recursively to the one-step probability yields the closed-form expression for the path probability P_L(ω_L; Δt|(x₀, $v_{0}$ ⁠)), shown in Eq. (22).

APPENDIX C: SOLVING THE DOUBLE INTEGRAL

We compute the integral

P_{L} (x_{k + 1}; Δ t | (x_{k}, v_{k})) = \int_{- \infty}^{\infty} d η_{k} P (η_{k}) P_{L} (x_{k + 1}; Δ t | (x_{k}, v_{k}, η_{k}))

(C1)

from Eq. (B5). First, we replace P(η_k) according to Eq. (B4). Second, we substitute P_L(x_k+1; Δt|(x_k, $v$ _k, η_k)), which is a δ-function [Eq. (B3)], with its Fourier transform

δ (z - z^{'}) = \int_{- \infty}^{+ \infty} \frac{d w}{2 π} \exp (i w (z - z^{'})),

(C2)

where z = x_k+1 and z′ is equal to the right-hand side of Eq. (20). This yields a double integral whose outer integral is with respect to $w$ ⁠, while the inner integral is with respect to η_k,

\begin{matrix} P_{L} (x_{k + 1}; Δ t | (x_{k}, v_{k})) \\ = \int_{- \infty}^{+ \infty} \frac{d w}{2 π} \int_{- \infty}^{+ \infty} \frac{d η_{k}}{N} \exp (- \frac{η_{k}^{2}}{2}) \\ \times \exp (i w [x_{k + 1} - x_{k} - \exp (- ξ Δ t) v_{k} Δ t)) \\ ((+ [1 - \exp (- ξ Δ t)] \frac{\nabla V (x_{k})}{ξ m} Δ t)) \\ ((- \sqrt{\frac{k_{B} T}{m} [1 - \exp (- 2 ξ Δ t)]} η_{k} Δ t]) \\ = \int_{- \infty}^{+ \infty} \frac{d w}{2 π} \exp (i w B) \int_{- \infty}^{+ \infty} \frac{d η_{k}}{N} \exp (- \frac{η_{k}^{2}}{2} - i w R η_{k}), \end{matrix}

(C3)

where we moved all terms that do not depend on η_k out of the inner integral and defined the abbreviations,

\begin{matrix} B & = & [x_{k + 1} - x_{k} - \exp (- ξ Δ t) v_{k} Δ t) \\ (+ [1 - \exp (- ξ Δ t)] \frac{\nabla V (x_{k})}{ξ m} Δ t], \\ R & = & Δ t \sqrt{\frac{k_{B} T}{m} [1 - \exp (- 2 ξ Δ t)]} . \end{matrix}

(C4)

Both integrals in Eq. (C3) can be solved with the completing-the-square technique for Gaussian integrals. The goal of this technique is to expand and rearrange the inner integral such that we can use the analytic solution

\int_{- \infty}^{\infty} d x \exp (- a {(x \pm b)}^{2}) = \sqrt{\frac{π}{a}} for a, b \in R .

(C5)

This can be achieved by a systematic step-to-step procedure that can be applied to all Gaussian integrals of this type,

\begin{matrix} \int_{- \infty}^{+ \infty} \frac{d η_{k}}{N} \exp (- \frac{η_{k}^{2}}{2} - i w R η_{k}) \\ = \int_{- \infty}^{+ \infty} \frac{d η_{k}}{N} \exp (- \frac{1}{2} [η_{k}^{2} + 2 i w R η_{k} \overset{= 0}{\overset{︷}{+ i^{2} w^{2} R^{2} - i^{2} w^{2} R^{2}}}]) \\ = \exp (- \frac{w^{2} R^{2}}{2}) \int_{- \infty}^{+ \infty} \frac{d η_{k}}{N} \exp (- \frac{1}{2} {(η_{k} + i w R)}^{2}) \\ = \exp (- \frac{w^{2} R^{2}}{2}) \frac{1}{N} \sqrt{2 π} \\ = \exp (- \frac{w^{2} R^{2}}{2}) . \end{matrix}

(C6)

In the first line, we isolate $η_{k}^{2}$ by factoring out $- \frac{1}{2}$ and complete the first binomial formula by adding a zero. Then, we separate the exponent into the binomial formula and the term $\exp (- \frac{w^{2} R^{2}}{2})$ ⁠, which can be moved in front of the integral because it does not depend on η_k. In the third line, we solve the remaining integral using Eq. (C5), which can be further simplified by inserting the normalization constant of the Gaussian distribution: $N = \sqrt{2 π}$ ⁠.

Inserting Eq. (C6) into Eq. (C3) yields the outer integral

\int_{- \infty}^{+ \infty} \frac{d w}{2 π} \exp (i w B) \exp (- \frac{w^{2} R^{2}}{2}) = \int_{- \infty}^{+ \infty} \frac{d w}{2 π} \exp (- \frac{w^{2} R^{2}}{2} + i w B),

which is solved using the same procedure,

\begin{matrix} \int_{- \infty}^{+ \infty} \frac{d w}{2 π} \exp (- \frac{w^{2} R^{2}}{2} + i w B) \\ = \int_{- \infty}^{+ \infty} \frac{d w}{2 π} \exp (- \frac{R^{2}}{2} [w^{2} + \frac{2 i w B}{R^{2}} \overset{= 0}{\overset{︷}{+ \frac{i^{2} B^{2}}{R^{4}} - \frac{i^{2} B^{2}}{R^{4}}}}]) \\ = \exp (- \frac{B^{2}}{2 R^{2}}) \int_{- \infty}^{\infty} \frac{d w}{2 π} \exp (- \frac{R^{2}}{2} {(w + \frac{i B}{R^{2}})}^{2}) \\ = \exp (- \frac{B^{2}}{2 R^{2}}) \frac{1}{2 π} \sqrt{\frac{2 π}{R^{2}}} \\ = \sqrt{\frac{1}{2 π R^{2}}} \exp (- \frac{B^{2}}{2 R^{2}}) . \end{matrix}

(C7)

Inserting the expressions for the constants R and B [Eq. (C4)] yields

\begin{matrix} P_{L} (x_{k + 1}; Δ t | (x_{k}, v_{k})) \\ = \sqrt{\frac{m}{2 π k_{B} T Δ t^{2} (1 - \exp (- 2 ξ Δ t))}} \\ \times \exp (- \frac{m {(x_{k + 1} - x_{k} - \exp (- ξ Δ t) v_{k} Δ t + (1 - \exp (- ξ Δ t)) \frac{\nabla V (x_{k})}{ξ m} Δ t)}^{2}}{2 k_{B} T (1 - \exp (- 2 ξ Δ t)) Δ t^{2}}) . \end{matrix}

(C8)

This is inserted into Eq. (B5) to yield Eq. (B6).

APPENDIX D: PROOF OF EQ. (31)

\begin{matrix} \frac{{(1 - e^{- x})}^{2}}{x \cdot (1 - e^{- 2 x})} = \frac{1}{2} - \frac{x^{2}}{24} + \frac{x^{4}}{240} \pm O (x^{5}), \\ {(1 - e^{- x})}^{2} = \frac{x}{2} \cdot (1 - e^{- 2 x}) - \frac{x^{2}}{24} \cdot x \cdot (1 - e^{- 2 x}) \\ \pm O (x^{4}) \cdot x \cdot (1 - e^{- 2 x}), \\ {(1 - e^{- x})}^{2} = \frac{x}{2} \cdot (1 - e^{- 2 x}) - O (x^{4}) . \end{matrix}

(D1)

The first line shows the Taylor expansion of the expression on the right-hand side. To obtain the second line, we multiplied by $x \cdot (1 - e^{- 2 x})$ ⁠. In the third line, we used the fact that the leading term of the Taylor expansion of $x \cdot (1 - e^{- 2 x})$ is 2x², thus yielding an error of $O (x^{4})$ ⁠. Substituting x = ξΔt yields

\begin{gathered} {(1 - e^{- ξ Δ t})}^{2} = \frac{ξ Δ t}{2} \cdot (1 - e^{- 2 ξ Δ t}) - O (ξ^{4} Δ t^{4}), \\ {(1 - e^{- ξ Δ t})}^{2} \approx \frac{ξ Δ t}{2} \cdot (1 - e^{- 2 ξ Δ t}), \end{gathered}

(D2)

and multiplying by $\frac{1}{k_{B} T ξ^{2} m (1 - e^{- 2 ξ Δ t})} {(\nabla U (x_{k}))}^{2}$ yields

\begin{gathered} \frac{1}{k_{B} T ξ^{2} m} \frac{{(1 - e^{- ξ Δ t})}^{2}}{1 - e^{- 2 ξ Δ t}} {(\nabla U (x_{k}))}^{2} \approx \frac{Δ t}{2 k_{B} T ξ m} {(\nabla U (x_{k}))}^{2}, \\ Δ η_{L, k}^{2} \approx Δ η_{o, k}^{2} . \end{gathered}

(D3)

Thus, the difference between $Δ η_{L, k}^{2}$ [Eq. (26)] and $Δ η_{o, k}^{2}$ [Eq. (17)] is of order $O (ξ^{4} Δ t^{4})$ ⁠. Equation (D3) is Eq. (31) squared.□

REFERENCES

1.

I. V.

Gopich

, “

Multisite reversible association in membranes and solutions: From non-Markovian to Markovian kinetics

,”

J. Chem. Phys.

152

,

104101

(

2020

).

https://doi.org/10.1063/1.5144282

Google Scholar

Crossref

PubMed

2.

T. J.

Lane

,

D.

Shukla

,

K. A.

Beauchamp

, and

V. S.

Pande

, “

To milliseconds and beyond: Challenges in the simulation of protein folding

,”

Curr. Opin. Struct. Biol.

23

,

58

(

2013

).

https://doi.org/10.1016/j.sbi.2012.11.002

Google Scholar

Crossref

PubMed

3.

R. O.

Dror

,

R. M.

Dirks

,

J. P.

Grossman

,

H.

Xu

, and

D. E.

Shaw

, “

Biomolecular simulation: A computational microscope for molecular biology

,”

Annu. Rev. Biophys.

41

,

429

(

2012

).

https://doi.org/10.1146/annurev-biophys-042910-155245

Google Scholar

Crossref

PubMed

4.

E. P.

Barros

,

L.

Casalino

,

Z.

Gaieb

,

A. C.

Dommer

,

Y.

Wang

,

L.

Fallon

,

L.

Raguette

,

K.

Belfon

,

C.

Simmerling

, and

R. E.

Amaro

, “

The flexibility of ACE2 in the context of SARS-CoV-2 infection

,”

Biophys. J.

120

,

1

(

2020

).

https://doi.org/10.1016/j.bpj.2020.10.036

Google Scholar

Crossref

PubMed

5.

T. J.

Harpole

and

L.

Delemotte

, “

Conformational landscapes of membrane proteins delineated by enhanced sampling molecular dynamics simulations

,”

Biochim. Biophys. Acta

1860

,

909

(

2018

).

https://doi.org/10.1016/j.bbamem.2017.10.033

Google Scholar

Crossref

6.

Z.

Cournia

,

B.

Allen

, and

W.

Sherman

, “

Relative binding free energy calculations in drug discovery: Recent advances and practical considerations

,”

J. Chem. Inf. Model.

57

,

2911

(

2017

).

https://doi.org/10.1021/acs.jcim.7b00564

Google Scholar

Crossref

PubMed

7.

M.

Badaoui

,

A.

Kells

,

C.

Molteni

,

C. J.

Dickson

,

V.

Hornak

, and

E.

Rosta

, “

Calculating kinetic rates and membrane permeability from biased simulations

,”

J. Phys. Chem. B

122

,

11571

(

2018

).

https://doi.org/10.1021/acs.jpcb.8b07442

Google Scholar

Crossref

PubMed

8.

J.-O.

Joswig

,

J.

Anders

,

H.

Zhang

,

C.

Rademacher

, and

B. G.

Keller

, “

Molecular mechanism of the pH-dependent calcium affinity in Langerin

,” bioRxiv:986851 (

2020

).

Google Scholar

9.

A. S. J. S.

Mey

,

B.

Allen

,

H. E. B.

Macdonald

,

J. D.

Chodera

,

M.

Kuhn

,

J.

Michel

,

D. L.

Mobley

,

L. N.

Naden

,

S.

Prasad

,

A.

Rizzi

,

J.

Scheen

,

M. R.

Shirts

,

G.

Tresadern

, and

H.

Xu

, “

Best practices for alchemical free energy calculations

,”

Living J. Comput. Mol. Sci. Living J. Comp. Mol. Sci. ASAP Version

, pages

2

,

1

, available at https://www.livecomsjournal.org/article/18378-best-practices-for-alchemical-free-energy-calculations-article-v1-0.

10.

M.

Tuckerman

,

Monte Carlo Statistical Mechanics: Theory and Molecular Simulation

(

Oxford University Press, Inc.

,

New York

,

2010

), pp.

300

–

304

.

Google Scholar

11.

D.

Frenkel

and

B.

Smit

,

Understanding Molecular Simulation: From Algorithms to Applications

, 1st ed. (

Academic Press

,

San Diego, San Francisco, NY, Boston, London, Sydney, Tokyo

,

2002

).

Google Scholar

12.

C. A. F.

de Oliveira

,

D.

Hamelberg

, and

J. A.

McCammon

, “

Estimating kinetic rates from accelerated molecular dynamics simulations: Alanine dipeptide in explicit solvent as a case study

,”

J. Chem. Phys.

127

,

175105

(

2007

).

https://doi.org/10.1063/1.2794763

Google Scholar

Crossref

PubMed

13.

P.

Tiwary

and

M.

Parrinello

, “

From metadynamics to dynamics

,”

Phys. Rev. Lett.

111

,

230602

(

2013

).

https://doi.org/10.1103/physrevlett.111.230602

Google Scholar

Crossref

PubMed

14.

O.

Valsson

,

P.

Tiwary

, and

M.

Parrinello

, “

Enhancing important fluctuations: Rare events and metadynamics from a conceptual viewpoint

,”

Annu. Rev. Phys. Chem.

67

,

159

(

2016

).

https://doi.org/10.1146/annurev-physchem-040215-112229

Google Scholar

Crossref

PubMed

15.

R.

Casasnovas

,

V.

Limongelli

,

P.

Tiwary

,

P.

Carloni

, and

M.

Parrinello

, “

Unbinding kinetics of a p38 MAP kinase type II inhibitor from metadynamics simulations

,”

J. Am. Chem. Soc.

139

,

4780

(

2017

).

https://doi.org/10.1021/jacs.6b12950

Google Scholar

Crossref

PubMed

16.

H.

Wu

,

A. S. J. S.

Mey

,

E.

Rosta

, and

F.

Noé

, “

Statistically optimal analysis of state-discretized trajectory data from multiple thermodynamic states

,”

J. Chem. Phys.

141

,

214106

(

2014

).

https://doi.org/10.1063/1.4902240

Google Scholar

Crossref

PubMed

17.

A. S. J. S.

Mey

,

H.

Wu

, and

F.

Noé

, “

xTRAM: Estimating equilibrium expectations from time-correlated simulation data at multiple thermodynamic states

,”

Phys. Rev. X

4

,

041018

(

2014

).

https://doi.org/10.1103/physrevx.4.041018

Google Scholar

Crossref

18.

H.

Wu

,

F.

Paul

,

C.

Wehmeyer

, and

F.

Noé

, “

Multiensemble Markov models of molecular thermodynamics and kinetics

,”

Proc. Natl. Acad. Sci. U. S. A.

113

,

E3221

(

2016

).

https://doi.org/10.1073/pnas.1525092113

Google Scholar

Crossref

PubMed

19.

L. S.

Stelzl

,

A.

Kells

,

E.

Rosta

, and

G.

Hummer

, “

Dynamic histogram analysis to determine free energies and rates from biased simulations

,”

J. Chem. Theory Comput.

13

,

6328

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00373

Google Scholar

Crossref

PubMed

20.

D. J.

Bicout

and

A.

Szabo

, “

Electron transfer reaction dynamics in non-Debye solvents

,”

J. Chem. Phys.

109

,

2325

(

1998

).

https://doi.org/10.1063/1.476800

Google Scholar

Crossref

21.

E.

Rosta

and

G.

Hummer

, “

Free energies from dynamic weighted histogram analysis using unbiased markov state model

,”

J. Chem. Theory Comput.

11

,

276

(

2014

).

https://doi.org/10.1021/ct500719p

Google Scholar

Crossref

PubMed

22.

L.

Donati

,

M.

Heida

,

B. G.

Keller

, and

M.

Weber

, “

Estimation of the infinitesimal generator by square-root approximation

,”

J. Phys.: Condens. Matter

30

,

425201

(

2018

).

https://doi.org/10.1088/1361-648x/aadfc8

Google Scholar

Crossref

PubMed

23.

S.

Kieninger

,

L.

Donati

, and

B. G.

Keller

, “

Dynamical reweighting methods for Markov models

,”

Curr. Opin. Struct. Biol.

61

,

124

(

2020

).

https://doi.org/10.1016/j.sbi.2019.12.018

Google Scholar

Crossref

PubMed

24.

D. M.

Zuckerman

and

T. B.

Woolf

, “

Dynamic reaction paths and rates through importance-sampled stochastic dynamics

,”

J. Chem. Phys.

111

,

9475

(

1999

).

https://doi.org/10.1063/1.480278

Google Scholar

Crossref

25.

T. B.

Woolf

, “

Path corrected functionals of stochastic trajectories: Towards relative free energy and reaction coordinate calculations

,”

Chem. Phys. Lett.

289

,

433

(

1998

).

https://doi.org/10.1016/s0009-2614(98)00427-8

Google Scholar

Crossref

26.

D. M.

Zuckerman

and

T. B.

Woolf

, “

Efficient dynamic importance sampling of rare events in one dimension

,”

Phys. Rev. E

63

,

016702

(

2000

).

https://doi.org/10.1103/physreve.63.016702

Google Scholar

Crossref

27.

C.

Xing

and

I.

Andricioaei

, “

On the calculation of time correlation functions by potential scaling

,”

J. Chem. Phys.

124

,

034110

(

2006

).

https://doi.org/10.1063/1.2159476

Google Scholar

Crossref

PubMed

28.

A. B.

Adib

, “

Stochastic actions for diffusive dynamics: Reweighting, sampling, and minimization

,”

J. Phys. Chem. B

112

,

5910

(

2008

).

https://doi.org/10.1021/jp0751458

Google Scholar

Crossref

PubMed

29.

I. V.

Girsanov

, “

On transforming a certain class of stochastic processes by absolutely continuous substitution of measures

,”

Theory Probab. Appl.

5

,

285

(

1960

).

https://doi.org/10.1137/1105027

Google Scholar

Crossref

30.

B.

Øksendal

,

Stochastic Differential Equations: An Introduction with Applications

, 6th ed. (

Springer Verlag

,

Berlin

,

2003

).

Google Scholar

Crossref

31.

L.

Onsager

and

S.

Machlup

, “

Fluctuations and irreversible processes

,”

Phys. Rev.

91

,

1505

(

1953

).

https://doi.org/10.1103/physrev.91.1505

Google Scholar

Crossref

32.

J.-H.

Prinz

,

J. D.

Chodera

,

V. S.

Pande

,

W. C.

Swope

,

J. C.

Smith

, and

F.

Noé

, “

Optimal use of data in parallel tempering simulations for the construction of discrete-state Markov models of biomolecular dynamics

,”

J. Chem. Phys.

134

,

244108

(

2011

).

https://doi.org/10.1063/1.3592153

Google Scholar

Crossref

PubMed

33.

C.

Schütte

,

A.

Nielsen

, and

M.

Weber

, “

Markov state models and molecular alchemy

,”

Mol. Phys.

113

,

69

(

2015

).

https://doi.org/10.1080/00268976.2014.944597

Google Scholar

Crossref

34.

L.

Donati

,

C.

Hartmann

, and

B. G.

Keller

, “

Girsanov reweighting for path ensembles and Markov state models

,”

J. Chem. Phys.

146

,

244112

(

2017

).

https://doi.org/10.1063/1.4989474

Google Scholar

Crossref

PubMed

35.

L.

Donati

and

B. G.

Keller

, “

Girsanov reweighting for metadynamics simulations

,”

J. Chem. Phys.

149

,

072335

(

2018

).

https://doi.org/10.1063/1.5027728

Google Scholar

Crossref

PubMed

36.

W.

Huisinga

,

C.

Schütte

, and

A. M.

Stuart

, “

Extracting macroscopic stochastic dynamics: Model problems

,”

Commun. Pure Appl. Math.

56

,

234

(

2003

).

https://doi.org/10.1002/cpa.10057

Google Scholar

Crossref

37.

W. C.

Swope

,

J. W.

Pitera

, and

F.

Suits

, “

Describing protein folding kinetics by molecular dynamics simulations. 1. Theory

,”

J. Phys. Chem. B

108

,

6571

(

2004

).

https://doi.org/10.1021/jp037421y

Google Scholar

Crossref

38.

N.-V.

Buchete

and

G.

Hummer

, “

Coarse master equations for peptide folding dynamics

,”

J. Phys. Chem. B

112

,

6057

(

2008

).

https://doi.org/10.1021/jp0761665

Google Scholar

Crossref

PubMed

39.

B.

Keller

,

X.

Daura

, and

W. F.

van Gunsteren

, “

Comparing geometric and kinetic cluster algorithms for molecular simulation data

,”

J. Chem. Phys.

132

,

074110

(

2010

).

https://doi.org/10.1063/1.3301140

Google Scholar

Crossref

PubMed

40.

J.-H.

Prinz

,

H.

Wu

,

M.

Sarich

,

B. G.

Keller

,

M.

Senne

,

M.

Held

,

J. D.

Chodera

,

C.

Schütte

, and

F.

Noé

, “

Markov models of molecular kinetics: Generation and validation

,”

J. Chem. Phys.

134

,

174105

(

2011

).

https://doi.org/10.1063/1.3565032

Google Scholar

Crossref

PubMed

41.

J.-H.

Prinz

,

B.

Keller

, and

F.

Noé

, “

Probing molecular kinetics with Markov models: Metastable states, transition pathways and spectroscopic observables

,”

Phys. Chem. Chem. Phys.

13

,

16912

(

2011

).

https://doi.org/10.1039/c1cp21258c

Google Scholar

Crossref

PubMed

42.

B. E.

Husic

and

V. S.

Pande

, “

Markov state models: From an art to a science

,”

J. Am. Chem. Soc.

140

,

2386

(

2018

).

https://doi.org/10.1021/jacs.7b12191

Google Scholar

Crossref

PubMed

43.

J. D.

Chodera

,

W. C.

Swope

,

F.

Noé

,

J.-H.

Prinz

,

M. R.

Shirts

, and

V. S.

Pande

, “

Dynamical reweighting: Improved estimates of dynamical properties from simulations at multiple temperatures

,”

J. Chem. Phys.

134

,

244107

(

2011

).

https://doi.org/10.1063/1.3592152

Google Scholar

Crossref

PubMed

44.

W. F.

van Gunsteren

and

H. J. C.

Berendsen

, “

Algorithms for Brownian dynamics

,”

Mol. Phys.

45

,

637

(

1981

).

https://doi.org/10.1080/00268978200100491

Google Scholar

Crossref

45.

A.

Brünger

,

C. L.

Brooks

, and

M.

Karplus

, “

Stochastic boundary conditions for molecular dynamics simulations of ST2 water

,”

Chem. Phys. Lett.

105

,

495

(

1984

).

https://doi.org/10.1016/0009-2614(84)80098-6

Google Scholar

Crossref

46.

G.

Stoltz

, “

Path sampling with stochastic dynamics: Some new algorithms

,”

J. Comput. Phys.

225

,

491

(

2007

).

https://doi.org/10.1016/j.jcp.2006.12.006

Google Scholar

Crossref

47.

G.

Bussi

and

M.

Parrinello

, “

Accurate sampling using Langevin dynamics

,”

Phys. Rev. E

75

,

056707

(

2007

).

https://doi.org/10.1103/physreve.75.056707

Google Scholar

Crossref

48.

M.

Ceriotti

,

G.

Bussi

, and

M.

Parrinello

, “

Langevin equation with colored noise for constant-temperature molecular dynamics simulations

,”

Phys. Rev. Lett.

102

,

020601

(

2009

).

https://doi.org/10.1103/physrevlett.102.020601

Google Scholar

Crossref

PubMed

49.

J. A.

Izaguirre

,

C. R.

Sweet

, and

V.

Pande

, “

Multiscale dynamics of macromolecules using normal mode Langevin

,”

Pac. Symp. Biocomput.

15

,

240

(

2010

).

https://doi.org/10.1142/9789814295291_0026

Google Scholar

Crossref

50.

N.

Goga

,

A. J.

Rzepiela

,

A. H.

de Vries

,

S. J.

Marrink

, and

H. J. C.

Berendsen

, “

Efficient algorithms for Langevin and DPD dynamics

.”

J. Chem. Theory Comput.

8

,

3637

(

2012

).

https://doi.org/10.1021/ct3000876

Google Scholar

Crossref

PubMed

51.

B.

Leimkuhler

and

C.

Matthews

, “

Robust and efficient configurational molecular sampling via Langevin dynamics

,”

J. Chem. Phys.

138

,

174102

(

2013

).

https://doi.org/10.1063/1.4802990

Google Scholar

Crossref

PubMed

52.

D. A.

Sivak

,

J. D.

Chodera

, and

G. E.

Crooks

, “

Time step rescaling recovers continuous-time dynamical properties for discrete-time Langevin integration of nonequilibrium systems

,”

J. Phys. Chem. B

118

,

6466

(

2014

).

https://doi.org/10.1021/jp411770f

Google Scholar

Crossref

PubMed

53.

J.

Fass

,

D.

Sivak

,

G.

Crooks

,

K.

Beauchamp

,

B.

Leimkuhler

, and

J.

Chodera

, “

Quantifying configuration-sampling error in Langevin simulations of complex molecular systems

,”

Entropy

20

,

318

(

2018

).

https://doi.org/10.3390/e20050318

Google Scholar

Crossref

PubMed

54.

P.

Eastman

,

J.

Swails

,

J. D.

Chodera

,

R. T.

McGibbon

,

Y.

Zhao

,

K. A.

Beauchamp

,

L.-P.

Wang

,

A. C.

Simmonett

,

M. P.

Harrigan

,

C. D.

Stern

,

R. P.

Wiewiora

,

B. R.

Brooks

, and

V. S.

Pande

, “

OpenMM 7: Rapid development of high performance algorithms for molecular dynamics

,”

PLoS Comput. Biol.

13

,

1

(

2017

).

https://doi.org/10.1371/journal.pcbi.1005659

Google Scholar

Crossref

55.

N.

Bou-Rabee

, “

Time integrators for molecular dynamics

,”

Entropy

16

,

138

(

2014

).

https://doi.org/10.3390/e16010138

Google Scholar

Crossref

56.

C. C.

Chow

and

M. A.

Buice

, “

Path integral methods for stochastic differential equations

,”

J. Math. Neurosci.

5

,

1

(

2015

).

https://doi.org/10.1186/s13408-015-0018-5

Google Scholar

Crossref

PubMed

57.

P. C.

Bressloff

,

Stochastic Processes in Cell Biology

, 1st ed. (

Springer

,

New York

,

2014

).

Google Scholar

Crossref

58.

P. H.

Hünenberger

, “

Thermostat algorithms for molecular dynamics simulations

,” in

Advanced Computer Simulation

(

Springer

,

Berlin, Heidelberg

,

2005

).

Google Scholar

Crossref

59.

J. E.

Basconi

and

M. R.

Shirts

, “

Effects of temperature control algorithms on transport properties and kinetics in molecular dynamics simulations

,”

J. Chem. Theory Comput.

9

,

2887

(

2013

).

https://doi.org/10.1021/ct400109a

Google Scholar

Crossref

PubMed

60.

J.

Wang

,

R. M.

Wolf

,

J. W.

Caldwell

,

P. A.

Kollman

, and

D. A.

Case

, “

Development and testing of a general amber force field

,”

J. Comput. Chem.

25

,

1157

(

2004

).

https://doi.org/10.1002/jcc.20035

Google Scholar

Crossref

PubMed

61.

A.

Onufriev

,

D.

Bashford

, and

D. A.

Case

, “

Exploring protein native states and large-scale conformational changes with a modified generalized born model

,”

J. Comput. Chem.

55

,

383

(

2004

).

https://doi.org/10.1002/prot.20033

Google Scholar

Crossref

62.

See http://docs.openmm.org/latest/api-python/generated/simtk.openmm.openmm.CustomIntegrator.html for information about the CustomIntegrator Class of the simulation package OpenMM; accessed 25 January 2021.

63.

M.

Bause

,

T.

Wittenstein

,

K.

Kremer

, and

T.

Bereau

, “

Microscopic reweighting for nonequilibrium steady-state dynamics

,”

Phys. Rev. E

100

,

060103

(

2019

).

https://doi.org/10.1103/physreve.100.060103

Google Scholar

Crossref

PubMed

64.

G.

Pérez-Hernández

,

F.

Paul

,

T.

Giorgino

,

G.

De Fabritiis

, and

F.

Noé

, “

Identification of slow molecular order parameters for Markov model construction

,”

J. Chem. Phys.

139

,

015102

(

2013

).

https://doi.org/10.1063/1.4811489

Google Scholar

Crossref

PubMed

65.

F.

Nüske

,

B. G.

Keller

,

G.

Pérez-Hernández

,

A. S. J. S.

Mey

, and

F.

Noé

, “

Variational approach to molecular kinetics

,”

J. Chem. Theory Comput.

10

,

1739

–

1752

(

2014

).

https://doi.org/10.1021/ct4009156

Google Scholar

Crossref

PubMed

66.

O.

Lemke

and

B. G.

Keller

, “

Density-based cluster algorithms for the identification of core sets

,”

J. Chem. Phys.

145

,

164104

(

2016

).

https://doi.org/10.1063/1.4965440

Google Scholar

Crossref

PubMed

67.

L. T.

Chong

,

A. S.

Saglam

, and

D. M.

Zuckerman

, “

Path-sampling strategies for simulating rare events in biomolecular systems

,”

Curr. Opin. Struct. Biol.

43

,

88

–

94

(

2017

).

https://doi.org/10.1016/j.sbi.2016.11.019

Google Scholar

Crossref

PubMed

68.

G.

Grazioli

and

I.

Andricioaei

, “

Advances in milestoning. I. Enhanced sampling via wind-assisted reweighted milestoning (WARM)

,”

J. Chem. Phys.

149

,

084103

(

2018

).

https://doi.org/10.1063/1.5029954

Google Scholar

Crossref

PubMed

69.

P. D.

Dixit

,

J.

Wagoner

,

C.

Weistuch

,

S.

Pressé

,

K.

Ghosh

, and

K. A.

Dill

, “

Perspective: Maximum caliber is a general variational principle for dynamical systems

,”

J. Chem. Phys.

148

,

010901

(

2018

).

https://doi.org/10.1063/1.5012990

Google Scholar

Crossref

PubMed

70.

E. K.

Peter

,

J.-E.

Shea

, and

A.

Schug

, “

CORE-MD, a path correlated molecular dynamics simulation method

,”

J. Chem. Phys.

153

,

084114

(

2020

).

https://doi.org/10.1063/5.0015398

Google Scholar

Crossref

PubMed

71.

See https://github.com/openmm/openmm/blob/master/platforms/cpu/src/CpuLangevinDynamics.cpp for information about the CpuLangevinDynamics Class of the simulation package OpenMM accessed 15 November 2020.

72.

C. W.

Gardiner

,

Handbook of Stochastic Methods for Physics, Chemistry and the Natural Sciences

, 2nd ed. (

Springer Verlag

,

Berlin, Heidelberg

,

1983

).

Google Scholar

Crossref

2021

Author(s)

Path probability ratios for Langevin dynamics—Exact and approximate

I. INTRODUCTION

II. PATH REWEIGHTING

III. OVERDAMPED LANGEVIN DYNAMICS

A. Equation of motion and integration scheme

B. Path probability and path probability ratio

C. Random numbers and random number probability ratio

IV. LANGEVIN DYNAMICS

A. Equation of motion and integration scheme

B. Path probability and path probability ratio

C. Random numbers and random number probability ratio

V. COMPARING LANGEVIN AND OVERDAMPED LANGEVIN DYNAMICS

A. Test system

B. From random numbers η to paths ωo and ωL

C. From a path ω to random numbers ηo and ηL

D. Path reweighting

VI. APPROXIMATE PATH PROBABILITY RATIO

A. Derivation and numerical results

B. Intuition

VII. MOLECULAR EXAMPLE: BUTANE

VIII. METHODS

A. Simulations of the test system

B. Butane—Direct simulations

C. Butane—Path reweighting

IX. CONCLUSION AND OUTLOOK

SUPPLEMENTARY MATERIAL

DEDICATION

ACKNOWLEDGMENTS

DATA AVAILABILITY

APPENDIX A: LANGEVIN LEAPFROG AND THE ISP SCHEME

APPENDIX B: PATH PROBABILITY FOR LANGEVIN DYNAMICS

APPENDIX C: SOLVING THE DOUBLE INTEGRAL

APPENDIX D: PROOF OF EQ. (31)

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

B. From random numbers η to paths ω_o and ω_L

C. From a path ω to random numbers η_o and η_L