Girsanov reweighting for path ensembles and Markov state models

TL;DR

This paper introduces a Girsanov reweighting method to estimate how molecular dynamics path ensembles and Markov state models respond to potential energy perturbations, enabling efficient sensitivity analysis and trajectory reweighting.

Contribution

The paper develops a Girsanov reweighting approach for path ensembles and MSMs, extending stochastic analysis tools to molecular dynamics simulations for the first time.

Findings

Effective reweighting of path ensembles demonstrated on test systems.

Method enables on-the-fly reweighting during simulations.

Applicable to both implicit and explicit solvent models.

Abstract

The sensitivity of molecular dynamics on changes in the potential energy function plays an important role in understanding the dynamics and function of complex molecules.We present a method to obtain path ensemble averages of a perturbed dynamics from a set of paths generated by a reference dynamics. It is based on the concept of path probability measure and the Girsanov theorem, a result from stochastic analysis to estimate a change of measure of a path ensemble. Since Markov state models (MSM) of the molecular dynamics can be formulated as a combined phase-space and path ensemble average, the method can be extended toreweight MSMs by combining it with a reweighting of the Boltzmann distribution. We demonstrate how to efficiently implement the Girsanov reweighting in a molecular dynamics simulation program by calculating parts of the reweighting factor "on the fly" during the simulation, and we benchmark the method on test systems ranging from a two-dimensional diffusion process to an artificial many-body system and alanine dipeptide and valine dipeptide in implicit and explicit water. The method can be used to study the sensitivity of molecular dynamics on external perturbations as well as to reweight trajectories generated by enhanced sampling schemes to the original dynamics.

Tables

Table 1. Table 1: Overview of the notation

	State: $𝐱$	Path: $\omega$
Space	$\Gamma \subset {𝐑}^{6N}$	${\Omega }_{\tau ,𝐱}={\Gamma }^n\subset {𝐑}^{6N\cdot n}$
Subsets	$A_i\subset \Gamma$	$𝒜=\left({\prod }_{i=1}^n{A}_i\right)\subset {\Omega }_{\tau ,𝐱}$
Probability density	${{\mu }_{\pi }\left(𝐱\right)|}_{t\to \infty }$	${\mu }_P\left(\omega \right)$
Probability of a state/path	$\pi \left(A\right)={\int }_A{\mu }_{\pi }\left(𝐱\right)d𝐱$	$P\left(𝒜\right)={\int }_{𝒜}{\mu }_P\left(\omega \right)d\omega$
Expectation value	${𝔼}_{\pi }\left[f\left(𝐱\right)\right]={\int }_{A}f\left(𝐱\right){\mu }_{\pi }\left(𝐱\right)d𝐱$	${𝔼}_P\left[f\left(\omega \right)\right]={\int }_{𝒜}f\left(\omega \right){\mu }_P\left(\omega \right)d\omega$
Absolute continuity	$\tilde{\pi }\left(A\right)=0\Rightarrow \pi \left(A\right)=0$	$\tilde{P}\left(𝒜\right)=0\Rightarrow P\left(𝒜\right)=0$
Radon-Nikodyn derivative	$\frac{d\tilde{\pi }}{d\pi }$, see eq. 31	$\frac{d\tilde{P}}{dP}$, see eq. 36

Equations

$$M\frac{\mathrm{d}\mathbf{v}(t)}{\mathrm{d}t}$$

$$\mathbf{v}(t)$$

$$\sigma =\sqrt{2k_{B}T\gamma M}$$

$${\mu }_{\pi }(\mathbf{x})=\frac{\exp \left[-\beta \mathcal{H}(\mathbf{x})\right]}{Z}$$

$$\pi (A)=\mathbb{P}(\mathbf{x}\in A)={\int }_A{\mu }_{\pi }(\mathbf{x})\mathrm{d}\mathbf{x},\forall A\subset \Gamma$$

$${\mathbb{E}}_{\pi }\left[a\right]$$

$$P(\mathcal{A})$$

$$=$$

$${\mu }_P(\omega )$$

$$P(\mathcal{A})$$

$${\mathbb{E}}_P[f(\omega )]$$

$$=$$

$$p(\mathbf{x},\mathbf{y},\tau )={\int }_{\Gamma }{\int }_{\Gamma }\cdots {\int }_{\Gamma }{\mu }_P({\mathbf{x}}_1,{\mathbf{x}}_2,...,\mathbf{y})\mathrm{d}{\mathbf{x}}_1\mathrm{d}{\mathbf{x}}_2...\mathrm{d}{\mathbf{x}}_{n-1}.$$

$$1_{B_i}(\mathbf{x}):=\{\begin{array}{ll}1& \text{if }\mathbf{x}\in B_i,\\ 0& \text{mbox}otherwise.\end{array}$$

$$C_{ij}(\tau )$$

$$T_{ij}(\tau )=\frac{C_{ij}(\tau )}{{\sum }_{j=1}^s{C}_{ij}(\tau )}.$$

$${\int }_{\Gamma }p(\mathbf{x},\mathbf{y};\tau )1_{B_j}(\mathbf{y})\mathrm{d}\mathbf{y}$$

$$=$$

$$C_{ij}(\tau )$$

$$=$$

$${\mathbb{E}}_{P,\pi }[f(\omega )]$$

$$C_{ij}(\tau )$$

$$=$$

$$M\frac{\mathrm{d}\mathbf{v}(t)}{\mathrm{d}t}$$

$$\mathbf{v}(t)$$

$${\mu }_{\tilde{\pi }}(\mathbf{x})=\frac{\exp \left[-\beta \tilde{\mathcal{H}}(\mathbf{x})\right]}{\tilde{Z}},$$

$$\displaystyle\widetilde{\pi}(A)=\int_{A}\mu_{\widetilde{\pi}}(\mathbf{x})\,\mathrm{d}\mathbf{x}=0\$$

$$g(\mathbf{x})=\frac{\mathrm{d}\tilde{\pi }}{\mathrm{d}\pi }=\frac{{\mu }_{\tilde{\pi }}(\mathbf{x})}{{\mu }_{\pi }(\mathbf{x})}=\frac{Z}{\tilde{Z}}\exp \left(-\beta U(\mathbf{x})\right).$$

$$\tilde{\pi }(A)={\int }_A{\mu }_{\tilde{\pi }}(\mathbf{x})\mathrm{d}\mathbf{x}={\int }_{A}g(\mathbf{x}){\mu }_{\pi }(\mathbf{x})\mathrm{d}\mathbf{x}.$$

$${\mathbb{E}}_{\tilde{\pi }}\left[a\right]$$

$$=$$

$$\tilde{P}(\mathcal{A})={\int }_{\mathcal{A}}{\mu }_{\tilde{P}}(\omega )\mathrm{d}\omega =0$$

$$M_{\tau ,\mathbf{x}}(\omega )=\frac{{\mu }_{\tilde{P}}(\omega )}{{\mu }_P(\omega )}=\exp \left\{-\sum _{i=1}^{3N}\left[\sum _{k=0}^n\frac{{\nabla }_{i}U({\mathbf{r}}_k)}{\sigma }{\eta }_k^i\sqrt{\Delta t}-\frac{1}{2}\sum _{k=0}^n{\left(\frac{{\nabla }_{i}U({\mathbf{r}}_k)}{\sigma }\right)}^2\Delta t\right]\right\}$$

Full text

Markov State Models, Dynamical Reweighting, Path Ensembles, Molecular Dynamics

††preprint: AIP/123-QED

I Introduction

Molecular dynamics (MD) simulations with explicit solvent molecules are routinely used as efficient importance sampling algorithms for the Boltzmann distribution of molecular systems. From the conformational snapshots created by MD simulations, one can estimate phase-space ensemble averages and thus interpret experimental data or thermodynamic functions in terms of molecular conformations.

In recent years, the scope of MD simulations has considerably widened, and the method has been increasingly used to construct models of the conformational dynamics Bolhuis2002 ; Faradjian2004 ; Best2005 ; Hegger2009 ; Faccioli2010 . Most notably Markov state models Schuette1999 ; Schuette1999b ; Deuflhard2000 ; Swope2004 ; Chodera2007 ; Buchete2008 ; Prinz2011 ; Keller2010 ; Keller2011 , in which the conformational space is discretized into disjoint states and the dynamics is modeled as a Markov jump process between these states, have become a valuable tool for the analysis of complex molecular dynamics Voelz2010b ; Stanley2014 ; Bowman2015 ; Plattner2015 ; Zhang2016 .

For the construction of dynamic models, one has to estimate path ensemble averages. For example in MSMs, the transition probabilities between a pair of states $B_{i}$ and $B_{j}$ is estimated by considering a set of paths $S_{\tau}=\{\omega_{1},\omega_{2},...\omega_{n}\}$, each of which has length $\tau$, counting the number of paths which start in $B_{i}$ and end in $B_{j}$, and comparing this number to the total number of paths in the set (vertical line of blue boxes in Fig. 1).

Suppose, one would like to compare the dynamics in a reference potential energy function $V(\mathbf{r})$ to the dynamics in a series of perturbed potential energy functions $V(\mathbf{r})+U(\mathbf{r},\kappa)$, where $U(\mathbf{r},\kappa)$ represents the perturbation and $\kappa$ is a tunable parameter, e.g. a force constant. While for phase-space ensemble averages, numerous methods exist to reweight the samples of the reference conformational ensemble to yield ensemble averages for the perturbed systems Zwanzig1954 ; Bennett1976 ; Kumar1992 ; Kaestner2005 , similar reweighting schemes have not yet been developed for path ensemble averages of explicit-solvent simulations. This means that currently one would have to re-simulate the dynamics at each parameter value (vertical lines of green boxes in Fig. 1) and then construct a MSM for each simulation separately. This is computationally extremely costly.

An alternative would be to reweight the path ensemble average at the reference potential energy function to obtain path ensemble averages for the perturbed systems (reweighting box in Fig. 1). From measure theory it is well known that a reweighting factor is given as the ratio between the probability measure associated to the reference potential energy function and the probability measure associated to the perturbed potential energy function. This applies to reweighting phase space ensemble averages as well as to reweighting path ensemble averages. Fig. 2 illustrates the idea of a path ensemble reweighting method. The figure shows two sets of paths, one generated by a Brownian dynamics simulation without drift, $S_{\tau}$ (Fig. 2.A), and one generated by a Brownian dynamics simulation with drift, $\widetilde{S}_{\tau}$ (Fig. 2.B). Both simulations sample the same path space $\Omega_{\tau,\mathbf{x}}$ but the probability with which a given path is realized differs in the two simulations. In Fig. 2.A and 2.B, the sets of paths, $S_{\tau}$ and $\widetilde{S}_{\tau}$, are colored according to their respective path probability density $\mu_{P}(\omega)$ and $\mu_{\widetilde{P}}(\omega)$. Fig. 2.C and 2.D shows again $S_{\tau}$, this time however we colored the paths according to the probability density $\mu_{\widetilde{P}}(\omega)$ with which they would have been generated by a Brownian dynamics with drift. For Brownian dynamics this probability can be calculated directly (Fig. 2.C), for other types of dynamics a reweighting method has to be used (Fig. 2.D). To estimate a path ensemble average for the Brownian dynamics with drift from $S_{\tau}$, the contribution of each path $\omega$ to the estimate is multiplied by the ratio $M_{\tau,\mathbf{x}}(\omega)=\mu_{\widetilde{P}}(\omega)/\mu_{P}(\omega)$ (Fig. 1).

Path ensemble reweighting schemes have initially been developed in the field of importance sampling for stochastic differential equations LPS98 ; PR09 ; HS12 . For Langevin dynamics, the Girsanov theorem Girsanov1960 ; Oksendal2003 provides us with an expression for the probability ratio, and thus reweighting path ensemble averages becomes possible for this type of dynamics (section II.4). It has recently been demonstrated that the theorem can be applied to reweight Markov state models of Brownian dynamics in one- and two-dimensional potential energy functions Schuette2015b .

Here we demonstrate how the Girsanov reweighting scheme can be applied to explicit-solvent all-atom MD simulations. For this, we need to address to critical pillars on which the Girsanov reweighting scheme rests:

• •

The equation of motion need to contain a stochastic term which generates random forces drawn from a normal distribution (white noise).

• •

To calculate the reweighting factor, the random forces need to be accessible for each MD simulation step.

In all-atom MD simulations, the system is propagated by the Newton equations of motion which do not contain a stochastic term. We will discuss, how the Girsanov theorem can nonetheless be applied to this type of simulations (section IV.2). The second point, in principle, requires that the forces are written out at every MD simulation step, i.e. at a frequency of femtoseconds rather than the usual output rate of several picoseconds. This quickly fills up any hard disc and slows the simulation by orders of magnitudes. We will present a computational efficient implementation of the reweighting scheme (section IV.1).

When applying the Girsanov theorem to reweight an MSM, an additional difficulty arises

• •

The degrees of freedom which are affected by the perturbation might not be part of the relevant subspace of the MSM.

We found that the reweighting becomes problematic in this case and propose to project the perturbed degrees of freedom onto the relevant subspace during the estimation of the ratio of probability measures (section II.5).

The Girsanov reweighting method is demonstrated and benchmarked on several systems, ranging from two-dimensional diffusion processes (section IV.3), over molecular model systems which follow a Langevin dynamics (section IV.4) to all-atom MD simulations of alanine and valine dipeptides in explicit and implicit solvent (section IV.5).

II Theory

II.1 Molecular dynamics

Consider a molecular system with $N$ particles, which evolves in time $t$ according to the Langevin equation

[TABLE]

where $M$ is the mass matrix, $\mathbf{r}(t)$ and $\mathbf{v}(t)\in\mathbb{R}^{3N}$ are the position vector and the velocity vector. $V(\mathbf{r})$ is the potential energy function. The interaction with the thermal bath is modelled by the friction coefficient $\gamma$, and an uncorrelated Gaussian white noise $\eta(t)\in\mathbb{R}^{3N}$ which is scaled by the volatility $\sigma$

[TABLE]

where $k_{B}$ is the Boltzmann constant and $T$ is the temperature of the system.

The phase-space vector $\mathbf{x}(t)=\{\mathbf{r}(t),\mathbf{v}(t)\}\in\Gamma$ fully represents the state of the system at time $t$, where $\Gamma=\mathbb{R}^{6N}$ denotes the phase space of the system. The dynamics in eq. 1 is associated to an equilibrium probability density

[TABLE]

where $\beta=\frac{1}{k_{B}T}$, $\mathcal{H}(\mathbf{x})=\frac{1}{2}\mathbf{v}^{\top}M\mathbf{v}+V(\mathbf{r})$ is the classical Hamiltonian of the system, and $Z=\int_{\Gamma}\exp\left[-\beta\mathcal{H}(\mathbf{x})\right]\,\mathrm{d}\mathbf{x}$ is the partition function. The function $\mu_{\pi}(\mathbf{x})$ is associated to the probability measure

[TABLE]

where $\pi(A)$ represents the equilibrium probability of finding the system in a subset $A$ of the phase space $\Gamma$. The expectation value of a function $a(\mathbf{x}):\Gamma\rightarrow\mathbb{R}$ with respect to the probability density $\mu_{\pi}(\mathbf{x})$ is given as

[TABLE]

where $S_{n}=\{\mathbf{x}_{1},...\mathbf{x}_{n}\}$ is a set of states distributed according to eq. 4. When the phase space vectors $\mathbf{x}_{i}$ are generated by numerically integrating eq. 1, the second equality only holds if the sampling is ergodic. Eq. 6 defines a phase-space ensemble average. The subscript $\pi$ indicates the measure for which the expectation value is calculated.

II.2 Path ensembles and MSMs

A path $\omega=\{\mathbf{x}(t=0)=\mathbf{x}_{0},\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}(\tau)=\mathbf{x}_{n}\}$ is a time-discretized realization of the dynamics $\mathbf{x}(t)$ on the time interval $[0,\tau=n\cdot\Delta t]$ starting at particular point $\mathbf{x}_{0}\in\Gamma$, where $\Delta t$ is the time step and $n\in\mathbb{N}$ is the number of time steps. The associated path space is denoted $\Omega_{\tau,\mathbf{x}}=\mathbb{R}^{6N\cdot n}$. A subset of the path space $\mathcal{A}$ is constructed as a product of subsets $A_{i}\subset\Gamma$ of the state space $\mathcal{A}=A_{1}\times A_{2}...\times A_{n}$, where the subset $A_{i}$ represents the phase space volume in which $\mathbf{x}_{i}$ may be found. The probability that by integrating eq. 1 one obtains a path $\omega$ which belongs to the subset $\mathcal{A}\subset\Omega_{\tau,\mathbf{x}}$ is given as

[TABLE]

The function $p(\mathbf{x}_{i},\mathbf{x}_{i+1};\,\Delta t)$ is the transition probability density , i.e. the conditional probability to be in $\mathbf{x}_{i+1}$ after a time $\Delta t$ given the initial state $\mathbf{x}_{i}$. The function $P$ is a path probability measure and is the analogon to $\pi$ in phase space ensemble averages (eq. 5). The path probability measure is associated to the path probability density function:

[TABLE]

and hence the formal analogon to eq. 5 in the path space is

[TABLE]

where the integration over $\mathrm{d}\omega$ is defined by eq. 7.

Let $f:\Omega_{\tau,\mathbf{x}}\rightarrow\mathbb{R}$ be an integrable function, which assigns a real number to each path. The expectation value of this function is

[TABLE]

where we again assumed that the paths have a common initial state $\mathbf{x}(t=0)=\mathbf{x}_{0}$, and $S_{\tau,\mathbf{x},m}=\{\omega_{1},\omega_{2},...\omega_{m}\}$ corresponds to a set of paths of length $\tau$ generated by numerically integrating eq. 1. When the paths are extracted from a single long trajectory, the last equality only holds if the sampling is ergodic. Eq. 13 defines a path ensemble average. The subscript $P$ indicates that the expectation value is calculated with respect to a path probability measure.

For Markov processes, one can define a transition probability density $p(\mathbf{x},\mathbf{y};\tau)$, i.e. the conditional probability to be in $\mathbf{x}_{n}=\mathbf{y}$ after a time $\tau$ given that the path started in $\mathbf{x}_{0}=\mathbf{x}$, by integrating the path probability density over all intervening states and applying recursively the Chapman-Kolmogorov equation

[TABLE]

Markov processes can be approximated by Markov state models Schuette1999 ; Schuette1999b ; Deuflhard2000 ; Swope2004 ; Chodera2007 ; Buchete2008 ; Prinz2011 ; Keller2011 . In these models, the phase space is discretized into disjoint sets (or microstates) ${B_{1},B_{2},...B_{s}}$ with $\cup_{i=1}^{m}B_{i}=\Gamma$, where the indicator function of the $i$th state is given by

[TABLE]

The associated cross-correlation function is

[TABLE]

and the transition probability between set $B_{i}$ and set $B_{j}$ is

[TABLE]

$T_{ij}(\tau)$ are the elements of the transition matrix whose dominant eigenvectors and eigenvalues represent the slow dynamic processes of the system Schuette1999 ; Schuette1999b ; Deuflhard2000 ; Swope2004 ; Chodera2007 ; Buchete2008 ; Prinz2011 ; Keller2011 .

Because of eq. 14, one can regard the inner integral of the cross-correlation function $C_{ij}(\tau)$ (eq. 16) as a path ensemble average (eq. 13) which depends on the initial state $\mathbf{x}_{0}=\mathbf{x}$ of the path ensemble

[TABLE]

This expected value is a linear operator and it is called backward transfer operator (more details are given in appendix). The outer integral in eq. 16 is a phase space ensemble average, thus we have

[TABLE]

where the combined phase space and ensemble average is defined as

[TABLE]

Eq. 24 extends eqs. 13 to path ensembles with arbitrary initial states. The elements $C_{ij}(\tau)$ can be estimated from a set of paths of length $\tau$, $S_{\tau,m}=\{\nu_{1},\nu_{2}...\nu_{m}\}$, in which the initial states are no longer fixed but are distributed according to $\mu_{\pi}(\mathbf{x})$,

[TABLE]

where $[\nu_{k}]_{t=i}=[\mathbf{x}_{i}]_{k}$ denotes the $i$th time step of the $k$th path.

II.3 Dynamical reweighting

We alter the reference dynamics (eq. 1) by adding a perturbation $U(\mathbf{r}):\Gamma\rightarrow\mathbb{R}$ to the potential energy function $V(\mathbf{r})$. Thus, $\widetilde{V}(\mathbf{r})=V(\mathbf{r})+U(\mathbf{r})$ is the perturbed potential energy function and the Langevin equations of motion are

[TABLE]

The perturbed dynamics is associated to a peturbed stationary probability density

[TABLE]

where $\widetilde{\mathcal{H}}(\mathbf{x})=\frac{1}{2}\mathbf{v}^{\top}M\mathbf{v}+\widetilde{V}(\mathbf{r})$ is the Hamiltonian of the perturbed system, and $\widetilde{Z}$ is its partition function. The perturbation of the potential energy function also changes the transition probability density $\widetilde{p}(\mathbf{x}_{i},\mathbf{x}_{i+1};\,\Delta t)$, which gives rise to a perturbed path probability density $\mu_{\widetilde{P}}(\omega)$ (eq. 9) and a perturbed path measure $\widetilde{P}$ (eq. 7) .

Reweighting methods compare the probability measure of the perturbed systems to the probability measure of a reference system. We first review the derivation of a reweighting scheme for phase space probability measures before discussing path space probability measures. The perturbed phase-space probability measure $\widetilde{\pi}$ is said to be absolutely continuous with respect to the reference phase space probability measure $\pi$ if

[TABLE]

This condition is sufficient and necessary to define the likelihood ratio between probability measures

[TABLE]

The function $g(\mathbf{x})$ is also called Radon-Nikodym derivative (Radon-Nikodym theorem Oksendal2003 ) and can be used to construct the phase-space probability measure of the perturbed system, from the phase-space probability of the reference system:

[TABLE]

As consequence, if $g(\mathbf{x})$ can be calculated, one can estimate a phase space ensemble average (eq. 6) for the perturbed dynamics (eq. 27) from a set of states $S_{n}=\{\mathbf{x}_{1},...\mathbf{x}_{n}\}$ which has been generated by the reference dynamics (eq. 1)

[TABLE]

The notion of absolute continuity is valid also for path probability densities:

[TABLE]

Thus we can reweight a path ensemble average, by using the likelihood ratio between the path probability density $\mu_{\widetilde{P}}(\omega)$ and $\mu_{P}(\omega)$. For diffusion processes like (1) and (27), the likelihood ratio is given as

[TABLE]

where $\eta^{i}_{k}$ are the random numbers, along dimension $i$ at time step $k$, generated to integrate equation 1 of the reference dynamics and $\nabla_{i}U(\mathbf{r}_{k})$ is the gradient of the perturbation along dimension $i$ measured a the position $\mathbf{r}_{k}$. Note that to evaluate eq. 36 one needs the positions and the random numbers for every time step of the time-discretized trajectory. Eq. 36 is derived in appendix C. We remark that the quantity $M_{\tau,\mathbf{x}}(\omega)$ exists also for continuous paths ($\Delta t\rightarrow 0$). In this case the existence of the Radon-Nikodym derivative is guaranteed by the Girsanov theoremGirsanov1960 ; Oksendal2003 that states the conditions under which a perturbed path probability density $\widetilde{P}$ can be defined with respect to a reference path probability density $P$. The differences between time-continuous and time-discrete paths are discussed in the appendix A.

Analogous to the reweighting of phase-space ensemble averages (eq. 34), we can use $M_{\tau,\mathbf{x}}(\omega)$ to reweight path ensemble averages

[TABLE]

The second equality shows how to estimate the path ensemble average $\mathbb{E}_{\widetilde{P}}[f(\omega)]$ at the perturbed dynamics (eq. 27) from a set of paths $S_{\tau,\mathbf{x},m}=\{\omega_{1},\omega_{2},...\omega_{m}\}$ which has been generated by the reference dynamics (eq. 1).

II.4 Reweighting MSMs

When reweighting a MSM, we estimate the cross-correlation function $\widetilde{C}_{ij}(\tau)$ for the perturbed dynamics (eq. 27), from a set of paths of length $\tau$, $S_{\tau,m}=\{\nu_{1},\nu_{2}...\nu_{m}\}$, which has been generated by the reference dynamics (eq. 1). The cross-correlation function $\widetilde{C}_{ij}(\tau)$ is a combined phase space and path ensemble average (eq. 24). Thus, both averages have to be reweighted with the appropriated probability ratio

[TABLE]

With $a(\mathbf{x}_{0})=\mathbf{1}_{B_{i}}(\mathbf{x}_{0})$ and $f(\omega)=\mathbf{1}_{B_{j}}(\mathbf{x}_{n})$, we obtain

[TABLE]

which can be estimated from a set of paths $S_{\tau,m}=\{\nu_{1},\nu_{2}...\nu_{m}\}$ as

[TABLE]

where $[\mathbf{x}_{i}]_{k}$ is the $i$th time step of the $k$th path. As in eq. 16, the initial states of the paths are not fixed but are distrubuted according to the equilibrium distribution $\mu_{\pi}(\mathbf{x})$ of the unperturbed dynamics. Finally, the transition probability between set $B_{i}$ and set $B_{j}$ for the perturbed dynamics is obtained as

[TABLE]

The Radon-Nikodym derivative for path ensembles $g(\mathbf{x})$ contains the ratio of the partition functions $\nicefrac{{Z}}{{\widetilde{Z}}}$ as a multiplicative factor. Since this factor appears both in the numerator and the denominator of eq. 43, it cancels, and the partition functions do not have to be calculated.

$T_{ij}(\tau)$ (and analogously $\widetilde{T}_{ij}(\tau)$) is an element of the $s\times s$ MSM transition matrix $\mathbf{T}(\tau)$, where $s$ is the number of disjoint sets (microstates). We characterize the MSM by plotting and analyzing the dominant left and right eigenvectors of the transitiom matrix

[TABLE]

where $\mathbf{l}_{i}^{\top}$ denotes the transpose of vector $\mathbf{l}_{i}$. We assess the approximation quality of the MSM by checking wether the implied timescales

[TABLE]

are constant Swope2004 ; Prinz2011 .

II.5 Projection

We now consider a perturbation $U(\cdot)$ that does not directly affect the relevant coordinates used to construct the MSM. In such a situation the perturbation acts mainly on the coordinates directly perturbed and has a minor effect on the other degrees of freedom, in particular on the relevant coordinates that do not capture the full effect of the perturbation. Thus the reweighting may become problematic because the reweighting formula (36) is dominated by large, fluctuating gradients, therefore a rescaling of the force needs to be introduced. To address this issue, we propose to project the gradient of the perturbation onto the coordinates used to construct the MSM. Let’s assume that the MSM has been built on a combination of $d$ coordinates $\chi_{1},...\chi_{d}$, then the equation 36 is rewritten as:

[TABLE]

with

[TABLE]

III Methods

III.1 Two dimensional system

The Brownian dynamics on a two-dimensional potential energy function $V(x,y)$

[TABLE]

have been solved using the Euler-Maruyama scheme(b9, ) with an integration time step of $\Delta t=0.001$. The term $B_{t}^{i}$ denotes a standard Brownian motion in direction $i=x,y$, $\sigma=1$ is the volatility, and the random variables $\eta^{i}$ were drawn from a standard Gaussian distribution. The reference potential energy function was

[TABLE]

and the perturbed potential energy function was $\widetilde{V}(x,y,)=V(x,y)+U(y)$ with

[TABLE]

For both potential energy functions, trajectories of $8\times 10^{7}$ time-steps were produced. In both simulations, the path probability ratio $M_{\tau,\mathbf{x}}$ was calculated using eq. 36. The MSM have been constructed by discretizing each dimension $x$ and $y$ into 40 bins, yielding 1600 microstates. The chosen lagtime was $\tau=400$ time-steps.

III.2 Many-body system in three-dimensional space

We designed a six-particle system, in which five particles form a chain while the sixth particle branches the chain at the central atom (Fig. 5.A). The position of the $i$th particle is $\mathbf{r}_{i}\in\mathbb{R}^{3}$. The potential energy between two directly bonded atoms (blue lines in Fig. 5.A) was

[TABLE]

with $\mathbf{r}_{ij}=\mathbf{r}_{j}-\mathbf{r}_{i}$. The bond potential energy function is a tilted double well potential. This ensures that the potential energy function of the complete system has multiple minima with varying depths. No non-bonded interactions were applied. Thus, the reference potential energy function of the complete system was

[TABLE]

The Langevin dynamics (eq. 1) of this system have been solved using the BBK integrator (b5, ) with an integration time step of $\Delta t=0.001$. The masses $M$ of the particles, temperature $T$, the friction coefficient $\gamma$, and the Boltzmann constant $k_{B}$ were all set to one. The perturbed potential energy function was $\widetilde{V}(\mathbf{r})=V(\mathbf{r})+U(\mathbf{r})$ with

[TABLE]

The perturbation is a harmonic potential energy along the through-space distance between atoms $(2,4)$ and $(3,4)$, respectively (green dashed line in Fig. 5.A). For both potential energy functions, trajectories of $3.2\times 10^{8}$ time-steps were produced. The MSMs were constructed on the two-dimensional space spanned by the bond-vectors $\mathbf{r}_{12}$ and $\mathbf{r}_{13}$, i.e. on two coordinates which were not directly perturbed. Each dimension $x$ and $y$ has been discretized into 40 bins, yielding 1600 microstates. The chosen lagtime was $\tau=400$ time-steps. In both simulations, the path probability ratio $M_{\tau,\mathbf{x}}$ was calculated by projecting the gradient vectors $\nabla U\left(\mathbf{r}_{24}\right)$ and $\nabla U\left(\mathbf{r}_{34}\right)$ on the vectors $\mathbf{r}_{12}$ and $\mathbf{r}_{13}$ and subsequently evaluating eq. 47.

III.3 Alanine and Valine dipeptide

We performed all-atom MD simulations of acetyl-alanine-methylamide (Ac-A-NHMe, alanine dipeptide) in implicit and explicit water and of acetyl-valine-methylamide (Ac-V-NHMe, valine dipeptide) in implicit water. All simulations were carried out with the OPENMM 7.01 simulation package(b7, ), in an NVT ensemble at 300 K. Each system was simulated with the force field AMBER ff-14sb (b6, ). The water model was chosen according to the simulation, i.e. the GBSA model (Onufriev2004, ) for implicit solvent simulation and the TIP3P model (b8, ) for explicit solvent simulation. For each of these setups, the aggregated simulation time was 1 $\mu$s and we printed out the positions every nstxout=100 time steps, corresponding to 0.2 ps. A Langevin thermostat has been applied to control the temperature and a Langevin leapfrog integratorIzaguirre2010 has been used to integrate eq. 1. For implicit solvent simulations, interactions beyond 1 nm are truncated. For explicit solvent simulations, periodic boundary conditions are used with the Particle-Mesh Ewald (PME) algorithm to estimate Coulomb interactions.

In the alanine dipeptide simulations, we have perturbed the potential energy function of the backbone dihedral angles $\phi$ and $\psi$. The reference potential energy functions were

[TABLE]

[TABLE]

where the parameters have been extracted from force field files and $\pi$ denotes the mathematical constant. The perturbed potential energy function was a harmonic potential along each dihedral angle degree of freedom

[TABLE]

where $\kappa_{\phi}$ and $\kappa_{\psi}$ are the force constants, which could be adjusted after the simulation (Fig. 6.A and 6.B). The gradient $\nabla_{\mathbf{r}_{i}}U(\cdot)$ in eq. 36 is defined with respect to the Cartesian coordinates. Thus, applying the chain rule, the path probability ratio for the perturbation of the backbone dihedral angles is given as

[TABLE]

In the valine dipeptide simulation, we have perturbed the $\chi_{1}$ side-chain dihedral angle. The reference potential energy function was

[TABLE]

where the parameters have been extracted from force field files and $\pi$ denotes the mathematical constant. The perturbed potential energy function was a harmonic potential

[TABLE]

where $\kappa_{\chi}$ is the force constant, which could be adjusted after the simulation (Fig. 8.A). Thus, the path probability ratio is given as

[TABLE]

The MSM, for both alanine dipeptide and valine dipeptide simulations, have been constructed by discretizing the dihedral angles $\phi$ and $\psi$ into 36 bins each, yielding 1296 microstates. The chosen lagtime was 20 ps for both the implicit and explicit solvent simulation.

IV Results and discussion

IV.1 Efficient implementation

To estimate a MSM at a perturbed potential energy function $V(\mathbf{x})+U(\mathbf{x},\kappa)$ with the dynamical reweighting method, one simulates a long trajectory $\mathbf{x}(t)$ at a reference potential energy function $V(\mathbf{x})$ using an integration time step $\Delta t$. From this trajectory $m$ short paths of length $\tau=n\Delta t$ are extracted, yielding a set of paths $S_{\tau,m}=\{\nu_{1},\nu_{2},...\nu_{m}\}$, which can subsequently be used to evaluate eq. 42, where $g([\mathbf{x}_{0}]_{k})$ is given by eq. 31 and $M_{\mathbf{x},\tau}(\nu_{k})$ is given by eq. 36 or eq. 47. As discussed in section II.4, the factor $\nicefrac{{Z}}{{\widetilde{Z}}}$ cancels for the estimate of $\widetilde{T}_{ij}(\tau)$ (eq. 43) and the partition functions do not need to be calculated.

Let us assume that the simulation integrates the Langevin equation of motion for the reference potential energy function $V(\mathbf{x})$ (eq. 1). To estimate a MSM with transition probabilities $\widetilde{T}_{ij}(\tau)$ for the dynamics in the perturbed potential energy function from a set of paths $S_{\tau,m}=\{\nu_{1},\nu_{2},...\nu_{m}\}$, we need to know the value of the perturbation $U([\mathbf{x}_{t}]_{k},\kappa)$ at every time step $t$ and for each path $\nu_{k}$, the gradient of the perturbation $\nabla U([\mathbf{x}_{t}]_{k},\kappa)$ at every time step $t$ and for each path $\nu_{k}$, the random numbers $\eta_{k}$ generated at every time step of the simulation of each path $\nu_{k}$, and the volatility $\sigma$. The volatility is determined by the temperature and the friction coefficient (eq. 3), both of which are input parameters for the simulation algorithm. To calculate the other three properties one needs to know the positions and the random numbers at every integration time step of the simulation.

In a naive implementation of the reweighting method one would hence write out the positions and random number at every time step and calculate $g([\mathbf{x}_{0}]_{k})$ and $M_{\mathbf{x},\tau}(\nu_{k})$ in a post-analysis step. The advantage of this approach is that the set of paths can be reweighted to the path probability measure of any perturbation $U(\mathbf{x},\kappa)$, as long as the absolute continuity is respected (eq. 30). On the other hand, this approach is hardly practical because writing out a trajectory at every integration time step quickly fills up any hard disc and slows down the simulation considerably.

We therefore decided to compute the probability ratios $g([\mathbf{x}_{0}]_{k})$ and $M_{\mathbf{x},\tau}(\nu_{k})$ “on the fly” during the simulation. In practice, the lag time $\tau$ of a MSM can only assume values integer multiples of the frequency nstxout at which the positions are written to file, i.e. of $\tau=n\Delta t=A\cdot\mbox{{nstxout}}\cdot\Delta t$ with $A\in\mathbb{N}$. The discretized It$\hat{\mathrm{o}}$ integral and the discretized Riemann integral in eq. 36 are sums from time step $k=0$ to $k=n$, which can be broken down into $A$ sums of size nstxout

[TABLE]

Thus, we calculate the terms

[TABLE]

and

[TABLE]

“on the fly” and write out the results at the same frequency nstxout as the positions. The path probability ratio is reconstructed after the simulation as

[TABLE]

The potential energy of the perturbation $U(\mathbf{x}_{t})$ is written out at the frequency nstxout and the complete weight $g([\mathbf{x}_{0}]_{k})\cdot M_{\mathbf{x},\tau}(\nu_{k})$ is calculated during the construction of the MSM. The lag-time $\tau$ can be chosen and varied after the simulation. The modification of the MD integrator can be readily implemented within the MD software package OpenMM (b7, ). An example script is provided in the supplementary material.

The approach requires that the perturbation potential energy function $U(\mathbf{r},\kappa)$ is chosen prior to the simulation. Note however that, if the perturbation potential energy function depends linearly on the parameter $\kappa$, i.e. if $\kappa$ is a force constant $U(\mathbf{r},\kappa)=\kappa\cdot U(\mathbf{r})$, then so do the two integrals $I(a,U(\mathbf{r},\kappa))=\kappa\cdot I(a,U(\mathbf{r}))$ and $R(a,U(\mathbf{r},\kappa))=\kappa\cdot R(a,U(\mathbf{r}))$. Thus, it is sufficient to calculate $I(a,U(\mathbf{r}))$ and $R(a,U(\mathbf{r}))$ “on the fly” and to scale the integrals after the simulation to any desired value of $\kappa$. A single simulation is sufficient to allow for reweighting to a whole series of perturbation potential energy functions. Also, the integrator can be modified such that the It$\hat{\mathrm{o}}$ integrals $I_{1}(a),I_{2}(2)...$ and the Riemann integral $R_{1}(a),R_{2}(2)...$ of several functionally different perturbation energy functions $U_{1}(\mathbf{x}),U_{2}(\mathbf{x})...$ are calculated. Thus, using a single reference simulation one can reweight to several functionally different perturbations and scale these perturbations by an arbitrary force constant.

The computional cost of adding the calculation of $I_{i}(a)$ and $R_{i}(a)$ for a perturbation energy function $U_{i}(\mathbf{x})$ is modest. The blue line in Fig. 3 shows the computational costs for simulating alanine dipeptide in implicit water with $n=0,1,2,3,4,5$ different perturbations as the number of days required to obtain a trajectory of 1 $\mu$s for each perturbation on a small workstation. The red line in Fig. 3 shows the computational cost of implementing the same $n$ perturbations into a single reference simulation. For a single perturbation 30% of the computational cost is saved by dynamical reweighting, whereas for 5 perturbations more than half of the computational cost is saved. The gain is even greater, if the dependence on a force constant is to be studied. Moreover, in simulation boxes of larger systems with explicit solvent, the cost of evaluating of the perturbation potential energy is small compared to the normal force evaluations, such that for these systems the computional cost of calculating the $I_{i}(a)$ and $R_{i}(a)$ “on the fly” becomes negligible. We remark that if $I_{i}(a)$ and $R_{i}(a)$ are too large, the Girsanov reweighting method might become numerically intractable.

IV.2 Stochastic forces

The dynamical reweighting method has been derived by using the Langevin equation of motion as starting point. The concept of the path probability is inextricably linked to the presence of Gaussian random forces in the equation of motion (see appendix A). However from a physical perspective, the Langevin dynamics only approximates the Hamiltonian dynamics of the complete system, by splitting the system in a subsystem $S$ and a heat bath $B$ and replacing the interaction of the subsystem with the heat bath by a friction term and a stochastic process. In most MD simulations, thermostatted versions of the Hamiltonian dynamics are propagated in which no stochastic forces are generated. In principle, one could split such a simulation box into two subsystem $S$ and $B$ and use the forces which subsystem $B$ exerts on subsystem $S$ as substitute for the random forces. Unfortunately, even if the heat bath $B$ consists of Lennard-Jones particles, the forces on $S$ deviate too much from a Gaussian random process to be used in the dynamical reweighting method (data not shown). We therefore decided to use the Langevin leapfrog integratorIzaguirre2010 to integrate the Langevin equation of motion for both the implicit and explicit solvent simulations and used the random forces generated by the integrator to reweight the path ensemble.

IV.3 Two dimensional system

As a first application, we consider the Brownian dynamics of a particle moving on a two-dimensional potential energy function (eq. 49). The reference potential energy function $V(x,y)$ (eq. 50, Fig. 4.A) has two minima at $(-1,-1)$ and $(1,1)$ which are connected by a transition state at $(0,0)$. We added at a perturbation (eq. 51) which tilts the energy function $\widetilde{V}(x,y)$ along the direction $y$, such that the mininum at $(1,1)$ becomes much deeper than the minimum at $(-1,-1)$. Fig. 4.B and 4.E show the dominant left MSM eigenvectors of the two systems (direct MSMs). The first eigenvector corresponds to the equilibrium distribution. In both cases, the second eigenvector represents the transition between the two wells, while the third eigenvector corresponds to an exchange of probability density between the transition state region and the two wells.

Fig. 4.C shows the MSM eigenvectors obtained by reweighting the simulation in $V(x,y)$ to the perturbed potential energy function $V(x,y)+U(y)$ (reweighted MSM). Both eigenvectors (eq. 45) and implied timescales (eq. 46) are in perfect agreement with Fig. 4.E. This confirms that reweighting MSMs using the Girsanov formula works well for low-dimensional Brownian dynamicsSchuette2015b . Note that the simulation in the reference potential energy function $V(x,y)$ exhibits frequent transitions between the two minima, and thus eq. 35 is certainly fulfilled. By contrast in the perturbed system $\widetilde{V}(x,y)$, the simulation sampled considerably fewer transitions between the minima and only a small fraction of the simulation time was spent in the minimum in the lower left corner. Fig. 4.F shows the dominant eigenvectors of a MSM of the reference dynamics constructed by reweighting the perturbed path ensemble. The fact that Fig. 4.C is in excellent agreement with Fig. 4.B demonstrates that the reweighting method yields accurate results even for path ensembes with low numbers of transitions across the largest barriers of the system and thus suggests that it can be applied to high-dimensional dynamics on rugged potential energy surfaces.

IV.4 Many-body system in three-dimensional space

As a second example, we studied a six-particle system, in which five particles form a chain while the sixth particle branches the chain at the central atom Fig. 5.A . We constructed the MSM on the central bonds $\mathbf{r}_{12}$ and $\mathbf{r}_{13}$, but applied the perturbation potential energy function along the through-space distances $\mathbf{r}_{24}$ and $\mathbf{r}_{34}$. Thus, the perturbation was projected onto the reaction coordinates $\mathbf{r}_{12}$ and $\mathbf{r}_{13}$ during the reweighting. The reference potential energy function for each bond is a tilted double well potential, such that the reference system exhibits four metastable states (first row of Fig. 5.B). The six-particle system and the reference potential energy function is symmetric. Thus the MSM of the reference system has two degenerate dominant eigenfunction with implied timescales of $1.1\cdot 10^{4}\Delta t$. with a relative error due to sampling of 1.6 %.

The perturbation contracts the bonds, thereby stabilizing the metastable state at the lower left corner in the 1st eigenvector (second row of Fig. 5.B). Its effect is to break the symmetry and to accelerate the dynamics, yielding implied timescales of $6.7\cdot 10^{4}\Delta t$ and $5.9\cdot 10^{4}\Delta t$. The third row of Fig. 5.B shows the eigenvectors and the implied time scales obtained by reweighting the simulation at the reference potential energy function to the perturbed potential energy function. The relative error of the implied timescale of the second and third eigenvector is 5.2% and 0.7%, respectively.

IV.5 Alanine dipeptide and valine dipeptide

Figure  6 shows the results for alanine dipeptide (Ac-A-NHMe) in implicit water. The MSM has been constructed on the $\phi$ and $\psi$ backbone dihedral angles and the slow eigenvectors of the unperturbed system are shown in Fig. 6.C. The first eigenvector shows the typical equilibrium distribution in the Ramachandran plane Vitalini2015 . The second eigenvector represents a torsion around the $\phi$ angle and corresponds to a kinetic exchange between the Lα-minimum ($\phi>0$) and the $\alpha$-helix and $\beta$-sheet minima ($\phi<0$). The associated timescale is 2.8 ns. The green arrows in Fig. 6.C represent the frequency of the transitions, with the transition L${}_{\alpha}\leftrightarrow\beta$-sheet conformation occuring more frequently than the transition L${}_{\alpha}\leftrightarrow\alpha$-helical conformation. The third eigenvector represents a transition $\beta$-sheet $\longleftrightarrow\alpha$-helical conformation, i.e. a torsion around $\psi$, and is associated to a timescale of 27 ps.

We perturbed the dynamics by adding a harmonic potential to the dihedral angle potentials of the $\phi$- and $\psi$-angle (Fig. 6.A and 6.B, eq. 57 with $\kappa_{\phi}=0.5$ and $\kappa_{\psi}=0.5$). The $\alpha$-helical region is somewhat stabilized by the perturbation but otherwise the dominant eigenvectors are very similar to the unperturbed system (second row in Fig. 6.C). However, the perturbation changes the relative frequency of the two possible transitions (green arrows) in the second eigenvector, resulting in an increased implied timescale of 4.5 ns. The third row of Fig. 6.C shows the dominant eigenvector of the perturbed system obtained by reweighting the reference simulations. The results are in excellent agreement with the direct simulation of the perturbed systems. The relative error of the implied timescale associated to the second eigenvector is 4.1%.

Fig. 7.A shows the implied timescale test (eq. 46) for alanine dipeptide in implicit water. All three systems (reference, perturbed and reweighted) show constant implied timescales, indicating that the MSM are well converged. Moreover, the graph shows that the reweighting method can recover the implied timescales of the perturbed system over a large range of lag times $\tau$. We have repeated the alanine dipeptide simulations in explicit water. The eigenvectors are very similar to those of alanine dipeptide in implicit water (data not shown), but the associated implied timescales differ from the implicit solvent simulations (Fig. 7.B). In the unperturbed system, the implied timescale of the second eigenvector was 3.1 ns and the implied timescale of the third eigenvector was 75 ps. Perturbing the dihedral angle potentials, slightly increased the implied timescale of the second eigenvector to 3.5 ns and left the implied timescale of the third eigenvector unaffected. The implied timescales of the reference simulation and the direct simulation of the perturbed system were constant, whereas we noticed a slight drift in the implied timescale of the second eigenvector for the reweighted MSM. At $\tau=45$ ps, the measured implied timescale is 3.4 ns which corresponds to a relative error of 2.8 %.

We also tested the dynamical reweighting method on valine dipeptide (Fig. 8). Here, however, we perturbed the potential energy function of the $\chi$-side chain dihedral angle (Fig. 8.A) by adding a harmonic potential, while constructing the MSM on the $\phi$ and $\psi$ backbone dihedral angles. Thus, the perturbation did not directly act on the variables of the MSM. The perturbation of the $\chi$ angle had no affect on the eigenvectors of the MSM (Fig. 8.B), but it did change the implied time scales. It caused a decrease of the implied time scale of the second eigenvector from 1.3 ns in the reference simulation to 1.0 ns in the perturbed simulation and a slight increase of the implied time scale of the third eigenvector from 159 ps in the reference simulation to 170 ps in the perturbed simulation. Reweighting the reference simulation to the perturbed potential energy function recovered the results of the direct simulation of the perturbed system. The implied timescales obtained by the reweighting calculation were 1.0 ns (relative error: 4.9% before rounding to ns) for the second eigenvector and 179 ps (relativ error: 5.3%) for the third eigenvector. This shows that the dynamical reweighting method also works, when the perturbation acts on degrees of freedom which are part of the relevant coordinates on which the MSM is constructed.

Fig. 9 illustrates how to use the dynamical reweighting method to study the influence of a force constant on the molecular dynamics. It shows the implied timescale of the second and third MSM eigenvectors of alanine dipeptide in explicit water as a function of the force constant $\kappa_{\phi}$, where the perturbation potential energy is given by eqs. 55 and 57. The scan has been repeated with different values of the force constant $\kappa_{\psi}$ for the potential energy function of the $\psi$-backbone dihedral angle (eqs. 56 and 57). The dynamics is more sensitive to a change in the value of $\kappa_{\phi}$ than to a change of $\kappa_{\psi}$, but the overall effect of the perturbation is moderate. It is important to point out that MSMs which are summarized in Fig. 9 have been constructed from a single simulation at the reference potential energy function. During this simulation, the It$\hat{\mbox{o}}$ integral $I(a)$ and the Riemann integral $R(a)$ have been calculated for $\kappa_{\phi}=1$ and $\kappa_{\psi}=1$, and the force constants have been scaled after the simulation during the construction of the MSM, as described in section IV.1.

V Conclusion

We have presented the Girsanov reweighting scheme, which is a method to study the dynamics of a molecular system subject to an (external) perturbation $U(\kappa,\mathbf{x})$ of the (reference) potential energy function $V(\mathbf{x})$. It allows for the estimation of a dynamical model, e.g. a MSM, of the perturbed system from a simulation at the reference potential energy function. The underlying assumption is that the equation of motion generates a path ensemble and that we can define a probability measure on this ensemble. A perturbation of the potential energy function causes a modification of the probability measure. The Girsanov theorem guarantees that the probability ratio between these two measures exists (under certain conditions) and leads to an analytical expression for this ratio (eq. 36). By reformulating the MSM transition probabilities as path ensemble averages, we can apply the Girsanov reweighting scheme to obtain the transition probabilities of the perturbed system from a set of paths generated at the reference potential energy function. The method can be extended to the variational approaches Nuske2014 ; Vitalini2015b , milestoning approaches Schuette2011 ; Lemke2016 or tensor approaches Nuske2016 to molecular dynamics, since in each of these methods the molecular transfer operator is discretized and the resulting matrix elements are estimated as path ensemble averages.

Calculating the path probability ratio requires the knowledge of the random forces at each integration time step. For the explicit solvent simulations, we have introduced stochastic forces by using a Langevin thermostat. In an efficient implementation of the method, two terms which are needed to calculate the probability ratio should be calculated “on the fly” during the simulation. We have demonstrated this using the MD simulation toolkit OPENMM b7 .

Two other dynamical reweighting schemes for MSMs have been published in recent years. In the reweighting scheme for parallel tempering simulations Prinz2011c ; Chodera2011 , the path probability density is defined for time discretized paths at a reference temperature and then reweighted to different temperatures using statistically optimal estimators Shirts2008 ; Minh2009 . Like the Girsanov reweighting scheme, this method relies on the random forces at each integration time step and reweights the contribution of each path to the estimate of the transition probability individually. We remark that the reweighting to different thermodynamic states cannot be extended to the limiting case of continuous paths because for different volatilities the path probability ratio relative to the Wiener process cannot be defined (see appendix A). This however seems to be of little practical importance. In the TRAM method Mey2014 ; Wu2014 ; Wu2016 , rather than reweighting the probabilities of each individual path, the MSM transition probabilities $T_{ij}(\tau)$ are directly reweighted to a different potential energy function or a different thermodynamic state using a maximum likelihood estimator for the transition counts. This becomes possible, if one additionally assumes that the dynamics is in local equilibrium within each microstate of the MSM, which possibly renders the method more sensitive to the MSM discretization than path based reweighting methods.

We have tested the Girsanov reweighting method on several systems, ranging from diffusion in a two-dimensional potential energy surface to alanine dipeptide and valine dipetide in implicit and explict water. Importantly, the direct simulations of the perturbed potential energy function (Figs. 5, 6, 8) are only included as a validation for the method. In an actual application one would only simulate the system at the reference potential energy function and then reweight to the perturbed potential energy function, thus saving the computational time of the direct simulation of the perturbed system.

Girsanov reweighting could be useful in several areas of research. First, one can very efficiently test the influence of a change in the potential energy function on the dynamics of the molecule. The influence of a change in the force constant on the dynamics is particularly easy to study. The Girsanov reweighting method can therefore be applied to improve the dynamical properties of force fields Vitalini2015 , by for example tuning the force constants to match an experimentally measured correlation time. Similarly, one can use Girsanov reweighting to understand the influence of restraining potentials Cesari2016 ; Keller2007 on the dynamics of the system. Second, the method can be used to understand which degrees of freedom have the largest influence on the slow modes of the molecule Tsourtis2015 ; Arampatzisa2016 . For example, for alanine dipeptide, we showed that the slow dynamic modes are more sensitive to a force field variation in the $\phi$-backbone dihedral angle than they are to a variation in the $\psi$-backbone dihedral angle. Last but not least, Girsanov reweighting can be used to account for the effect of any external potential which has been added to the simulation in order to enhance the sampling. Thus, one can for example estimate MSMs from metadynamics simulationsHuber1994 ; Laio2002 , Hamilton replica exchange simulations Sugita2000 or umbrella sampling simulations Torrie1977 .

VI Supplementary material

Example script for the implementation of the dynamical reweigthing method with OpenMM b7 .

Acknowledgements.

This research has been funded by Deutsche Forschungsgemeinschaft (DFG) through grant CRC 1114 “Scaling Cascades in Complex Systems”, Project B05 “Origin of the scaling cascades in protein dynamics”.

Appendix A Discussion over path probability measure of time-discrete and continuos paths

In the theory section, we assumed to work with time-discrete trajectories and we have defined the corresponding concepts of path ensemble and path probability measures. In this section, we show why we cannot extend these concepts to continuous paths.

Let’s consider a diffusion process $x_{t}\in\Gamma\subset\mathbb{R}$ solution of the stochastic differential equation

[TABLE]

where $a(\cdot)$ is a drift, $B_{t}$ is a Brownian motion, $\sigma$ is a volatility and $\tau$ is the total time. The equation can be discretized in time according to several schemes. We consider the Euler-Maruyana method:

[TABLE]

where $a_{k}=a(x_{t})$, $\Delta t$ is a time-step and $\eta$ is a random number drawn from a standard Gaussian distribution.

We now recall the same definitions given in the theory section. A time-discrete path is the set of the points $x_{k}$ generated by the equation 71 $\omega=\{x(t=0)=x_{0},x_{1},x_{2},...,x_{n}\}$ where $n\in\mathbb{N}$ is the total number of time steps and $\tau=n\Delta t$.

If we consider an ensemble of paths starting from the same point $x_{0}$ and length $\tau$, we can define a path space $\Omega_{\tau,x}=\Gamma^{n}$, a path probability measure $P$ and the associated path probability density $\mu_{P}(\omega)$ as defined respectively in eq. (7) and eq. (9)

The path probability density function is defined as in the theory section, but here we remark that, in a rigorous formalism, a density function is a derivative respect to the Lebesgue measure $\mathrm{d}\omega=\mathrm{d}x_{1}\mathrm{d}x_{2}...\mathrm{d}x_{n}$, i.e. $\mu_{P}(\omega)=\frac{\mathrm{d}P}{\mathrm{d}\omega}$.

For Brownian motion with drift, the transition probability density $p(x_{k-1},x_{k};\Delta t)$ that appear in eq. (7) and eq. (9) is well defined for time-discrete paths:

[TABLE]

Let’s see now what happen if we take $\Delta t\rightarrow 0$, i.e. if we consider continuous paths:

The normalization constant tends to infinity:

[TABLE] 2. 2.

The quantity in the exponential function becomes a derivative of the process $x_{t}$:

[TABLE]

Then, formally

[TABLE]

But the derivative $\frac{\mathrm{d}x_{t}}{\mathrm{d}t}$ is not defined, because even though Brownian motion is everywhere continuous, it is nowhere differentiable. 3. 3.

The continuous path is made by an infinite number of points. If $n\rightarrow+\infty$, then the infinitesimal volume space $\mathrm{d}\omega=\mathrm{d}x_{1}\mathrm{d}x_{2}...\mathrm{d}x_{\infty}$ becomes infinite-dimensional. However, the Lebesgue measure of an infinite-dimensional space does not exist.

In conclusion, we cannot define a path probability density for a continuous path respect to the Lebesgue measure.

On the other hand we can define the probability density of a path $\omega$ respect to a Wiener measure $W$. By computing explicitly the ratio between the probability density associated to a (discrete) brownian motion with drift and a probability density associated to a Brownian motion without drift, but with the same volatility, we get the following expression:

[TABLE]

We observe that if we take the limit $\Delta t\rightarrow 0$, the three problems described above are solved:

The normalization constant cancels. 2. 2.

The term $x_{k}-x_{k-1}$ is not divided by the time-step $\Delta t$, so there is not derivative respect to time. 3. 3.

The ratio between the two probability densities cancels also the infinitesimal volume $\mathrm{d}\omega$. The new probability density is defined respect to the Wiener measure, not respect the Lebesgue measure.

Thus in the continuous case the likelihood (76) is written as:

[TABLE]

This quantity denotes the probability density, in the continuos case, of a path $x_{t}$ solution of a diffusion process, respect to the Wiener measure. The complete computation of the formulas (76) and (77) is reported in the appendix C. The Girsanov theorem states the condition under which the likelihood (77) exists.

Appendix B The Girsanov theorem

The Girsanov theorem (Girsanov1960, ; Oksendal2003, ) is an important result of stochastic analysis that states the conditions under which a change of measure is possible for diffusion processes.

Let’s consider two diffusion processes $x_{t}\in\mathbb{R}$ and $y_{t}\in\mathbb{R}$ defined on the probability space $(\Omega,\mathcal{F},\{\mathcal{F}_{t}\}_{t\geq 0})$, respectively with probability measures $P$ and $Q$ on the filtration $\{\mathcal{F}_{t}\}_{t\geq 0})$, be solutions of the stochastic differential equations:

[TABLE]

where the functions $a:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ and $b:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ are measurable, $\sigma$ is the volatility, $B_{t}$ is an $n$-dimensional Brownian motion, $\tau\leq\infty$ and $x_{0}=y_{0}\in\mathbb{R}^{n}$.

Suppose that there exists a process $\xi_{t}$ such that

[TABLE]

and define

[TABLE]

If the Novikov condition

[TABLE]

is satisfied, then the process $M_{\tau}$ is a martingale [(Oksendal2003, ),pag.31] and we can define the probability measure $Q$ on $\mathcal{F}_{t}$ by the Radon-Nikodym derivative

[TABLE]

and the process

[TABLE]

is a new Brownian motion w.r.t. $Q$. As consequence, for any measurable functional $f(x)$,

[TABLE]

Remark

The measures $P$ and $Q$ on the filtration $\mathcal{F}_{t}$ represent the probability measures on a path space. In other words, given a set of trajectories, $P$ and $Q$ return the probability of a certain path. The concept of absolute continuity means that if a path is not allowed under $P$, then it is not allowed under $Q$:

[TABLE]

for every measurable set $F_{1},...,F_{n}\subseteq\mathcal{F}_{\tau}$ with $\tau>0$ fixed. Because $M_{\tau}$ is strictly positive, the probability measure $P$ is also absolutely continuous w.r.t $Q$:

[TABLE]

then the two measures $P$ and $Q$ are equivalent (or mutually absolutely continuous w.r.t each other) and we can write:

[TABLE]

Appendix C Derivation of the Girsanov formula

Let’s now consider two diffusion processes with different drift terms, solutions of the stochastic differential equations:

[TABLE]

with initial conditions $x_{0}=x\in\mathbb{R}$. We discretize the equations 88, 89 with the Euler-Maruyama method:

[TABLE]

where $a_{k}=a(x_{k})$ and $b_{k}=b(x_{k})$ are two different drift terms. The transition probability densities are defined respectively as

[TABLE]

for the first equation 88 and

[TABLE]

for the second equation 89.

The path probability density $\mu_{P_{a}}(\omega)$ associated to the equation 88 and $\mu_{P_{b}}(\omega)$ associated to the equation 89 are defined as in eq. 9. Thus, the ratio is given as

[TABLE]

If we take the limit $\Delta t\rightarrow 0$, the two sums in the exponential function converge respectively to the Ito integral and the Riemann integral. The jump $x_{k}-x_{k-1}$ is Wiener jump with drift $a_{k}\Delta t$:

[TABLE]

Finally we have:

[TABLE]

Thus we get the more general Girsanov formula valid also for two diffusion processes with different drifts.

[TABLE]