Non-equilibrium currents and fluxes

Non-equilibrium thermodynamics

From a physics perspective, the most obvious property of living systems is that they are not in thermodynamic equilibrium. Non-equilibrium systems are open systems, in contact with an environment from which energy is extracted and transformed, dissipating heat and producing entropy along the way. Hence non-equilibrium processes are irreversible and characterized by currents or fluxes, which persist even when the system is in a stationary state.

Great examples of non-equilibrium currents in living systems are the fluxes in metabolic pathways, linked series of chemical reactions where the product of one reaction acts as the substrate of the next.

Definition of entropy production

A good review on irreversibility and entropy production is Jarzynski’s article. Here I’ll keep things very simple.

Imagine watching a movie of a dynamic process and then watching the same movie backwards in time. If you cannot tell whether you’re watching the forward or backward movie, then the process is called reversible. If you can tell the difference, that is, can guess the direction of time when watching either movie, the process is irreversible. An irreversible process produces entropy and the amount of entropy produced per unit of time is a measure for the degree of irreversibility.

Consider a system whose dynamics can be modelled by a Markov process. For a sequence of times \(t_0=0<t_1<t\dots<t_N=T\), a trajectory is a sequence of states \(x_{[0,T]}=(x(0),x(t_1),\dots,x(T))\). The time-reversal of a trajectory is defined as the sequence of states \(\tilde{x}_{[0,T]}=(x(T),x(t_{N-1}),\dots,x(0))\)¹, that is, \(\tilde{x}(t_i)=x(T-t_i)\). In a non-equilibrium stationary state, each trajectory has a probability (density) of being observed, described by the forward path-wise distribution \(\mathbf{P}_{[0:T]}\). Similarly, the time-reversed distribution \(\tilde{\mathbf{P}}_{[0:T]}\) describes for a trajectory \(x_{[0,T]}\) the probability (density) of observing the time-reversed trajectory \(\tilde{x}_{[0,T]}\).

Mathematically, the entropy production rate is quantified by the relative entropy or Kulback-Leibler divergence between the forward and time-reversed path distributions:²

\[ e_P = \frac{1}{T}D_{KL}\bigl(\mathbf{P}_{[0,T]} \mid \tilde{\mathbf{P}}_{[0,T]}\bigr) \]

Non-equilibrium current in the multivariate Ornstein-Uhlenbeck process

The simplest mathematical model with a non-equilibrium current is a multi-variate Ornstein-Uhlenbeck process. The model is simple, because its properties can be expressed in terms of normal distributions. It is also in some sense generic, because it arises in the so-called linear noise approximation (LNA) as an approximate description of the fluctuations in more complex stochastic processes, incl. chemical reaction systems. The LNA can be considered as an extension of the central limit theorem for stochastic processes.

The definitive guide to the multi-variate Ornstein-Uhlenbeck process, and its non-equilibrium properties in particular, is the paper “Characterising the nonequilibrium stationary states of Ornstein-Uhlenbeck processes”. To make things easy, I’ll follow their notation, except for not bold-facing synbols for vectors and matrices. The model is defined by a set of coupled linear Langevin equations, in matrix-vector notation

\[ \frac{dx}{dt} = - B x(t) + \eta(t) \]

where \(\eta(t)\) is a vector of Gaussian white noises, such that

\[ \langle \eta(t)\eta^T(t')\rangle = 2 D \delta(t-t'). \]

\(B\) is called the drift matrix, and \(D\) the diffusion matrix

Starting from an initial position \(x(0)\), as \(t\to\infty\), the system reaches a Gaussian stationary state

\[ P_{st}(x) = \mathcal{N}(0,S) \]

with the covariance matrix \(S\) solving the continuous-time Lyapunov equation

\[ BS + SB^T = 2D \]

Using a linear transformation it is always possible to reparametrise the problem such that \(D=I\), the identity matrix, and we will hence assume that this is the case.

The OUP is reversible and satisfies detailed balance if

\[ BS = SB^T = D \]

Hence we can define

\[ B_{rev} = DS^{-1}=S^{-1} \]

where we used the assumption \(D=I\). For any antisymmetric matrix \(Q=-Q^T\), we can define

\[ B_{irr} = QS^{-1} \]

It follows immediately that \(S\) solves the Lyapunov equation for any \(B=B_{rev} + B_{irr}\). In other words, we can consider a family of OUPs which all have the same stationary state covariance \(S\), but varying degrees of irreversibility.

To justify this terminology, we consider the transition probabilities \(T(x,t\mid x_0,t_0)\) for the process to be at position \(x\) at time \(t\) given it was at position \(x_0\) at time \(t_0<t\). These transition probabilities solve a Fokker-Planck equation, which can be written in the form

\[ \frac{\partial P(x,t)}{\partial t} = \nabla \cdot J(x) \]

with initial condition \(P(x,t_0)=\delta(x-x_0)\), where \(J(x)\) is the probability current

\[ J(x) = -B_{irr} x P_{st}(x) \]

with \(B_{irr}\) as defined above and \(P_{st}\) the stationary distribution. To make the connection complete, it can be shown that the entropy production of the OUP is equal to

\[ e_P = \int \frac{J^T(x) D^{-1} J(x)}{P_{st}(x)} dx \]

In other words, both the entropy production and the probability current vanish if \(B_{irr}=0\), that is, if the process is reversible.

Although the maths is fairly straightforward, I struggled for a long time to have a clear mental image of what the probability current and entropy production in this process really mean. Enlightenment came after coming across the paper “The entropy production of stationary diffusions”, and especially their Figure 2. Since the paper is otherwise quite advanced, I’ll repeat their beautiful and simple simulations here.

We consider a 2-dimensional OUP, and make the further assumption that \(S\) is diagonal such that with two parameters \(\sigma^2\) and \(\theta\) we can define

\[ \begin{aligned} S &= \begin{pmatrix} \sigma^2 & 0\\ 0 & \sigma^2 \end{pmatrix}\\ B_{rev} = S^{-1} &= \begin{pmatrix} \frac{1}{\sigma^2} & 0\\ 0 & \frac{1}{\sigma^2} \end{pmatrix}\\ Q &= \begin{pmatrix} 0 & \theta\\ -\theta & 0 \end{pmatrix}\\ B_{irr} = QS^{-1} &= \begin{pmatrix} 0 & \frac{\theta}{\sigma^2} \\ -\frac{\theta}{\sigma^2} & 0 \end{pmatrix} \end{aligned} \]

To specify the OUP, we follow the DifferentialEquations.jl documentation and define drift and diffusion functions, where the two parameters are collected in a parameter vector \(p=(\sigma^2,\theta)\):

using StochasticDiffEq
function f!(du,u,p,t)
    du[1] = -(u[1] + p[2]*u[2])/p[1]
    du[2] = -(-p[2]*u[1] + u[2])/p[1]
end

function g!(du,u,p,t)
    du[1] = √2
    du[2] = √2
end

Note that the first function indeed computes \(-Bu\), while the second one indeed defines a vector of Gaussian white noises with covariance matrix \(2D=2I\).

We set parameter values, sample an initial position from the stationary distribution, and sample trajectories, one for the reversible system (\(\theta=0\)) and one for an irreversible system:

σ² = 0.25
θ = 5.0

u0 = √σ²*randn(2,1) # sample initial position from the stationary distribution
tspan = (0.0, 10.0)

prob_eq = SDEProblem(f!, g!, u0, tspan, (σ²,0.0))
sol_eq = solve(prob_eq, EM(), dt=0.01)

prob_noneq = SDEProblem(f!, g!, u0, tspan, (σ²,θ))
sol_noneq = solve(prob_noneq, EM(), dt=0.01)

Here is the all telling visualization of the trajectories:

using Plots
lo = @layout([a b;c d])
@gif for i in 2:length(sol_eq.u)
    p1 = plot(sol_eq[1:i], idxs = (1, 2),xlim=(-2,2),ylim=(-2,2), legend=false)
    scatter!([sol_eq.u[i][1]], [sol_eq.u[i][2]])
   
    p2 = plot(sol_noneq[1:i], idxs = (1, 2),xlim=(-2,2),ylim=(-2,2), legend=false)
    scatter!([sol_noneq.u[i][1]], [sol_noneq.u[i][2]])

    p3 = plot(sol_eq[1:i], legend=false, xlim=(tspan[1], tspan[2]),ylim=(-2,2))

    p4 = plot(sol_noneq[1:i], legend=false, xlim=(tspan[1], tspan[2]),ylim=(-2,2))

    plot(p1, p2, p3, p4; layout=lo)
end

[ Info: Saved animation to /tmp/jl_ngD0f6xaZY.gif

Now it is clear that the probability current is a circular motion in phase space (counterclockwise if \(\theta>0\))! The equilibrium process (\(\theta = 0\)) is reversible, because if we play the trajectory forwards or backwards in time, it would be impossible to tell which is which. For the irreversible process (\(\theta \neq 0\)) this is clearly not the case, because the direction of rotation informs immediately on the direction of time.

With hindsight, it is also easy to see that the probability current corresponds to a motion along the contours of the stationary distribution. Indeed consider a one-parameter flow \(x(s)\) that has \(J(x)\) as its tangent vector:

\[ \frac{dx(s)}{ds} = J(x(s)) \]

Then

\[ \begin{aligned} \frac{d}{ds} P_{st}(x(s)) &= \left(\nabla P_{st}(x(s))\right) \cdot \frac{dx(s)}{ds}\\ &= P_{st}(x(s))^2 \bigl(S^{-1}x(s))\cdot \bigl(Q S^{-1}x(s))\\ &= 0 \end{aligned} \]

Where the last equality frollows from the fact that \(Q\) is asymmetric, such that for any vector \(v\)

\[ v\cdot (Qv) = v^TQv = v^TQ^Tv = - v^TQv= - v\cdot (Qv) \]

that is, \(v\cdot (Qv)=0\).

In other words, along the flow \(x(s)\), \(P_{st}(x(s))\) is constant, that is, the current flows along the contours of \(P_{st}\).

Reconciling the model and thermodynamics

We have seemingly arrived at a contradiction: the thermodynamic description of irreversible processes and entropy production emphasizes that non-equilibrium systems must be open systems that exchange particles or energy with the environment, while the OUP process describes a system in two (or more) variables that is apparently closed, that is, there is no description of (exchange with) an environment in the model, and yet produces entropy. The resolution of this apparent contradiction lies in the asymmetry of the reaction rates (the non-zero parameter \(\theta\) in the two-dimensional model). In a closed physical system that is not in contact with an environment, reaction rates will always balance after some time, that is, the system will reach detailed balance. This is the very definition of reaching equilibrium. In other words, to maintain asymmetric and constant reaction rates, the system must be in contact with an environment that supplies the necessary energy. Hence the assumption that \(\theta\neq 0\) and constant over time implicitly means that we are modelling an open system. The connection can be made explicit via the principle of local detailed balance, which (roughly) states that the asymmetry in reaction or transition rates is balanced by the change in entropy of the mediating reservoir. A simple example of local detailed detailed balance can be found in Box 1 of this paper, and a more complete description here.

Footnotes

1. I’m assuming for simplicity that \(x\) does not contain any momentum variables that change sign under time-reversal and glossing over a lot of other details too.

2. I’m following roughly the definition in this paper.