Intrinsic and extrinsic thermodynamics for stochastic population processes with multi-level large-deviation structure

Eq. (4) is decomposed in terms of the quantities in equations (5–7) as

The total is non-negative and the third line is independently non-negative, as already mentioned. The sum of the first two lines is also non-negative, a result that can be proved as a fluctuation theorem for what is normally called “total entropy change” Esposito:fluct_theorems:10 .¹⁰¹⁰10More detailed proofs for a decomposition of the same sum will be given below.

Here we encounter the first property that makes a decomposition of a thermal system “natural”. The term in ${\dot{D}}^{e\mid s}_{{{\rm n}}_{s}}\!\left({\rho}^{e\mid s}\parallel{\underline{\rho}}^{e\mid s}\right)$ is not generally considered, and it does not need to be considered, because thermal relaxation at the microscale makes transitions in ${{\rm n}}_{e}$ at fixed ${{\rm n}}_{s}$ conditionally independent of transitions that change ${{\rm n}}_{s}$. Total relative entropy changes as a sum of two independently non-decreasing contributions.

The first and second lines in Eq. (9) are not likewise independently non-negative. Negative values of the second or first line, respectively, describe phenomena such as randomization-driven endothermic reactions, or heat-driven information generators. To the extent that they do not use thermalization in the microscale to make the system and environment conditionally independent, as the third term in Eq. (9) is independent, we say they do not provide a natural system/environment decomposition.

Remarkably, a natural division does exist, based on the housekeeping heat introduced by Hatano and Sasa Hatano:NESS_Langevin:01 . The decomposition uses the solution ${\bar{\rho}}^{s}$ to ${\mathbb{T}}^{s}{\bar{\rho}}^{s}=0$, which would be the stationary marginal distribution for the system $s$ at the instantaneous value of ${\rho}^{e\mid s}$. As for the whole-system stationary distribution $\underline{\rho}$, we restrict to cases in which ${\bar{\rho}}^{s}$ exists and is normalizable.¹¹¹¹11This can be a significant further restriction when we wish to study limits of sequences of environments to model chemostats. For systems with unbounded state spaces, such as arise in polymerization models, it is quite natural for a chemostat-driven system to possess no normalizable steady-state distributions. In such cases other methods of analysis must be used Esposito:copol_eff:10 .

Treating ${\bar{\rho}}^{s}$ as fixed and considering only the dynamics of ${\rho}^{s}$ with ${\bar{\rho}}^{s}$ as a reference, we may consider a time derivative $\dot{D}\!\left({\rho}^{s}\parallel{\bar{\rho}}^{s}\right)$ in place of Eq. (5). Note that the rates $w^{s}_{{{\rm n}}_{s}{{\rm n}}^{\prime}_{s}}$ in the transition matrix ${\mathbb{T}}^{s}$ no longer need satisfy any simplified balance condition in relation to ${\bar{\rho}}^{s}$, such as detailed balance. Non-negativity of $-\dot{D}\!\left({\rho}^{s}\parallel{\bar{\rho}}^{s}\right)$ was proved by Schnakenberg Schnakenberg:ME_graphs:76 by an argument that applies to discrete population processes with normalized stationary distributions of the kind assumed here. We will derive a slightly more detailed decomposition proving this result for the case of CRNs, in a later section.

The dissipation term that complements the change in $\dot{D}\!\left({\rho}^{s}\parallel{\bar{\rho}}^{s}\right)$ is obtained by shifting the “environmental” entropy change (8) by $\log{\bar{\rho}}^{s}$ to obtain the housekeeping heat¹²¹²12These terms are used because this is how they are known. No energy interpretation is assumed here, so a better term would be “housekeeping entropy change”.

Non-negativity of Eq. (10) is implied by a fluctuation theorem Speck:HS_heat_FT:05 , and a time-local proof with interesting further structure for CRNs will be given below.

The total change in relative entropy (4) is then the sum

Each line in Eq. (11) is now independently non-negative. The first measures a gain of entropy within the system $s$, conditionally independent of changes in the environment given the marginal transition matrix ${\mathbb{T}}^{s}$. The third measures a gain of entropy in the environment $e$ independent of any changes in the system at all. The second, housekeeping entropy rate, measures a change of entropy in the environment that is conditionally independent of changes of entropy within the system, given ${\mathbb{T}}^{s}$ as represented in ${\bar{\rho}}^{s}$. Any of the terms may be changed, holding the conditioning data ${\bar{\rho}}^{s}$ or ${{\rm n}}_{s}$, or omitted, without changing the limits for the others. They respect the conditional independence created by thermalization in the microscale, and by that criterion constitute a natural decomposition of the system.

We take ${\mathbb{T}}^{s}$ and $D\!\left({\rho}^{s}\parallel{\bar{\rho}}^{s}\right)$ as a specification of the intrinsic thermodynamics of the system $s$, analogous to the role of intrinsic curvature of a manifold in differential geometry. The vector (indexed by pairs of system indices) of housekeeping entropy differentials, $\left\{{\sigma}^{\rm HK}_{{{\rm n}}_{s}{{\rm n}}^{\prime}_{s}}\right\}$, correspondingly defines the way $s$ is thermodynamically embedded in the mesoscale system, analogous to the role of components of an embedding curvature for a submanifold within larger manifold.

The natural decomposition (11) differs from the information/heat decomposition (9) in what is regarded as inherent to the system versus the environment. The attribution of Shannon entropy as a “system” property follows from the fact that it involves only ${\rho}^{s}_{{{\rm n}}_{s}}$, and its change counts only actual transitions ${{\rm n}}_{s}^{\prime}\rightarrow{{\rm n}}_{s}$ with rate $w^{s}_{{{\rm n}}_{s}{{\rm n}}_{s}^{\prime}}$. Likewise, the “environment” heat (8) is a function only of the actual distribution ${\rho}^{e\mid s}$ and the realized currents.

Those events that occur within $s$ or $e$, however, fail to capture the additional features of ${\mathbb{T}}^{s}$: that specific transitions are coupled as a system, and that they have the dependence on ${\rho}^{e\mid s}$ of Eq. (2). The stationary distribution ${\bar{\rho}}^{s}$ is the function of ${\mathbb{T}}^{s}$ reflecting its status as a system.

The differences between the three entropy-change terms (7,8,10) are differences of the Hartley informations Hartley:information:28 , respectively for ${\underline{\rho}}^{s}$ or ${\bar{\rho}}^{s}$. In the natural decomposition (11), they are not acted upon by the time derivative, but rather define the tangent plane of zero change for $\log\rho$ terms that are acted upon by the transition matrix, and resemble connection coefficients specifying parallel transport in differential geometry.

The $P$-dimensional Laplace transform of a distribution ${\rho}_{{\rm n}}$ on discrete population states gives the moment-generating function (MGF) for the species-number counts ${\left\{{{\rm n}}_{p}\right\}}_{p\in 1,\ldots,P}$. For time-dependent problems, it is convenient to work in the formal Doi operator algebra for generating functions, in which a vector of raising operators $a^{\dagger}\equiv\left[a^{\dagger}_{p}\right]$ are the arguments of the MGF, and a conjugate vector of lowering operators $a\equiv\left[a_{p}\right]$ are the formal counterparts to $\partial/\partial a^{\dagger}$. MGFs are written as vectors $\left|\rho\right)\equiv\sum_{{\rm n}}{\rho}_{{\rm n}}{\left(a^{\dagger}\right)}^{{\rm n}}\left|0\right)\equiv\sum_{{\rm n}}{\rho}_{{\rm n}}\left|{\rm n}\right)$ in a Hilbert space built upon a ground state $\left|0\right)$, and the commutation relations of the raising and lowering operators acting in the Hilbert space are $\left[a_{p},a^{\dagger}_{q}\right]={\delta}_{pq}$.

The basis vectors corresponding to specific population states are denoted $\left|{\rm n}\right)$. They are eigenvectors of the number operators $a^{\dagger}_{p}a_{p}$ with eigenvalues ${{\rm n}}_{p}$:

Through by-now-standard constructions Smith:LDP_SEA:11 , the master equation $\dot{\rho}=\mathbb{T}\rho$ is converted to a Liouville equation for time evolution of the MGF,

in which the Liouville operator $\mathcal{L}\!\left(a^{\dagger},a\right)$ is derived from the elements of the matrix $-\mathbb{T}$.

To define exponential families and a cumulant-generating function (CGF), it is convenient to work with the Laplace transform with an argument that is a vector $z\equiv\left[z_{p}\right]$ of complex coefficients. The corresponding CGF, $-\Gamma\!\left(\log z\right)$,¹⁴¹⁴14We adopt the sign for $\Gamma$ corresponding to the free energy in thermodynamics. Other standard notations, such as $\psi$ for the CGF Amari:inf_geom:01 , are unavailable because they collide with notations used in the representation of CRNs below. for which the natural argument is $\log z$, is constructed in the Doi Hilbert space as the inner product with a variant on the Glauber norm,

$z$ will be called the tilt of the exponential family, corresponding to its usage in importance sampling Siegmund:IS_seq_tests:76 .

An important quantity will be the vector of expectations of the number operators in the tilted distribution

Time evolution of the CGF follows from Eq. (13), as

Under the saddle-point or leading-exponential approximation that defines the large-deviation limit (developed further below), the expectation of $\mathcal{L}$ in Eq. (16) is replaced by the same function at classical arguments

To establish notation and methods, consider first distributions ${\rho}_{{\rm n}}$ that are convex with an interior maximum in ${\rm n}$. Then the gradient of the CGF

gives the mean (15) in the tilted distribution.

The stochastic effective action is the Legendre transform of $-\Gamma$, defined as

For distributions over discrete states, to leading exponential order, $S_{\rm eff}$ is a continuously-indexed approximation to minus the log-probability: $S_{\rm eff}\!\left(n\right)\sim-{\left.\log{\rho}_{{\rm n}}\right|}_{{\rm n}\approx n}$. Its gradient recovers the tilt coordinate $\theta$,

and the CGF is obtained by inverse Legendre transform.

The time evolution of $S_{\rm eff}$ can be obtained by taking a total time derivative of Eq. (19) along any trajectory, and using Eq. (18) to cancel the term in $\dot{\theta}$. The partial derivative that remains, evaluated using Eq. (17), gives

Eq. (22) is of Hamilton-Jacobi form, with $-\mathcal{L}$ filling the role of the Hamiltonian. (The reason for this sign correspondence, which affects nothing in the derivation, will become clear below.)

Partial derivatives with respect to $t$ and $n$ commute, so the relation (20), used to evaluate $\partial/\partial n$ of Eq. (22), gives the relation

The dual construction for the time dependence of $n$ from Eq. (18) gives

The evolution equations (23,24) describe stationary trajectories of an extended-time Lagrange-Hamilton action functional which may be written in either coherent-state or number-potential coordinates, as¹⁶¹⁶16The same action functionals are arrived at somewhat more indirectly via 2-field functional integral constructions such as the Doi-Peliti method.

From the form of the first term in either line, it is clear that the two sets of coordinates relate to each other through a canonical transformation Goldstein:ClassMech:01 .

A special case, which is illustrative for its simplicity and which arises for an important sub-class of stochastic CRNs, is the case when $\underline{\rho}$ is a coherent state, an eigenvector of the lowering operator $a$. Coherent states are the generating functions of product-form Poisson distributions, or cross-sections through such products if the transitions in the population process satisfy conservation laws. They are known Anderson:product_dist:10 to be general solutions for CRN steady states satisfying a condition termed complex balance, and the fixed points associated with such stationary distributions are also known to be unique and interior (no zero-expectations for any ${{\rm n}}_{p}$) Feinberg:notes:79 .

Let $\underline{n}$ be the eigenvalue of the stationary coherent state: $a\left|{\rho}^{\left(\underline{n}\right)}\right)=\underline{n}\left|{\rho}^{\left(\underline{n}\right)}\right)$. Then the mean in the tilted distribution lies in a simple exponential family, $\bar{n}\equiv n_{\underline{\rho}}\!\left(\theta\right)=e^{\theta}\underline{n}$ (component-wise), and the tilted macrostate ${\rho}^{\left(\bar{n}\right)}$ is also a coherent state: $a\left|{\rho}^{\left(\bar{n}\right)}\right)=\bar{n}\left|{\rho}^{\left(\bar{n}\right)}\right)$.

The logarithm of a product-form Poisson distribution in Stirling’s approximation is given by

$D_{f}\!\left({\rm n}\parallel\bar{n}\right)$, known as the $f$-divergence, is a generalization of the Kullback-Leibler divergence to measures such as $n$ which need not have a conserved sum. The Lyapunov function from Eq. (4) reduces in the same Stirling approximation to

giving ${\underline{S}}_{\rm eff}\!\left(\bar{n}\right)$ in Eq. (29). On coherent states, the Kullback-Leibler divergence on distributions, which may be of arbitrarily large dimension, reduces to the $f$-divergence on their extensive state variables which have dimension $P$.

The coherent states play a much more general role than their role as exact solutions for the restricted case of complex-balanced CRNs. In the Doi-Peliti 2-field functional integral formalism Doi:SecQuant:76 ; Doi:RDQFT:76 ; Peliti:PIBD:85 ; Peliti:AAZero:86 for generating functionals over discrete-state stochastic processes, the coherent states form an over-complete basis in the Peliti representation of unity. The saddle-point approximation on trajectories, which yields the classical actions (25) and the resulting Hamilton-Jacobi equations, approximates expectations in the exact distribution by those in the nearest coherent-state basis element. Observables in macrostates are thus mapped to observables in coherent states, although in cases when the coherent state is not an exact solution, the saddle-point condition may be sensitive to which observable is being evaluated.

Systems with multiple metastable fixed points correspond to non-convex $\underline{\rho}$ and thus multiple modes. For these, monotone decrease of $\dot{D}\!\left(\rho\parallel\underline{\rho}\right)$ in Eq. (4) does not entail monotonicity of the $f$-divergence in Eq. (31). In such systems, first passages between basins of attraction are solutions to the Hamiltonian equations (23,24) with momentum coordinate $\theta\neq 0$. Along these ${\underline{S}}_{\rm eff}\!\left(\bar{n}\right)$ increases, and that increase is what is sometimes termed the “violation of the 2nd law”.

For unimodal $\underline{\rho}$, the $\theta\neq 0$ large-deviation trajectories have a separate use and interpretation from the relaxation trajectories at $\theta=0$ that give the classical 2nd law in Eq. (49). For multi-modal $\underline{\rho}$ a special sub-class of $\theta\neq 0$ trajectories, those known as instantons and responsible for first-passages between fixed-points Cardy:Instantons:78 ; Coleman:AoS:85 , must be used to refine the interpretation of classical relaxation trajectories. That refinement relates the transient increases in the large-deviation function to the deterministic 2nd law (4) that continues to apply.

This section briefly introduces the Legendre duality that defines first-passage probabilities in metastable systems, arriving at the chain rule for entropy that separates the roles of classical and instanton trajectories. Let $\underline{n}$ be a fixed point of the Hamiltonian equations for $\underline{\rho}$, and denote by $n_{\underline{\rho}}\!\left(\theta\mid\underline{n}\right)$ the values of classical state variables $\bar{n}$ obtained along $\theta\neq 0$ trajectories from $\underline{n}$. Call these escape trajectories. The set of all $\bar{n}$ is partitioned among basins of repulsion from fixed points. Saddle points and escape separatrices are limit points of escapes from two or more basins.

Within one such basin, we may construct ${\underline{S}}_{\rm eff}$ as a Legendre transform of a summand ${\Gamma}_{\underline{n}}$ in the overall CGF, as

$\theta$ ranges only over the values that arise on escape trajectories from $\underline{n}$, which generally are bounded Smith:CRN_CTM:20 , and within that range

Next, let the function $\underline{n}\!\left(\bar{n}\right)$¹⁷¹⁷17Note that the maps $\underline{n}\!\left(\bar{n}\right)$ and $n_{\underline{\rho}}\!\left(\theta\mid\underline{n}\right)$ need not be reflexive. That is, we may have $n_{\underline{\rho}}\!\left(\theta\mid\underline{n}\!\left(\bar{n}\right)\right)\neq\bar{n}$ for any $\theta$, because escape and relaxation separatrices may differ. denote the fixed point to which a trajectory with $\theta\equiv 0$ relaxes, starting from $\bar{n}$. From the large-deviation identification of $S_{\rm eff}\!\left(n\right)\sim-{\left.\log{\rho}_{{\rm n}}\right|}_{{\rm n}\approx n}$, recognize that the large-deviation function for the stationary distribution is given by ${\underline{S}}_{\rm eff}\!\left(\underline{n}\right)=D\!\left(1_{\underline{n}}\parallel\underline{\rho}\right)$, where $1_{\underline{n}}$ denotes the special macrostate for which the stationary value is a fixed-point value $\underline{n}$.

Then the expression (29) for the Kullback-Leibler divergence appearing in Eq. (4) may be approximated to leading exponential order as

Eq. (34) is the chain rule for relative entropy. ${\underline{S}}_{\rm eff}\!\left(\underline{n}\!\left(\bar{n}\right)\right)-{\underline{S}}_{\rm eff}\!\left(\bar{n}\right)$ is the conditional entropy of the macrostate $\bar{n}$ relative to the macrostate $1_{\underline{n}\left(\bar{n}\right)}$ to which it is connected by a $\theta=0$ Hamiltonian trajectory. $-D\!\left(1_{\underline{n}\left(\bar{n}\right)}\parallel\underline{\rho}\right)$ is the unconditioned entropy of $1_{\underline{n}\left(\bar{n}\right)}$ relative to $\underline{\rho}$.

Relaxation along $\theta=0$ trajectories describes a classical 2nd law only for the conditional part of the relative entropy. Deterministic relaxation of the unconditioned entropy is derived from the refinement of the (as it turns out, only apparent) classical trajectory with an instanton sum. The general method is described in depth in Cardy:Instantons:78 and Coleman:AoS:85  Ch.7, using path integrals that require too much digression to fit within the scope of this paper. The structure of the instanton sum, and in particular the way it creates a new elementary stochastic process at the macroscale for which $D\!\left(1_{\underline{n}\left(\bar{n}\right)}\parallel\underline{\rho}\right)$ is the Lyapunov function, will be explained in Sec. IV.2.2 following Smith:LDP_SEA:11 , after the behavior of ${\underline{S}}_{\rm eff}$ along $\theta=0$ and $\theta\neq 0$ trajectories has been characterized.