Microscopic reweighting for non-equilibrium steady states dynamics

TL;DR

This paper develops a method to reweight and analyze non-equilibrium steady states in simulations by extending equilibrium statistical reweighting techniques using a maximum path entropy approach and stochastic thermodynamics.

Contribution

It introduces a novel framework for reweighting non-equilibrium steady states based on path entropy and Markovian pathways, bridging a gap in simulation analysis.

Findings

Reweighting in and out of equilibrium using path probabilities.

Identification of an invariant quantity analogous to density of states.

Systematic construction of pathways through Markovian transitions.

Abstract

Computer simulations generate trajectories at a single, well-defined thermodynamic state point. Statistical reweighting offers the means to reweight static and dynamical properties to different equilibrium state points by means of analytic relations. We extend these ideas to non-equilibrium steady states by relying on a maximum path entropy formalism subject to physical constraints. Stochastic thermodynamics analytically relates the forward and backward probabilities of any pathway through the external non-conservative force, enabling reweighting both in and out of equilibrium. We avoid the combinatorial explosion of microtrajectories by systematically constructing pathways through Markovian transitions. We further identify a quantity that is invariant to dynamical reweighting, analogous to the density of states in equilibrium reweighting.

Figures

Equations

$$\mathcal{C}[p(x)]=-\int \mathrm{d}xp(x)\ln \left(\frac{p(x)}{q(x)}\right).$$

$${\mathcal{C}}_{\text{equ}}=-\int \mathrm{d}xp(x)\ln \left(\frac{p(x)}{q(x)}\right)-\zeta \left(\int \mathrm{d}xp(x)-1\right)$$

$$-\beta \left(\int \mathrm{d}xp(x)E(x)-⟨E⟩\right),$$

$$p(x)=\frac{q(x)}{\tilde{Z}(\beta )}\exp \left(-\beta E(x)\right),$$

$$\mathcal{C}=-\sum _{i,j}{\pi }_i{p}_{ij}\ln \frac{p_{ij}}{q_{ij}},$$

$$\frac{\mathcal{P}[\Gamma (+t)|f(+t)]}{\mathcal{P}[\bar{\Gamma }(-t)|\bar{f}(-t)]}=\exp \left(-\beta Q|\Gamma (+t)|f(+t)]\right),$$

$$⟨\Delta S_{ij}⟩=\ln \frac{p_{ij}}{p_{ji}},$$

$${\mathcal{C}}_{\text{dyn}}=\underset{\text{Caliber}}{\underbrace{-\sum _{i,j}{\pi }_i{p}_{ij}\ln \frac{p_{ij}}{q_{ij}}}}+\underset{\text{Normalization}}{\underbrace{\sum _i{\mu }_i\left(\sum _j{p}_{ij}-1\right)+\zeta \sum _i({\pi }_i-1)}}+\underset{\text{Global Balance}}{\underbrace{\sum _j{\nu }_j\left(\sum _i{\pi }_i{p}_{ij}-{\pi }_j\right)}}+\underset{\text{Local Balance}}{\underbrace{\sum _{ij}^{i<j}{\pi }_i{\alpha }_{ij}\left(\ln \left(\frac{p_{ij}}{p_{ji}}\right)-\Delta S_{ij}\right)}},$$

$$p_{ij}=q_{ij}\exp \left(\frac{c_i+c_j}{2}+\frac{\Delta S_{ij}-\Delta S_{ij}^q}{2}+\zeta \right),$$

$$1=\sum _j{p}_{ij}=\sum _j\sqrt{q_{ij}{q}_{ji}}\exp \left(\frac{c_i+c_j+2\zeta }{2}+\frac{\Delta S_{ij}}{2}\right).$$

$$0=-\frac{\partial U(x)}{\partial x}-\gamma \dot{x}+\sqrt{2\gamma k_{\mathrm{B}}T}R(t)+f,$$

$$\Delta S[\Gamma (t)]=\int \mathrm{d}t\frac{\mathbf{F}\cdot \dot{\mathbf{\Gamma }}}{T},$$

$$\Delta S_{ij}(\{{\mathbf{x}}_k\})\approx \frac{U(x_i)-U(x_j)+{\sum }_k({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)\cdot \mathbf{f}}{k_{\mathrm{B}}T}.$$

$$\Delta S_{ij}\approx \frac{U(x_i)-U(x_j)+(x_j-x_i)f}{k_{\mathrm{B}}T}.$$

$$p_{ij}=\frac{I_{ij}}{Z_{ij}}\exp \left(\frac{\Delta S_{ij}}{2}\right),$$

Full text

Many chemical and biological processes are influenced by external driving forces and operate away from equilibrium—examples include colloidal particles, biopolymers, enzymes, or molecular motors Seifert (2008). Despite our current lack of a universal theory for statistical mechanics off equilibrium Dougherty (1994), computer simulations can complement experiments by providing microscopic insight into these complex processes. Unfortunately current computational power often prevents molecular simulations from reaching the experimentally-relevant time scales, or alternatively, obliges them to operate at artificially-large driving forces Perilla et al. (2015). The latter motivates a formalism to reweight dynamics across off-equilibrium conditions.

When dealing with systems in equilibrium, Ferrenberg and Swendsen introduced a statistical-reweighting procedure to infer information about a system when sampled at another state point Ferrenberg and Swendsen (1988, 1989). It requires microscopic information at fixed thermodynamic conditions, e.g., temperature, collected by computer simulations or experiments. A probability associated with each microstate is reweighted according to physical relationships linking the initial and final thermodynamic conditions. Reweighting can be conducted arbitrarily far from the initial state, provided it is sufficiently sampled.

Equilibrium reweighting has led to a number of developments in the field, from estimating accurate free energies Kumar et al. (1992) to building more robust Markov state models Chodera et al. (2011); Wu et al. (2014). In this Letter, we generalize reweighting to dynamical processes by replacing microstates with microtrajectories. The proposed methodology, valid for non-equilibrium steady state (NESS) systems, employs a maximum path entropy formalism while generalizing the standard detailed balance relation.

Jaynes’ maximum entropy approach offers a general variational principle to understand macroscopic phenomena from microscopic knowledge of a statistical system Jaynes (1957). This information-theoretic method regards entropy as the measure of uncertainty of the system. Consider a coordinate $x$ of a system with unknown probability distribution, $p(x)$. We further define another distribution, $q(x)$, used as a prior on $p(x)$. The most likely representation of $p(x)$ can be found by minimizing the cross-entropy functional

[TABLE]

This quantity was shown to fulfill the axioms for an uncertainty measure Shore and Johnson (1980). Setting a uniform prior (i.e., $q(x)=\text{const.}$) reduces to the well-known Shannon entropy and minimizing the cross entropy in this case is equivalent to maximizing the Shannon entropy.

According to Jaynes, a system would maximize the number of microscopic realizations compatible with a certain macroscopic state, linking the two scales via constraints. For instance, working in the canonical ensemble will lead to a constraint on the average energy $\langle E\rangle$

[TABLE]

where $\zeta$ and $\beta$ are Lagrangian multipliers, controlling the normalization of probabilities and the fixed average energy. Minimization of $\mathcal{C_{\text{equ}}}$ with respect to $p(x)$ yields

[TABLE]

where the partition function $\tilde{Z}(\beta)=\int\mathop{}\!\mathrm{d}x\,q(x)\exp\left(-\beta E(x)\right)$ becomes a normalization constant and the Lagrange multiplier $\beta$ is identified with the inverse temperature, $\beta^{-1}=k_{\rm B}T$. This approach naturally lends itself to reweighting: Given a reference distribution $q(x)$ sampled at inverse temperature $\beta^{\prime}$, microscopic information at a third inverse temperature $\beta^{\prime\prime}=\beta+\beta^{\prime}$ is inferred through the calculation of $p(x)$. This reweighting becomes exact under full knowledge of the density of states function $\Omega(x)=q(x)\exp\left(\beta^{\prime}E(x)\right)$.

The maximum entropy formalism has been generalized to the study of dynamical systems by working with microtrajectories—an approach called Maximum Caliber Jaynes (1980). It was shown to recover known off-equilibrium relations Dixit et al. (2018) and predict dynamical pathways in NESSs correctly when supplied with appropriate constraints Agozzino and Dill (2019). Conceptually the approach follows the same scheme as the previous derivation of equilibrium reweighting: The most likely microtrajectories maximize the path entropy function subject to physical constraints.

The following discusses a rich and relevant subset of non-equilibrium processes: non-equilibrium steady states. NESS correspond to the long-time limit under constant driving by an external reservoir Zia and Schmittmann (2007). As such, time symmetry is broken, but the fluxes within the system are time independent, and so are the distributions of microtrajectories.

Compared to equilibrium reweighting that focuses on the sampling of microstates, dynamical reweighting of NESS considers microtrajectories—collections of microstates—which become computationally intractable for all but the smallest of systems. To complicate things further, the length of time of a microtrajectory is a priori unknown. To resolve these issues, we map all trajectories to a first-order Markov process. This coarse-graining of microtrajectories leaves us with the easier task of sampling transition probabilities, and subsequently constructing microtrajectories out of the combination of individual microtransitions. Such an approach facilitates the sampling of microtrajectories.

Markov state models (MSM) discretize configurational space into so-called microstates (i.e., collection of microscopic states) as well as time in terms of steps of constant length $\tau$ (i.e., the lagtime), thereby mapping trajectories to a discrete-time Markov Chain Bowman et al. (2013). All observed transitions are collected to infer a transition probability matrix $p_{ij}(\tau)$, where $i$ and $j$ label microstates. An appropriate space and time discretization helps fulfill the Markovian assumption Prinz et al. (2011). Markov state models have proven powerful tools for reaching time scales that are unattainable by brute-force computer simulations Plattner et al. (2017).

Utilizing the Markovian assumption, the microtrajectories of the abovementioned cross-entropy functional (Eqn. 1) reduces to

[TABLE]

where $\mathbf{\pi}$ corresponds to the stationary distribution Lee and Pressé (2012). In the absence of constraints, the minimum of the cross entropy is its prior $q_{ij}$. Previous work has shown how to constrain the system according to microscopic and/or macroscopic constraints: ($i$) the matching of simulation and experimental data at equilibrium by enforcing detailed balance Rudzinski et al. (2016); Dixit and Dill (2018); ($ii$) inferring kinetic rates given variations in equilibrium populations Wan et al. (2016); or ($iii$) by using the stationary distribution and macroscopic constraints corresponding to a NESS experiment Dixit et al. (2015). Such macroscopic constraints are typically process dependent and not always known. In the current Letter, we propose to constrain dynamics to NESS by drawing microscopic balance constraints from stochastic thermodynamics and enforce them on the cross entropy (Eqn. 4). We consider several constraints: The conservation of probability flow through so-called global balance, $\sum_{i}\pi_{j}p_{ji}=\sum_{i}\pi_{i}p_{ij}$, allowing for global fluxes in the system; Normalization considerations imply $\sum_{j}p_{ij}=1$ and $\sum_{i}\pi_{i}=1$. Since all transition probabilities in NESSs are time-independent the existence of a steady-state distribution $\mathbf{\pi}$ is guaranteed Evans et al. (2010, 2016).

Furthermore, we want to constrain the system according to microscopic reversibility Crooks (2011, 1998)

[TABLE]

where $\mathcal{P}[\Gamma(+t)|f(+t)]$ denotes the probability of observing the time-forward trajectory $x(+t)$ under the external driving force $f(+t)$, $\mathcal{P}[\bar{\Gamma}(-t)|\bar{f}(-t)]$ points at the time-reversed trajectory, while $Q[\Gamma(+t)|f(+t)]$ refers to the amount of heat exchanged between the system and the reservoir along a given trajectory and acting forces. Critically, this links the probability of a forward trajectory with its time-reversed counterpart. In case of equilibrium dynamics, the relation becomes path independent and simplifies to detailed balance Zhang et al. (2012). For a more general expression, we integrate Eqn. 5 over the complete set of initial states $i$, target states $j$, as well as the set of all trajectories connecting them, to obtain the coarser expression (see derivation in S1)

[TABLE]

where $\Delta S_{ij}$ is called local entropy production, describing the amount of work an external reservoir has to perform on the system to transition between two states. This quantity naturally generalizes detailed balance Maes et al. (2007) and will be used as a microscopic constraint on the Caliber.

Having defined all constraints, the Caliber functional becomes

[TABLE]

where $\mu_{i}$, $\zeta$, $\nu_{j}$ and $\alpha_{ij}$ are Lagrange multipliers. Eqn. 7 modifies the equilibrium-reweighting Caliber (Eqn. 2) in several ways: ($i$) The entropy expression is discretized and replaced by the path entropy; ($ii$) transition probabilities are normalized; ($iii$) a global-balance condition ensures a steady state; and ($iv$) local entropy production is introduced as a NESS extension to detailed balance. The solution is obtained by minimizing with respect to the set of transition probabilities and the stationary distribution. Assuming that $\Delta S_{ij}$ is small (see derivation in S2), we obtain

[TABLE]

which only depends on $\Delta S^{q}_{ij}=\ln\left(q_{ij}/q_{ji}\right)$ from the reference simulation, $\Delta S_{ij}$, and the unknown constants $c_{i}$. We note that $\Delta S_{ij}$ corresponds to the local entropy production in the target state. Analogous to histogram reweighting, the unknown coefficients $c_{i}$ can be found via nonlinear relationships Bereau and Swendsen (2009)

[TABLE]

The set of equations is convex (see S4) and is solved by self iteration from a randomly-selected starting point until convergence.

An alternative formalism by means of the Girsanov theorem was introduced, which relies on single-trajectory probability reweighting and the Boltzmann distribution to estimate the change between equilibrium state points Donati et al. (2017). In another work, equilibrium transition rates are reweighted by a Maximum Caliber formalism enforcing the Boltzmann distribution Wan et al. (2016). In contrast, the present method reweights MSMs in NESS without prior knowledge of the steady-state distribution, but rather through the entropy production.

Application. The reweighting procedure is tested on a single particle driven by a non-conservative force $f$ along a periodic one-dimensional potential $U(x)$ (see Fig. 1). The non-conservative force may emerge from magnetic fields, mechanical flows, or mechanical dragging. An analogous setup was experimentally studied, using silica spheres on a tilted surface Ma et al. (2017). The overdamped equation of motion for the particle is given by

[TABLE]

where $T$ is the temperature of a canonical reservoir coupled to the system by friction constant $\gamma$. $R(t)$ is a $\delta$-correlated Gaussian process with mean [math]. Both the temperature and the potential energy are fixed, while the reweighting is performed over different ranges of non-conservative forces. We report results in reduced units, where the box size is set to $\mathcal{L}$, the mass of the particle is set to $\mathcal{M}$, and energy is measured in $\epsilon$. The temperature is chosen to be $T=1~{}\epsilon/k_{\rm B}$, the energy barriers shown in Fig. 1 are $2-4~{}k_{\rm B}T$ . Following our Markovian approximation, the model is separated in 60 microstates of equal size and a lagtime $\tau=8\cdot 10^{-4}\mathcal{T}$, where $\mathcal{T}=\mathcal{L}\sqrt{\mathcal{M}/{\epsilon}}$ is the unit of time. The integration time step was set to $\delta t=10^{-5}~{}\mathcal{T}$. The non-conservative force is varied between 0 and $9~{}\epsilon/\mathcal{L}$.

For the given model, we want to derive an analytic expression for the local entropy productions $\Delta S_{ij}$, required by the reweighting procedure (Eqn. 8). Given an underlying force field $\mathbf{F}=-\nabla U(x)+\mathbf{f}$, the local entropy production of a single continuous trajectory $\mathbf{\Gamma}(t)$ is given by

[TABLE]

where the quantity is integrated over time, $\mathbf{\dot{\Gamma}}$ is the velocity and $T$ is the temperature Seifert (2005). Assuming a constant non-conservative force $\mathbf{f}$ and making use of a numerically discretized trajectory, $\mathbf{\Gamma}(t)\approx\{\mathbf{x}_{k}\}$, we approximate $\Delta S$ between starting and target points $x_{i}$ and $x_{j}$, respectively,

[TABLE]

Because the entropy production of forward and backward steps directly cancel, the quantity is unaffected by path variations in one dimension. Still, the periodic boundary conditions permit two different results between $i$ and $j$, as indicated in Fig. 1: the shorter and longer paths (green and orange, respectively), such that Eqn. 12 has two solutions. By choosing the lagtime of the MSM reasonably short, we effectively scale down the longer paths to a negligible weight. As such, our expression for the local entropy production becomes

[TABLE]

We find excellent agreement between this expression and Eqn. 6 when directly sampled from an MSM: a weighted average error of 1%, which does not affect the quality of the reweighting upon insertion in Eqn. 8. A detailed comparison is illustrated in S3.

To assess our reweighting procedure, we monitor both static and kinetic properties: ($i$) the stationary distribution of the particle position and ($ii$) the first-passage-time distributions between metastable states. The metastable states are labeled $A$, $B$, and $C$ (Fig. 1).

Figure 1 shows the stationary distributions of the particle position both in equilibrium and under the influence of a driving force. We reweight the simulation data from equilibrium to the NESS and vice versa, demonstrating that the correct static distributions are recovered when reweighting both in and out of equilibrium—a result that holds for any pair of state points as described further below.

Turning to dynamical properties, Fig. 2 reports the first-passage-time distributions between two metastable states in equilibrium and under a strong driving force ($f=9~{}\epsilon/\mathcal{L})$. The change in the broadness of the distributions (Fig. 2a) and the shift in the peak position (Fig. 2b) suggest that the set of dominant trajectories change significantly under driving. This example shows that the external driving force changes both the timescale and corresponding processes of a transition. Here again, the reweighting procedure recovers the first-passage-time distributions in either directions: from equilibrium to NESS and vice versa.

This analysis is extended to other state points by comparing the mean, variance, and skewness of the first-passage-time distributions in Fig. 3. The equilibrium system is chosen as a reference and is continuously reweighted to off-equilibrium driven systems, even though any other reference state point could be selected. Additionally, the reweighting procedure is applied to equilibrium systems under variation of potential (see S6 for results). The reference potential energy surface (Fig. 1) is such that several pairs of processes show the same dynamics at equilibrium: The transition $A\rightarrow B$ and $C\rightarrow B$, but also $B\rightarrow A$ and $B\rightarrow C$ as well as $A\rightarrow C$ and $C\rightarrow A$.

Upon driving the system off-equilibrium, these symmetries break—a phenomenon captured by the reweighting procedure. Increasing $f$ strongly affects the mean-first-passage times (Fig. 3), thereby altering the nature of the slowest processes. We first analyze the metastable transitions situated along the direction of $f$. While the processes $A\rightarrow B$ and $B\rightarrow C$ speed up with increasing driving, $C\rightarrow A$ shows a non-monotonic behavior: it first slows down up to a driving force $f\approx 4~{}\epsilon/\mathcal{L}$ before speeding up. The variance reveals that a broader selection of trajectories becomes dominant before the threshold. Turning to metastable transitions that oppose to the driving force, only the process $C\rightarrow B$ slows down with driving, while $A\rightarrow C$ constantly speeds up, despite unfavorable driving. Process $B\rightarrow A$ shows non-monotonic behavior, similar to $C\rightarrow A$. This counterintuitive behavior can be explained by the growing number of long transitions from $B$ to $A$ via $C$.

Reweighting implies the existence of an invariant quantity, irrespective of the state point or driving force. Here we can isolate the following invariant $I_{ij}=\sqrt{q_{ij}q_{ji}}\exp\left(\frac{c_{i}+c_{j}}{4}+\zeta\right)$ (see S5 for derivation). We can thereby rewrite the abovementioned solution of the Caliber (Eqn. 8)) as

[TABLE]

using the normalization $Z_{ij}=\exp{\left(\frac{-c_{i}-c_{j}}{4}\right)}$. Note that $Z_{ij}$ depends on both $I_{ij}$ and $\Delta S_{ij}$ via the relation $1=\sum_{j}\frac{I_{ij}}{Z_{ij}}\exp{\left(\frac{\Delta S_{ij}}{2}\right)}$ and accounts for the interconnection of the states. We draw similarities with equilibrium reweighting in Eqn. 3: ($i$) The probability is proportional to the product of an invariant and an exponential function (the density of states and the Boltzmann factor in equilibrium); ($ii$) The partition function depends on the control variable ($T$ or $\Delta S_{ij}$); and ($iii$) The reweighting only depends on relative quantities, only requiring knowledge of temperature difference or changes in the local entropy production. Both procedures show striking similarities in their derivation, functional form, and properties.

The present reweighting method is a generalization of existing Likelihood and Maximum Caliber methods that have been applied to systems in and out of equilibrium with varying microscopic and macroscopic constraints. The microscopic expression for the local entropy production acts as a local constraint that generalizes detailed balance for NESS. We show that this choice governs static and dynamic properties of a NESS and enables us to reproduce these properties over a wide range of driving. The analytic expression for the relative entropy production allows us to continuously tune the external driving force and quantitatively reweight the stationary distribution and kinetic properties.

The Maximum Caliber formalism in combination with local entropy productions offer an analytic relation between NESSs. Dynamical data of a system can be gathered in a driving-invariant quantity and detailed kinetic information at any thermodynamic state point can be recovered. This idea allows to populate rare transition paths Bolhuis et al. (2002): A driving force may push the system to discover new paths, the reweighting procedure recovers detailed dynamical information of the sampled path at any another thermodynamic state point. In case equilibrium dynamics are of interest, the entropy productions only depend on the free energy. Low weight trajectories can thus be calculated with high accuracy and no further information. By tuning the relative local entropy productions of the system, the reweighting allows to study dynamical properties and pathways in NESS without further simulation.

Acknowledgments

We thank Joseph F. Rudzinski, Dominik Spiller and Johannes Zierenberg for insightful discussions. This work was supported in part by the Emmy Noether program of the Deutsche Forschungsgemeinschaft (DFG) to TB, the Graduate School of Excellence Materials Science in Mainz (MAINZ) to MB, and the European Research Council under the European Union’s Seventh Framework Programme (MOLPROCOMP) to KK.