Correcting Equilibrium Markov Models - American Chemical Society

Subscriber access provided by READING UNIV

Article 2

2

Caliber Corrected Markov Modeling (CM): Correcting Equilibrium Markov Models Purushottam D Dixit, and Ken Dill J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b01126 • Publication Date (Web): 11 Jan 2018 Downloaded from http://pubs.acs.org on January 11, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Caliber Corrected Markov Modeling (C2M2): Correcting Equilibrium Markov Models Purushottam D. Dixit∗,† and Ken A. Dill‡ Department of Systems Biology, Columbia University, New York, NY 10032, and Laufer Center for Quantitative Biology, Department of Chemistry, and Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY, 11790 E-mail: [email protected]

Abstract Rate processes are often modeled using Markov-State Models (MSM). Suppose you know a prior MSM, and then learn that your prediction of some particular observable rate is wrong. What is the best way to correct the whole MSM? For example, molecular dynamics simulations of protein folding may sample many microstates, possibly giving correct pathways through them, while also giving the wrong overall folding rate, when compared to experiment. Here, we describe Caliber Corrected Markov Modeling (C2 M2 ): an approach based on the principle of maximum entropy for updating a Markov model by imposing state- and trajectory-based constraints. We show that such corrections are equivalent to asserting positiondependent diffusion coefficients in continuous-time continuous-space Markov processes modeled by a Smoluchowski equation. We derive the functional form of the diffusion coefficient explicitly in terms of the trajectory-based constraints. We illustrate with examples of 2D particle diffusion and an overdamped harmonic oscillator. ∗

To whom correspondence should be addressed Department of Systems Biology, Columbia University, New York, NY 10032 ‡ Laufer Center for Quantitative Biology, Department of Chemistry, and Department of Physics and Astronomy, Stony Brook University, Stony Brook, NY, 11790 †

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 29

The problem: Correcting Markov Models from data Consider the following type of problem. You have a network of states a = 1, 2, . . .. You have a Markov model with known a priori transition probabilities Pab between all pairs of states. Now, you learn from data that some single global average rate quantity predicted by this model is incorrect. What is the ‘best’ way to correct the full transition matrix to bring it into consistency with the new limited information? This is a common problem. First, Markov models are ubiquitous. Among many other things, they are used to study folding, binding and mechanisms of action of biomolecules (1) , chemical and biochemical reaction networks (2,3) , and the evolutionary dynamics of organisms (4) . Second, such models often require many states, and yet are faced with limited experimental data, or limited physical insights that can constrain the model. Here’s an example. Computer simulations of proteins identify several different metastable conformational states. The simulated dynamics among these states can then be captured in Markov State Models (MSMs) (1) . But, the underlying forcefields are imperfect, so global rate quantities found from molecular dynamics simulations – such as their folding times, or rates along dominant reaction coordinates – are often found to be in error. Notably, estimation of kinetic rates may be inaccurate even with accurate forcefields owing to infrequent sampling of rare transitions (5) . Because the MSMs are likely to be mimicking the relative rates of microscopic processes, correcting the MSMs to agree with one or more experimental observables is likely to approximate well the full microscopic kinetics. Here, we describe a solution to this problem, ‘Caliber Corrected Markov Modeling’ (C2 M2 ), which employs the principle of Maximum Caliber (6) , the dynamical version of the Maximum Entropy principle of inference. Prior work solves a related problem of correcting an equilibrium distribution. Pitera and Chodera (7) developed an approach to fix equilibrium distributions based on experimental constraints. However, we need a different approach here, for two reasons: (1) There are infinitely many Markov models of the dynamics that are consistent with a given equilibrium distribution. And: (2) our focus is on the dynamics, because we are interested in biological mechanisms, not just the equilibrium states. For example, a drug’s efficacy is 2


Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


often determined by the dynamics of its unbinding from target proteins, not just its equilibrium binding strength (8,9) . Moreover, computational as well as experimental estimates of transition rates are much more error prone compared to equilibrium distributions. In the present work, we address the following general question: how do we update a detailed-balanced equilibrium Markov model so that it satisfies with user-imposed equilibrium and dynamical constraints? Relative entropy is commonly used to update models of biophysical systems (7,10,11) . Accordingly, we seek a Caliber Corrected Markov Model (C2 M2 ) that reproduces the imposed constraints and has the maximum relative path entropy (or a minimum Kullback-Leibler divergence) with respect to a prior Markov model. We first review the work of Pitera and Chodera (7) which serves as the first step in our theoretical development. Consider a system with discrete states {a, b, c, . . . } with Hamiltonian H at thermal equilibrium with its surroundings. The equilibrium distribution over states is given by xa ∝ e−βH(a) . Here β is the inverse tempearture. Conisder a state-dependent property f (a), for example, the end-to end distance of a peptide. The ensemble average hf ix is given by

hf ix =

X

xa f (a).

(1)

a

Imagine a situation where the model prediction hf ix does not agree with the corresponding experimentally measured ensemble average f¯. How do we then update the equilibrium distribution xa → ya (or equivalently the Hamiltonian H) such that the biased distribution {ya } reproduces f¯? Pitera and Chodera (7) seeked an updated equilibrium distribution {ya } that least deviated from the prior distribution {xa } while reproducing the ensemble average hf iy = f¯. They invoked the principle of maximum relative entropy (minimum Kullback-Leibler divergence). Briefly, one maximizes the relative entropy

S=−

X

ya log

a

3

ya xa


(2)


Page 4 of 29

subject to constraint X

ya f (a) = f¯

(3)

a

and

P

a

ya = 1. Carrying out the maximization using Lagrange multipliers (7) , ya ∝ xa e−λf (a) ∝ e−βH(a)−λf (a) .

(4)

The Lagrange multiplier λ dictates the deviation in the prediction of hf i between the unbiased ensemble {xa } and the biased ensemble {ya }. For example, if λ = 0, ya = xa . How do we impose similar biases in dynamics? Below, we develop our maximum entropy framework by updating a continuous time continuous space Smoluchowski equation for a particle diffusion on a one dimensional free energy landscape. Generalizations to continuous time discrete space and discrete time discrete space models are presented along the way.

Example of a particle diffusing with bias along one dimension Consider a particle diffusing in one dimension between a ∈ [−L, L] (see Fig. 1). If p(a, t) is the instantaneous probability distribution and xa ∝ exp (−βF (a)) is the equilibrium distribution, the dynamics of p(a, t) are described by the diffusion equation: ∂ ∂ p(a, t) = D0 ∂t ∂a

βF (a) −βF (a) ∂ e e p(a, t) ∂a

(5)

Figure 1: Diffusion of a particle in one dimension. Discretized scheme of a particle diffusing in one dimension between a = −L and a = L divided in 2n + 1 nodes labeled from −n to +n.

4


Page 5 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In Eq. 5, D0 is the diffusion coefficient and F (a) is the free energy surface. The diffusion coefficient sets the time scale of the system. In practice, it can be determined as the constant of proportionality between the ensemble average of the variance of the displacement a(t) − a(0) and the time t (12,13) . We can discretize the partial differential equation in steps of da and write a continuous-time discrete-state Markov process (13,14) , d p(a, t) = ωa−da,a p(a − da, t) + ωa+da,a p(a + da, t) dt − p(a, t) (ωa,a+da + ωa,a−da )

(6)

where the transition rates are given by (13,14) r ωa,a+da = D0

xa+da da2 . xa

(7)

Eq. 6 can be time-discretized using a small time interval dt. We write

p(a, t + dt) =

X

Pba × p(b, t)

(8)

b

where the transition probabilities are given by

Pab = ωab dt if b 6= a

(9)

and

Paa = 1 −

X

Pab

(10)

b6=a

Eq. 10 ensures that probabilities are conserved and normalized throughout the time evolution. We carry out the desired state- and trajectory-based biasing for the discrete time discrete state Markov model and then take appropriate continuous limits. First, we comment on the nature of trajectory-based observables. Consider a dynamical variable rab that is defined over individual transitions of a trajectory of the Markov process. An 5



Page 6 of 29

example of rab is the number of contacts formed/broken in a single time step by a polymer. The average rab over an ensemble of stationary state trajectories is given by (15–17)

hri =

X

xa Pab rab .

(11)

a,b

We want to modify the Markov model described by Eq. 9 such that the updated Markov model (with transition probabilities {kab }) has its equilibrium distribution equal to {ya } and reproduces the trajectory-ensemble average r¯ which is different than hri defined in Eq. 11. Previously, Wan et al. (18) have addressed the problem of updating Markov processes by updating their equilibrium distribution alone (see also Zhou et al. (19) ). Recently, Olsson et al. (20) also addressed the problem of fixing the equilibrium distribution within the Bayesian framework. We proceed by maximizing the relative entropy (18,21)

S=−

X

ya kab log

a,b

kab Pab

(12)

subject to constraint X

ya kab rab = r¯

(13)

a

and X

ya kab = ya ∀ a,

(14)

ya kab = yb ∀ b, and

(15)

ya kab = yb kba ∀ a and b.

(16)

b

X a

Eq. 14 ensures that state probabilities are conserved and normalized throughout the time evolution of the Markov process. Eq. 15 imposes {ya } as the stationary distribution of the Markov process. Finally, Eq. 16 explicitly imposes microscopic detailed balance. We write the Caliber C by incorporating these constraints using Lagrange multipli6


Page 7 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ers (6,15–18)

C = −

X

ya kab log kab +

a,b

X

ya kab log Pab

a,b

! +

X a

la

X

ya kab − ya

! +

X

b

mb

X

ya kab − yb

a

b

! +

X

ab (ya kab − yb kba ) + γ

a,b

X

ya kab rab − r¯ .

(17)

a,b

In Eq. 17, Lagrange multipliers {la } impose constraints in Eq. 14. {mb } impose constraints in Eq. 15. {ab } impose microscopic detailed balance. Finally, γ imposes the constraint of ensemble average r¯ of the dynamical variable rab (see Eq. 13). Differentiating with respect to kab and setting the derivative to zero and imposing the detailed balance constraint we have (see appendix A1)

kab =

µa µb Gab ya

(18)

where Pab rab + rba Gab = exp γ . xb 2

(19)

Note that the constraints in Eq. 14, Eq. 15, and Eq. 16 are not independent of each other. Indeed, as a result, the three sets of Lagrange multipliers {la }, {mb }, and {ab } that impose the constraints are also not independent of each other. Consequently, the total number of effective Lagrange multipliers {µa } is smaller than the number of Lagrange multipliers used to impose the constraints. For a specific value of the Lagrange multiplier γ, we determine the modified Lagrange multipliers µa by imposing the constraints given in Eq. 14. We have (15,16) X

kab = 1 ∀ a

(20)

b

⇒

X

Gab µb =

b

7

ya ∀ a. µa


(21)


Page 8 of 29

Eq. 21 can be reorganized by defining a non-linear operator D[¯ µ]a = ya /µa . We note that D[D[¯ µ]] = µ ¯. We have

G¯ µ = D[¯ µ] ⇒ D[G¯ µ] = µ ¯.

(22)

In other words, the vector µ ¯ of modified Lagrange multipliers can be numerically solved as a fixed point of Eq. 22. Note that the matrix G is symmetric since the transition probabilities {Pab } satisfy detailed balance with respect to the equilibrium distribution{xa }. Eq. 19 indicates that an important consequence of imposing detailed balance (15–17) is that the dynamical constraint rab appears in its symmetrized form (rab + rba )/2. From now onwards, for simplicity, we assume the constraint rab is already symmetrized in a and b; rab = rba for all a and b. In the limit dt → 0 in Eq. 9, the modified Lagrange multipliers {µa } can be solved analytically. In this limit, the transition probabilities kab of the updated Markov process are given by (see appendix A2) r kab = dt

yb ya

r

xa † ωab exp γrab if a 6= b xb

(23)

where (see appendix A2 for details) † rab =

rab + rba raa + rbb − . 2 2

(24)

Eq. 19 and Eq. 24 indicate that when we impose detailed balance and take the continuous time limit, we modify the dynamical constraint to make it symmetric in a and b and to have raa = 0 ∀ a. Mathematically, we perform the transformation given by Eq. 24. For brevity, from now onwards, we assume that this transformation has already been performed (unless specified otherwise). We drop the † superscript for simplicity. Before we proceed further, let us examine the consequences of the transformation in Eq. 24. Consider an antisymmetric dynamical quantity such that rab + rba = 0. A constraint which imposes a finite value of hrab i is clearly inconsistent with detailed balance.

8


Page 9 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


† = 0. Similarly, consider when Indeed, the transformation in Eq. 24 modifies rab to rab

rab = λ1 f (a) + λ2 g(b) (where λ1 and λ2 are constants and f and g are state-dependent functions) can be separated as a sum of two state-based constraints. Since we explicitly constrain the stationary distribution {ya }, we do not have additional freedom to constrain state-dependent quantities. Here too, the transformation modifies the constraint † = 0. to rab

From Eq. 23, the updated continuous time transition rates κab are given by kab κab = lim = dt→0 dt

r

yb ya

r

xa ωab exp (γrab ) xb

(25)

In Eq. 25, the transition rates κab describe an updated Markov process that is minimally biased with respect to the prior Markov process given by rates ωab (see Eq. 7) and a) has a prescribed equilibrium distribution {ya } and b) reproduces a prescribed dynamical average hrab i. Next we take the continuous space limit of Eq. 25 by substituting ωab given in Eq. 7. We have r κa,a+da =

ya+da D(a + da/2) da2 ya (26)

where

D(a + da/2) = D0 exp (γra,a+da )

(27)

is the updated diffusion coefficient at a + da/2 and ya ∝ e−βG(a) is the prescribed equilibrium distribution. We denote γra,a+da = γh(a). We have

D(a + da/2) = D0 exp (γh(a)) .

(28)

Below, we show how to explicitly derive D(a) from the functional form of the constraints.

9



Page 10 of 29

Comparing the ‘prior’ discretized transition rates {ωa,a+da } in Eq. 7 and the ‘updated’ transition rates {κa,a+da } in Eq. 26, we write the ‘updated’ Smoluchowski equation as ∂ ∂ p(a, t) = ∂t ∂a

βG(a) −βG(a) ∂ D(a)e e p(a, t) ∂a (29)

where

D(a) = D0 eγh(a)

(30)

is the position-dependent diffusion coefficient. Note that here we have used the C2 M2 approach to ‘update’ a ‘prior’ Smoluchowski equation rather than deriving the Smoluchowski equation from first principles as was done by Haken (22) . Before we further illustrate Eq. 29 with examples, we make a few observations. First, if we only update the equilibrium distribution (xa → ya ) and impose no additional dynamical constraint, the corresponding change in the Smoluchowski equation is simply changing its equilibrium distribution (see Eq. 46 and Eq. 29). In contrast, the imposition of trajectory-based constraints leads to a diffusion coefficient D(a) that depends on the position a. Second, a straightforward modification allows us to incorporate multiple dynamical constraints. Each constraint is associated with one Lagrange multiplier. In this case, the diffusion coefficient is given by

P

D(a) = D0 e

i

γi hi (a)

(31)

As an illustration for the recipe to calculate diffusion coefficients from trajectorybased constraints, let us also look at specific constraints. Consider rab = φ(a)φ(b) for some function φ of the position. The constraint hrab i represents the autocorrelation of the quantity φ along dynamical trajectories of the Markov process. After performing the

10


Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


transformation in Eq. 24, we have (omiting the † for brevity)

rab = −(φ(a) − φ(b))2 /2.

(32)

Thus from Eq. 30 0

ra,a+da = −1/2φ0 (a)da2 ⇒ D(a) = D0 e−γφ (a) ≈

2 /2

D0 . 1 + γφ0 (a)2 /2

(33)

The second approximation holds true when |γ| 1 or when φ(a) is slowly varying. Two notable examples of constraints are (1) a position-position autocorrelation function, and (2) a PMF-PMF autocorrelation along a stochastic trajectory. These constraints can be represented as habi and hF (a)F (b)i respectively. Here, a is the position coordinate and F (a) is the corresponding free energy. When we constrain the position-position autocorrelation, the updated diffusion coefficient does not depend on the position but simply takes a different value than the ‘prior’ diffusion coefficient. We have

D(a) = D0 e−γ .

(34)

In contrast, when we constrain the PMF-PMF correlation, we have

D(a) = D0 e−γF

0 (a)2 /2

D0 1 + γF 0 (a)2 /2 D0 = 1 + γfa2 /2

≈

(35) (36)

where we have identified fa = −F 0 (a) as the average force at position a. Position-dependent diffusion coefficients have been interpreted as effective corrections to lower dimensional projections of higher dimensional dynamics (23,24) . Specifically, Zwanzig (23) showed that the one dimensional diffusive dynamics along the length of a two dimensional channel with variable width w(a) is best described by a position-dependent

11



Page 12 of 29

diffusion coefficient

D(a) =

D0 . 1 + w0 (a)2 /12

(37)

From Eq. 33 and Eq. 37, this result can be interpreted within the maximum relative entropy framework as a position dependent diffusion coefficient arising from the dynamical constraint hrab i where rab = w(a)w(b) (width-width correlation along stochastic trajectories). Berezhkovskii and Szabo (24) considerably generalized the original work by Zwanzig and explicitly derived the formula for the position dependent diffusion coefficient for diffusion along a ‘slow’ dimension in a multi-dimensional system. In recent years, position-dependent diffusion coefficients have proved to be a very popular in studying lower-dimensional dynamics of complex biomolecules. For example, Best and Hummer (25,26) have studied the effective dynamics of protein folding along a one-dimensional reaction coordinates defined as the fraction of native contacts, Chodera and Pande (27) have studied unfolding of a DNA hairpin along the extension of the hairpin. In many such examples, the central goal is to infer the position-dependent diffusion coefficient from molecular dynamics data. Complementary to these studies, in this work we interpret position dependent diffusion coefficient as arising from trajectory-based constraints on Markovian dynamics.

Multidimensional problem: Illustration in two dimensions In the above development, we updated the Smoluchowski equation in one dimension. However, the method developed can be generalized to a multidimensional problem in a straightforward manner. We note that the rest of the manuscript can be read without this section. We illustrate the two dimensional derivation with a particle diffusing on a two dimensional landscape (see Fig. 2. For simplificty of notation, we assume that the en-

12


Page 13 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: Discretization of a two dimensional diffusion problem. A particle diffusing on a two dimensional lattice with spacing dx = dy in x and y directions. In a single time step, the particle can hop to one of its nearest neighbors shown in the figure.

ergy landscape and equilibrium probability distribution of the ‘prior’ process is flat; peq ([x, y]) = const. Let p([x, y]|t) denote the probability of observing the particle at position [x, y] at time t. The ‘prior’ dynamics of p([x, y]|t) is given by ∂2 ∂2 ∂ p([x, y]|t) = Dx 2 p([x, y]|t) + Dy 2 p([x, y]|t) ∂t ∂x ∂y (38)

In Eq. 38 Dx is the diffusion coefficient in the x direction and Dy is the diffusion coefficient in the y direction. Consider that we update the prior Markov model given by Eq. 38 by imposing a dynamical constraint hrab i as was done for the one dimensional case above. We also impose that the equilibrium distribution remains unchanged; peq ([x, y]) = const. For simplicity, as above, we assume that rab quantifies correlation in some quantity φ along dynamical trajectories. Mathematically, rab = φ(a)φ(b). As we show in the appendix A3, imposing a dynamical constraint introduces position dependent coefficient in the 2-dimensional problem as well. We have the updated diffusion

13



Page 14 of 29

coefficients Dxnew and Dynew :

Dxnew ([x, y]) = Dx exp −γ0   Dynew ([x, y]) = Dy exp −γ0

2 ∂ φ([x, y]) ∂x

!

2 2 

∂ φ([x, y]) ∂y

2

 

(39)

where γ0 is the modified Lagrange multiplier. Notably, updating diffusion coefficients based on dynamical constraints can introduce dynamical anisotropy. The update to the diffusion coefficient in the x− direction is different from the one in the y−direction.

Example application: overdamped oscillator

Figure 3: The equilibrium probability distribution of the overdamped harmonic oscillator. The equilibrium distribution of the Harmonic oscillator is given by a two dimensional normal distribution (see Eq. 41). We now illustrate an application of present method to an overdamped Harmonic oscillator. Consider a two-dimensional harmonic oscillator in equilibrium with its thermal surroundings and undergoing overdamped Langevin dynamics. The equations of motion of the 2 dimensional Harmonic oscillator are

p x(t) ˙ = Dx (−1.6x(t) + 0.4y(t)) + 2Dx η1 (t) p y(t) ˙ = Dx (0.4x(t) − 0.8y(t)) + 2Dy η2 (t)

(40)

where hηi (t)ηj (t0 )i = ∆ij δ(t − t0 ) for i, j = 1, 2. Here, ∆ij is the Kronecker delta function 14


Page 15 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and δ(t) is the Dirac delta function. The diffusion constants are Dx = 5 and Dy = 1 units and the inverse temperature is β = 1. The equilibrium distribution peq (x, y) is described by a two dimensional Gaussian distribution,

peq (x, y) ∝ e−0.8x

2 +0.4xy−0.4y 2

(41)

and is shown in Fig. 3. We simulate Eqs. 40 with a discretized Langevin dynamics scheme with dt = 0.001 units. In Fig. 4 we show the normalized autocorrelation function

C(t) =

hx(t)x(0)i − hx(0)i2 hx(0)2 i − hx(0)i2

(42)

of the one dimensional projection x(t) of the two dimensional dynamics. The autocorrelation decays with two time scales, a fast decay for t < 1000dt and a slower decay after t > 1000dt.

Figure 4: Autocorrelation function of the two dimensional dynamics reveals two time scales. The normalized autocorrelation function C(t). The dashed red line indicates the faster of the two time scales in the autocorrelation function. The error bars represent standard deviation in mean estimated from 500 independent calculations of C(t).

How do we describe the effective stochastic dynamics of x(t)? The marginal equilibR rium distribution peq (x) = peq (x, y)dy is (see Appendix ) 2

peq (x) ∝ e−0.7x .

15


(43)


Page 16 of 29

As a first guess, we write down the simplest Smoluchowski equation that relaxes to this equilibrium distribution. We have ∂ ∂ p(x, t) = Dx ∂t ∂x

h i −0.7x2 ∂ 0.7x2 e e p(x, t) . ∂x

(44)

It is well known that the autocorrelation function described by Eq. 44 decays exponentially with a single time constant (28) . As a result Eq. 44 cannot capture the essential features of x(t) dynamics. Can we model the dynamics with a position-dependent diffusion coefficient? We impose two kinetic constraints noted above (see Eq. 34 and Eq. 36). The first constraint corresponds to the position-position autocorrelation and the second constraint corresponds to the PMF-PMF autocorrelation. The corresponding position dependent diffusion coefficient is given by (see Eq. 33)

D(x) = Dx exp −γ1 − γ2 x2 .

(45)

Here, γ1 and γ2 are Lagrange multipliers that relate to the dynamical constraint hrab i = habi and hrab i = hF (a)F (b)i respectively. As discussed above (see Eq. 34 and Eq. 36), the Lagrange multiplier γ1 allows us to adjust the overall diffusion constant. The Lagrange multiplier γ2 > 0 slows down diffusion in regions of the x−space where PMF changes most rapidly. Specifically, diffusion coefficient gets smaller as |x| increases. The Smoluchowski equation with a position dependent diffusion coefficient is given by ∂ ∂ p(x, t) = ∂t ∂x

h i −0.7x2 ∂ 0.7x2 D(x)e e p(x, t) ∂x (46)

In Fig. 5 we plot the normalized autocorrelation function C(t) (red dots) as predicted by the stochastic dynamics described by Eq. 46 and compare it to the autocorrelation function shown in Fig. 4 (black line). We have used γ1 = −0.4, γ2 = 0.9. These two

16


Page 17 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5: Effective dynamics with a position dependent diffusion coefficient captures the autocorrelation function of the two dimensional harmonic oscillator. The comparison between the normalized autocorrelation function C(t) estimated using Eq. 46 (red) and Eq. 40 (black). The effective one dimensional dynamics of Eq. 46 captures the two time scales in x(t) dynamics. The error bars represent standard deviation in mean estimated from 500 independent calculations of C(t). The inset shows the position dependence of the diffusion coefficient D(x) on x. values of the Lagrange multipliers were found out by trial by matching both the average mean squared displacement (position-position correlation) and PMF-PMF correlation. We have a discretization time step of dt = 0.001 and used the Ito convention to simulate the position dependent diffusion coefficient (see appendix A3). The inset shows the dependence of the diffusion coefficient on x. Notably, incorporating a position dependent diffusion coefficient in Eq. 45 allows us to capture the two time scales observed in x(t)dynamics with sufficient accuracy. The effective one dimensional dynamics described by Eq. 46 can also predict other trajectory-based dynamical quantities without any adjustible parameters. In Fig. 6, we show the agreement between the mean first pasage time τp to reach x(t) = xf for the first time when starting from x(0) = 0 as a function of xf .

Conclusion We have described a method for updating a Smoluchowski equation based on observables captured in state- and path-dependent constraints. We showed how this can be expressed in terms of position-dependent diffusion coefficients. We illustrated by considering the effective one-dimensional dynamics of a two-dimensional overdamped harmonic oscilla-

17



Figure 6: Position dependent diffusion coefficient captures the mean first passage time. Comparison of the estimated mean first passage time τp to reach x(t) = xf for the first time when starting from x(0) = 0 as a function of xf in the 2 dimensional dynamics (black, Eq. 40) and the effective one dimensional dynamics (red circles, Eq. 46).

tor, with a position dependent diffusion coefficient D(q). The present approach is not limited to updating Markov models that are continuous time and continuous space; this approach can also handle discrete-time discrete-space models using Eq. 18. Similarly, Eq. 25 illustrates how to update a continuous time Markov process. Acknowledgments: KD appreciates support from the National Science Foundation (grant number 1205881)

References (1) Chodera, J. D.; Noé, F. Markov state models of biomolecular conformational dynamics. Curr. Opin. Struc. Biol. 2014, 25, 135–144. (2) Gillespie, D. T. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 1977, 81, 2340–2361. (3) Paulsson, J. Summing up the noise in gene networks. Nature 2004, 427, 415–418. (4) Dixit, P. D.; Pang, T. Y.; Maslov, S. Recombination-driven genome evolution and stability of bacterial species. Genetics 2017, genetics–300061. (5) Sawle, L.; Ghosh, K. Convergence of molecular dynamics simulation of protein native

18


Page 18 of 29

Page 19 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


states: Feasibility vs self-consistency dilemma. J. Chem. Th. Comp. 2016, 12, 861– 869. (6) Dixit, P. D.; Wagoner, J.; Weistuch, C.; Pressé, S.; Ghosh, K.; Dill, K. A. Perspective: Maximum caliber is a general variational principle for dynamical systems. J. Chem. Phys. 2018, 148, 010901. (7) Pitera, J. W.; Chodera, J. D. On the use of experimental observations to bias simulated ensembles. J. Chem. Th. Comp. 2012, 8, 3445–3451. (8) Copeland, R. A. The drug-target residence time model: a 10-year retrospective. Nature Reviews. Drug discovery 2016, 15, 87. (9) Tiwary, P.; Mondal, J.; Berne, B. How and when does an anticancer drug leave its binding site?. Science Advances 2017, 3, e1700014. (10) Shell, M. S. The relative entropy is fundamental to multiscale and inverse thermodynamic problems. J. Chem. Phys. 2008, 129, 144108. (11) Roux, B.; Weare, J. On the statistical equivalence of restrained-ensemble simulations with the maximum entropy method. J. Chem. Phys. 2013, 138, 02B616. (12) Woolf, T. B.; Roux, B. Conformational flexibility of o-phosphorylcholine and ophosphorylethanolamine: a molecular dynamics study of solvation effects. J. Am. Chem. Soc. 1994, 116, 5916–5926. (13) Hummer, G. Position-dependent diffusion coefficients and free energies from Bayesian analysis of equilibrium and replica molecular dynamics simulations. New J. Phys. 2005, 7, 34. (14) Bicout, D.; Szabo, A. Electron transfer reaction dynamics in non-Debye solvents. J. Chem. Phys. 1998, 109, 2325–2338. (15) Dixit, P. D.; Dill, K. A. Inferring microscopic kinetic rates from stationary state distributions. J. Chem. Th. Comp. 2014, 10, 3002–3005.

19



(16) Dixit, P. D.; Jain, A.; Stock, G.; Dill, K. A. Inferring transition rates of networks from populations in continuous-time markov processes. J. Chem. Th. Comp. 2015, 11, 5464–5472. (17) Dixit, P. D. Stationary properties of maximum-entropy random walks. Phys. Rev. E 2015, 92, 042149. (18) Wan, H.; Zhou, G.; Voelz, V. A. A maximum-caliber approach to predicting perturbed folding kinetics due to mutations. J. Chem. Th. Comp. 2016, 12, 5768–5776. (19) Zhou, G.; Pantelopulos, G. A.; Mukherjee, S.; Voelz, V. A. Bridging microscopic and macroscopic mechanisms of p53-MDM2 binding with kinetic network models. Biophys. J. 2017, 113, 785–793. (20) Olsson, S.; Wu, H.; Paul, F.; Clementi, C.; Noé, F. Combining experimental and simulation data of molecular processes via augmented Markov models. Proc. Natl. Acad. Sci. 2017, 114, 8265–8270. (21) Rached, Z.; Alajaji, F.; Campbell, L. L. The Kullback-Leibler divergence rate between Markov sources. IEEE Trans. Inf. Th. 2004, 50, 917–921. (22) Haken, H. A new access to path integrals and fokker planck equations via the maximum calibre principle. Zeit. Phys. B. Cond. Matt. 1986, 63, 505–510. (23) Zwanzig, R. Diffusion past an entropy barrier. J. Phys. Chem. 1992, 96, 3926–3930. (24) Berezhkovskii, A.; Szabo, A. Time scale separation leads to position-dependent diffusion along a slow coordinate. J. Chem. Phys. 2011, 135, 074108. (25) Best, R. B.; Hummer, G. Diffusive model of protein folding dynamics with Kramers turnover in rate. Phys. Rev. Lett. 2006, 96, 228104. (26) Best, R. B.; Hummer, G. Coordinate-dependent diffusion in protein folding. Proc. Natl. Acad. Sci. 2010, 107, 1088–1093.

20


Page 20 of 29

Page 21 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(27) Chodera, J. D.; Pande, V. S. Splitting probabilities as a test of reaction coordinate choice in single-molecule experiments. Phys. Rev. Lett. 2011, 107, 098102. (28) Zwanzig, R. Nonequilibrium statistical mechanics; Oxford University Press, 2001.

21



Page 22 of 29

Appendix A1: Imposing detailed balance in discrete time Markov processes We start with Eq. 17 in the main text. We have the Caliber,

C = −

X

ya kab log kab +

a,b

X

ya kab log Pab

a,b

! +

X

+

X

la

X

a

ya kab − ya

! +

b

X

mb

X

ya kab − yb

a

b

! ab (ya kab − yb kba ) + γ

X

a,b

ya kab rab − r¯ .

(47)

a,b

In Eq. 47, Lagrange multipliers {la } impose constraints in Eq. 14. {mb } impose constraints in Eq. 15. {ab } impose microscopic detailed balance. Finally, γ imposes the constraint of ensemble average r¯ of the dynamical variable rab (see Eq. 13). Differentiating with respect to kab and setting the derivative to zero,

log

kab = la − 1 + mb + δab + γrab Pab

(48)

where δab = ab − ba . We have

kab = τa λb ρab Wab

(49)

where τa = exp(la − 1), λb = exp(mb ), ρab = exp(δab ), and Wab = Pab exp (γrab ). We impose detailed balance, ya kab = yb kba , to evalulate ρab . We have

ya τa λb ρab Wab = yb τb λa ρba Wba s rba −rab yb τb λa Pba ⇒ ρab = × e(γ 2 ) ya τa λb Pab s rba −rab yb τb λa xa = × e(γ 2 ) ya τa λb xb

(50) (51) (52)

The last equality is a result of the fact that the Markov chain described by transition

22


Page 23 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


probabilities Pab obeys detailed balance with respect to the stationary distribution xa , xa Pab = xb Pba . Substituting ρab in Eq. 52 into Eq. 49, s

rba −rab yb τb λa xa × e(γ 2 ) × τa λb Pab eγrab ya τa λb xb p 1p Pab (γ rba +r ab ) 2 = e τa λa xa ya × τb λb xb yb ya xb

kab =

(53)

(54)

We substitute

√

τa λa xa ya = µa and

Pab (γ e xb

rba +rab 2

) = G and we obtain Eq. 18 in the main ab

text.

A2: Deriving transition rates for the continuous time Markov process In the main text, we claimed that the transition rates for the maximum entropy continuous time Markov process with a updated equilibrium distribution xa → ya and after imposing additional dynamical constraints is given by Eq. 55, kab κab = lim = dt→0 dt

r

yb ya

r

xa ωab exp (γrab ) . xb

(55)

Here, we prove this assertion. Let us consider the equation

D[G¯ µ] = µ ¯

(56)

where

G=

Pab exp (γrab ) xb

(57)

where rab is assumed to be symmetric in a and b; rab = rba ∀ a and b. Plugging the

23



Page 24 of 29

transition probabilities Pab in Eq. 9 in Eq. 19 and 22, we have ωab exp(γrab ) if a 6= b xb P 1 − dt b6=a ωab = exp(γraa ) xa

Gab = dt

(58)

Gaa

(59)

Thus,

G = J + dt∆

(60)

where J is a diagonal matrix with Jaa = exp(γraa )/xa and ωab exp(γrab ) if a 6= b xbP b6=a ωab = − exp(γraa ) xa

∆ab =

(61)

∆aa

(62)

Substituting Eq. 61 and 62 in Eq. 22, we have

D[(J + dt∆)¯ µ] = µ ¯

(63)

Note that if µ ¯ is a solution of Eq. 63, we have

kab = dt

µa µb Gab if a 6= b. ya

(64)

Thus, we need to find µ ¯ only till the zeroth order in dt as dt → 0. We have ya P = µa Jaa µa + dt b ∆ab µb

(65)

The solution for µ as dt → 0 to the zeroth order of Eq. 65 is given by r µa =

ya Jaa

24


(66)

Page 25 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Substituting this value of µa in Eq. 18, we have r kab = dt

yb ya

r

xa † ωab exp γrab xb

(67)

where † rab = rab −

raa + rbb 2

(68)

is a transformed version of the dynamical constraint rab such that raa = 0 ∀ a.

A3: Two dimensional diffusion We start with the prior dynamics ∂ ∂2 ∂2 p([x, y]|t) = Dx 2 p([x, y]|t) + Dy 2 p([x, y]|t) ∂t ∂x ∂y (69)

In Eq. 38 Dx is the diffusion coefficient in the x direction and Dy is the diffusion coefficient in the y direction. Discretizing the space in steps of dx = dy in the x and y direction respectively and discretizing time in steps of dt (omiting the time dependence for brevity), dp([x, y]) ≈ ω([x−dx,y]),([x,y]) p([x − dx, y]) dt + ω([x+dx,y]),([x,y]) p([x + dx, y]) + ω([x,y−dy]),([x,y]) p([x, y − dy]) + ω([x,y+dy]),([x,y]) p([x, y + dy]) − p([x, y])ω[x,y],[x,y]

25


(70)


Page 26 of 29

where

ω[x,y],[x,y] = ω([x,y]),([x−dx,y]) + ω([x,y]),([x+dx,y]) + ω([x,y]),([x,y−dy]) ω([x,y]),([x,y+dy])

(71)

In Eq. 70 ωa,b denotes the transition rate of going from a to b. We have ωa,b = 0 if a and b are not nearest neighbors on the lattice. From Eq. 38 in the main text, we can write ω[x,x±dx,y],[x,y] = Dx dx2 , P[x,y±dy],[x,y] = Dy dy 2 , and so on. Next, we impose a dynamical constraint hrab i = hφ(a)φ(b)i (see Eq. 32). Consider two points a = [x, y] and b = [w, u]. First, we explicitly carry out the transformation in Eq. 24. We write (omitting the † for brevity)

r[x,y],[w,u] = (φ([x, y]) − φ([w, u]))2

(72)

From Eq. 27 and Eq. 26 We can write the updated transition rates

κ[x,y],[x+dx,y] = Dx dx2 dt exp γr[x,y],[x+dx,y] and κ[x,y],[x,y+dy] = Dy dy 2 dt exp γr[x,y],[x,y+dy] .

(73)

Other transition probabilities can be written down similarly by recognizing that r[x,y],[u,w] = r[u,w],[x,y] . We can further simplify Eq. 73,

r[x,y],[x+dx,y] r[x,y],[x,y+dy]

2 ∂ φ([x, y]) and ∂x 2 dy 2 ∂ ≈ − φ([x, y]) 2 ∂y dx2 ≈ − 2

26


(74) (75)

Page 27 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Consequently,

κ[x,y],[x+dx,y] = Dx dx2 exp −γ0 

2 ∂ φ([x, y]) ∂x 2

 κ[x,y],[x,y+dy] = Dx dx2 exp −γ0

! ,

2 

∂ φ([x, y]) ∂y

2

  (76)

In Eq. 76 have recognized γ0 = γ/dx2 . Finally, the position dependent diffusion coefficients are given by

Dxnew ([x, y]) = Dx exp −γ0   Dynew ([x, y]) = Dy exp −γ0

2 ∂ φ([x, y]) ∂x

!

2 2 

∂ φ([x, y]) ∂y

2

 

(77)

Details of the Langevin dynamics Let us start with a Smoluchowski equation with a position dependent diffusion coefficient. ∂ ∂ p(x, t) = ∂t ∂x

h i 0.7x2 −0.7x2 ∂ e D(x)e p(x, t) ∂x (78)

There are multiple ways to map this Smoluchowski equation to a Langevin equation. The two popular approaches are the Ito approach and the Stratonovich approach. Both approaches lead to the same equilibrium distribution and have the same dynamics. The time-discretized Langevin equation with Ito convention is given by

x(t + dt) ≈ x(t) − 1.4D(x(t))x(t)dt + ρ dD(x) + |x=x(t) dt dx

p

2D(x(t))dt (79)

where ρ is a normally distributed random number with mean 0 and standard deviation

27



1. As above, we use dt = 0.001 units.

28


Page 28 of 29

Page 29 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23



Correcting Equilibrium Markov Models - American Chemical Society

Recommend Documents