Inferring Microscopic Kinetic Rates from Stationary State Distributions.

We present a principled approach for estimating the matrix of microscopic transition probabilities among states of a Markov process, given only its stationary state population distribution and a single average global kinetic observable. We adapt Maximum Caliber, a variational principle in which the path entropy is maximized over the distribution of all possible trajectories, subject to basic kinetic constraints and some average dynamical observables. We illustrate the method by computing the solvation dynamics of water molecules from molecular dynamics trajectories.

Introduction

We propose a method for inferring the kinetic rate matrix for stochastic systems for which the steady state popuations are known. The types of systems of interest include protein folding,[1−3] ion channels,[4] molecular motors,[5] the evolutionary dynamics of protein sequences,[6] the collective firing patterns of neurons,[7] or noisy gene expression.[8] In such studies, stationary probabilities {p} of N stable states {i} are often known or estimated from experimental data,[6,7] detailed molecular mechanics calculations,[3] coarse grained theories,[1,4] or from Maximum Entropy.[6,7,9−12] Moreover, suppose further that we know some average global dynamical quantity, such as the average mutation rate (number of amino acid changes per unit time), the average current through an ion channel, or the average number of neurons that change their state from spiking to resting (and vice versa) per unit time step. How can we infer the full distribution of microscopic state-to-state rates, given just the stationary-state populations and one or more average overall rate quantities?

More specifically, we are interested in a principled way to solve the following under-determined “inverse” kinetics problem. Consider a stationary and irreducible Markov process among i = 1,2,3,...,N states. Suppose you know the following: (a) the stationary state probability distribution, {p} of the occupancies of those states, and (b) the value ⟨w⟩ of some dynamical observable w averaged over the ensemble of stationary state trajectories. From these N + 1 quantities, we want to infer the N × N microscopic transition probabilities,between those states.

While the stationary distribution indeed constrains the transition probability matrix, it does not uniquely determine it. Here, we develop a procedure based on the principle of Maximum Caliber, a variant of the principle of Maximum Entropy, that is applicable to dynamical processes.[13,14] This variational principle allows us to uniquely infer ∼N2 transition probabilities that are consistent with ∼N imposed constraints.

First, we define the path entropy, S, over a given ensemble {Γ} of trajectories Γ ≡ ... → i → j → k → l... asMaximum Caliber is a variational principle that chooses a unique probability distribution {P(Γ)} over trajectories from all possible candidate distributions as the one that maximizes the path entropy while otherwise satisfying the relevant stationary and dynamical constraints.[13−16]

Consider an ensemble of stationary state trajectories {Γ} (see above) having a total time duration T. Restricting our attention to first-order Markov trajectories allows us to greatly simplify the path entropy and carry out our analysis in terms of trajectories of single steps, i → j.[2,17−19] The Markov property implies that the probability of any particular trajectory Γ can be expressed in terms of the transition probabilities {k},The path entropy of the above ensemble is directly proportional to the total duration T of the trajectory. The path entropy per unit time is given by[20,21]

The microscopic transition probabilities of any Markov process are subject to two types of constraints. First, from state i at time t, the system must land somewhere at time t + dt. Second, a system in state j at time t + dt must arrive from somewhere, so

Third, we require one additional constraint that is global, that is, averaged over the entire ensemble of trajectories. We fix the path ensemble average of some dynamical quantity w. The average ⟨w⟩Γ over any given stationary trajectory Γ is given by

The path ensemble average ⟨w⟩ isSince Γ is a stationary state trajectory, the path ensemble average ⟨w⟩ of eq 7 simplifies to

Maximization of the path entropy subject to these three constraints can be expressed equivalently in terms of maximization of a quantity called the Caliber :[14]where γ is the Lagrange multiplier associated with the constraint ⟨w⟩ and {a} and {l} enforce the to-somewhere constraint and the from-somewhere constraint, respectively.

To solve for the matrix k of transition probabilities, we take the derivative of the Caliber with respect to k and equate it to zero. This giveswhere we have made the substitutions: e = β/p and e = λ. The values k* are the transition probabilities that satisfy the constraints and otherwise maximize the caliber. For simplicity of notation, we drop the superscript * in the remainder of this paper, that is, k* ≡ k.

In this problem, the values of p are given. To compute the k values, we first must determine the values of the Lagrange multipliers β, λ, and γ. We do so by substituting the constraint relations mentioned above.

Determining the Lagrange Multipliers

For a given value of γ, the modified Lagrange multipliers β and λ are determined by satisfying the to-somewhere and from-somewhere conditions indicated above. From eqs 5where e–γ = W. Equation 11 can be simplified if we define a nonlinear operator D over column vectors x̅ = [x1,x2,...]T as D[x̅] = p/x. We havewhere λ̅ = [λ1,λ2,...]T and β̅ = [β1,β2,...]T are the column vectors of Lagrange multipliers.

For a particular value of the Lagrange multiplier γ, eqs 12 can be numerically and self-consistently solved for {β} and {λ}. In practice, we choose an appropriate γ by first constructing transition probabilities {k} for multiple values of γ (see eq 10) and choosing the value of γ that satisfieswhere ⟨w⟩ is the prescribed value of the ensemble average of the dynamical quantity w.

Detailed Balance

So far, the treatment above is generally applicable to nonequilibrium systems. However, if we are interested in restricting our attention to systems in thermodynamic equilbrium, we can impose an additional constraint that the system must satisfy detailed balance, pk = pk. In this case, the Caliber can be expressed asHere, d are Lagrange multipliers that impose the detailed balance condition. Differentiating with respect to k and setting the derivative to zero, we find (see eq 10)where δ = e = 1/δ.

Now we have to determine the modified Lagrange multipliers β, λ, δ, and γ from the imposed constraints. Let us first impose detailed balance to determine δ. We haveSince δ = 1/δ, we haveIn eq 17, identifying ρ = (βλ)1/2,where e = Wsym is the symmetrized form of W = e–γ (see above). It is easy to see that k in eq 18 satisfies detailed balance. The Lagrange multipliers ρ̅ = [ρ1,ρ2,...,ρ] and γ are determined by the same procedure as described above. We haveIn other words, ρ̅ is the solution of the nonlinear equationand γ is determined by adjusting the path ensemble average

An Illustration: Computing the Dynamics of a Solvation Shell from Simulated Populations

We now illustrate how the present MaxCal method can be used to take a stationary-state distribution and a global dynamical constraint and to infer microscopic kinetics. Consider a shell of solvating water molecules surrounding a single water molecule. The number, n(t), of water molecules in the hydration shell is a quantity that fluctuates with time t (see Figure Figure 1). We want to compute how fast the water molecules enter or exit the solvation shell. If the time interval dt is small, n(t) and n(t + dt) will be statistically correlated. Here, we construct a Markov process to model the time series {n(t)}. We will require the Markov process to reproduce (a) the stationary distribution p(n) that is observed in molecular dynamics simulations, and (b) the average change in occupancy Δ per time step of duration dt, a path ensemble average (see Appendix for details of molecular simulation). While the choice of dynamical constraint(s) remains somewhat arbitrary, it is validated only a posteriori. In other words, the constrain(s) is a modeling aspect of any maximum entropy method.[22] We havewhere n(t) = i and n(t + dt) = j.

Figure 1

Hydration shell (black circle) around a central water molecule (blue disc) is dynamically populated by other water molecules in the bulk solvent medium (red discs). The probability, p(n) that the hydration shell has exactly n water molecules is a key quantity in determining the solvation free energy of liquid water.[23,24]

Since the system should satisfy detailed balance, the transition probability k for a transition n(t) = i → n(t + dt) = j is given by (see eq 18),

For a given value of γ, we determine the Lagrange multipliers ρ from eqs 20 above. In order to determine the Lagrange multiplier γ which dictates the rate of transition between states, we first construct Markov processes for different values of γ. Panel B of Figure Figure 2 shows that the path ensemble average of the change in occupation number per unit time step Δ is exponentially decreasing with γ. From trajectories sampled at every 5 fs from the MD simulation, we find that experimental trajectory average Δexpt = ≈ 0.0629 which corresponds to γ ≈ 3.29. From here onward, we use γ = 3.29 and construct the transition probabilities {k} (see eq 23). Note that the path ensemble average Δ and consequently the Lagrange multiplier γ, depend on the time interval dt between two observation (dt = 5 fs here).

Figure 2

Panel A: The stationary distribution p(n) of the number of water molecules in the hydration shell of radius r = 3.2 Å of a reference water molecule. Panel B: The dependence of the ensemble average of change in water occupancy number Δ = ⟨|n(t + dt) – n(t)|⟩ on the Lagrange multiplier γ. We see that Δ depends exponentially on γ. A higher γ implies slower dynamics and vice versa. We choose γ = 3.29 to match the observed Δ ≈ 0.0629 in the molecular dynamics simulation.

From the Markov process constructed with γ = 3.29 (see above), we now compute various dynamical quantities: (a) the probability P of jump size d, (b) the occupancy autocorrelation ⟨δ⟩, and (c) the transition probabilities k, and we compare to those obtained directly from the MD simulation trajectory. In general, the MaxCal method will be of value when transitions are hard to simulate, such as for large kinetic barriers. Here, we are just illustrating with a toy problem for which we can determine the transition probabilities independently from the simulations.

From the long simulation trajectory, the probability P of jump size is estimated as the histogram of d = n(t + dt) – n(t). Here d could be both positive and negative. P is given byThe normalized occupancy autocorrelation is simply the joint probability that n(t) and n(t + τ) are equal. It is given bywhere K(τ) = kτ is the τth power of the matrix of transition probabilities {k}.

In Figure 3 we plot P, ⟨δ⟩, and {k} estimated from the molecular dynamics trajectory and compare them to our predictions from Markov modeling. Even though we constrained only the mean value Δ = ⟨|d|⟩ of P, the modeled Markov process captures the entire distribution P with high accuracy. Similarly the occupancy correlation ⟨δ⟩ is also reproduced with high accuracy even though we did not utilize any information about it when inferring the transition probabilities of the Markov process. Moreover, our modeling also accurately captures the individual transition probabilities {k} over 4 orders of magnitude.

Figure 3

Panel A: The probability P of jump size estimated from the trajectory derived from MD simulation (red squares) is compared to the one predicted using the transition probabilities of the Markov process (dashed black line). Panel B: The normalized occupancy autocorrelation ⟨δ⟩ as estimated from the MD trajectory and as predicted from the transition probabilities of the Markov process. Panel C: We directly compare the transition probabilities k for the probability of transition i → j empirically obtained from the MD trajectory to the ones predicted by using eq 23.

Discussion and Summary

We have presented here a variational approach that computes N × N microscopic transition probabilities of a Markov process, given only knowledge of a stationary state population distribution and one trajectory-averaged dynamical property. In this approach, we maximize the path entropy subject to constraints; that is, we maximize the Caliber. We show that this method gives correct values of dynamical quantities in an example of molecular dynamics simulations of a water solvation shell around a single water molecule. This method may be useful for analyzing single-molecule experiments such as on ion channels,[4] dynamics of neuron firing,[7] and the dynamics of protein-sequence evolution,[6] for example.