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

Eric Smith  Department of Biology, Georgia Institute of Technology, 310 Ferst
Drive NW, Atlanta, GA 30332, USA  Earth-Life Science Institute, Tokyo
Institute of Technology, 2-12-1-IE-1 Ookayama, Meguro-ku, Tokyo 152-8550,
Japan  Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA  Ronin
Institute, 127 Haddon Place, Montclair, NJ 07043, USA

###### Abstract

A set of core features is set forth as the essence of a thermodynamic
description, which derive from large-deviation properties in systems with
hierarchies of timescales, but which are _not_ dependent upon conservation
laws or microscopic reversibility in the substrate hosting the process. The
most fundamental elements are the concept of a macrostate in relation to the
large-deviation entropy, and the decomposition of contributions to
irreversibility among interacting subsystems, which is the origin of the
dependence on a concept of heat in both classical and stochastic
thermodynamics. A natural decomposition is shown to exist, into a relative
entropy and a housekeeping entropy rate, which define respectively the
intensive thermodynamics of a system and an extensive thermodynamic vector
embedding the system in its context. Both intensive and extensive components
are functions of Hartley information of the momentary system stationary state,
which is information _about_ the joint effect of system processes on its
contribution to irreversibility. Results are derived for stochastic Chemical
Reaction Networks, including a Legendre duality for the housekeeping entropy
rate to thermodynamically characterize fully-irreversible processes on an
equal footing with those at the opposite limit of detailed-balance. The work
is meant to encourage development of inherent thermodynamic descriptions for
rule-based systems and the living state, which are not conceived as reductive
explanations to heat flows.

##  I Introduction

The statistical derivations underlying most thermodynamic phenomena are
understood to be widely applicable, and are mostly developed in general terms.
Yet where thermodynamics is offered as an ontology to understand new patterns
and causes in nature – the thermodynamics of computation Szilard:MD:29 ;
Landauer:IHGCP:61 ; Bennett:TC:82 ; Wolpert:stoch_thermo_comp:19  or
stochastic thermodynamics Seifert:stoch_thermo_rev:12 ,111 Stochastic
thermodynamics is the modern realization of a program to create a non-
equilibrium thermodynamics that began with Onsager Onsager:RRIP1:31 ;
Onsager:RRIP2:31  and took much of its modern form under Prigogine and
coworkers Glansdorff:structure:71 ; Prigogine:MT:98 . Since the beginning its
core method has been to derive rules or constraints for non-stationary thermal
dynamics from dissipation of free energies defined by Gibbs equilibria, from
any combination of thermal baths or asymptotic reservoirs.
The parallel and contemporaneous development of thermodynamics of computation
can be viewed as a quasistatic analysis with discrete changes in the boundary
conditions on a thermal bath corresponding to logical events in algorithms.
Early stochastic thermodynamics combined elements of both traditions, with the
quantities termed “non-equilibrium entropies” corresponding to the information
entropies over computer states, distinct from quasistatic entropies associated
with heat in a locally-equilibrium environment, with boundary conditions
altered by the explicitly-modeled stochastic state transitions.
More recently stochastic thermodynamics incorporated time-reversal methods
originally developed for measures in dynamical systems Searles:fluct_thm:99 ;
Evans:fluct_thm:02 ; Gallavotti:dyn_ens_NESM:95 ; Gallavotti:dyn_ens_SS:95
leading to a variety of fluctuation theorems Jarzynski:fluctuations:08 ;
Chetrite:fluct_diff:08 ; Esposito:fluct_theorems:10  and nonequilibrium work
relations Crooks:NE_work_relns:99 ; Crooks:path_ens_aves:00 ; Kurchan:NEWRs:07
, still however relating path probability ratios either to dissipation of heat
or to differences in equilibrium free energies. or where these methods are
taken to define a foundation for the statistical physics of reproduction or
adaptation England:statphys_selfrepl:13 ; Perunov:adaptation:15  – these
problems are framed in terms of two properties particular to the domain of
mechanics: the conservation of energy and microscopic reversibility. Other
applications of the same mathematics, with designations such as intensity of
choice Luce:choice:59 ; McFadden:quantal_choice:76  tipping points, or early
warning signs Lenton:climate_tip_ew:12 , are recognized as analogies to
thermodynamics, on the understanding that they only become the “thermodynamics
of” something when they derive its causes from energy conservation connected
to the entropies of heat.

This accepted detachment of mathematics from phenomenology, with thermodynamic
phenomena interpreted in their historical terms and mathematics kept
interpretation-free, contrasts with the way statistical mechanics was allowed
to expand the ontological categories of physics at the end of the last
century. The kinetic theory of heat Joule:sci_papers:44 ;
Clausius:mech_the_heat:65 ; Boltzmann:second_law:86  was not established as an
analogy to gambling222The large-deviation rate function underlies the solution
of the gambler’s ruin problem of Pascal, Fermat, and Bernoulli
Hald:prob_pre_1750:90 . used to describe patterns in the fluid
caloric.333However, as late as 1957 Jaynes Jaynes:ITSM_I:57 ;
Jaynes:ITSM_II:57  needed to assert that the information entropy of Shannon
referred to the same quantity as the physical entropy of Boltzmann, and was
not merely the identical mathematical function. Even the adoption of “entropy
production” to refer to changes in the state-variable entropy by irreversible
transformations – along the course of which the state-variable entropy is not
even defined – is a retreat to a substance syntax; I would prefer the less
euphonious but categorically better expression “loss of large-deviation
accessibility”. Thermodynamics instead took the experienced phenomena
involving heat, and where there had formerly been only names and metaphors to
refer to them, it brought into existence concepts capturing their essential
nature, not dissolving their reality as phenomena Rota:lect_notes:08 , but
endowing it with a semantics. The distinction is between formalization in the
service of reduction to remain within an existing ontology, and formalization
as the foundation for discovery of new conceptual primitives. Through
generalizations and extensions that could not have been imagined in the late
19th century GellMann:RG:54 ; Wilson:RG:74 ; Weinberg:phenom_Lagr:79 ;
Polchinski:RGEL:84  (revisited in Sec. VII), essentially thermodynamic
insights went on to do the same for our fundamental theory of objects and
interactions, the nature of the vacuum and the hierarchy of matter, and the
presence of stable macro-worlds at all.

This paper is written with the view that the essence of a thermodynamic
description is not found in its connection to conservation laws, microscopic
reversibility, or the equilibrium state relations they entail, despite the
central role those play in the fields mentioned Landauer:IHGCP:61 ;
England:statphys_selfrepl:13 ; Perunov:adaptation:15 ;
Seifert:stoch_thermo_rev:12 . At the same time, it grants an argument that has
been maintained across a half-century of enormous growth in both statistical
methods and applications Fermi:TD:56 ; Bertini:macro_NEThermo:09 : that
thermodynamics should not be conflated with its statistical methods. The focus
will therefore be on the patterns and relations that make a phenomenon
essentially thermodynamic, which statistical mechanics made it possible to
articulate as concepts. The paper proposes a sequence of these and exhibits
constructions of them unconnected to energy conservation or microscopic
reversibility.

Essential concepts are of three kinds: 1) the nature and origin of
macrostates; 2) the roles of entropy in relation to irreversibility and
fluctuation; and 3) the natural apportionment of irreversibility between a
system and its environment, which defines an intrinsic thermodynamics for the
system and an extrinsic thermodynamics that embeds the system in its context,
in analogy to the way differential geometry defines an intrinsic curvature for
a manifold distinct from extrinsic curvatures that may embed the manifold in
another manifold of higher dimension.444This analogy is not a reference to
natural information geometries Amari:inf_geom:01 ; Ay:info_geom:17  that can
also be constructed, though perhaps those geometries would provide an
additional layer of semantics to the constructions here.

All of these concepts originate in the large-deviation properties of
stochastic processes with multi-level timescale structure. They will be
demonstrated here using a simple class of stochastic population processes,
further specified as Chemical Reaction Network models as more complex
relations need to be presented.

Emphasis will be laid on the different roles of entropy as a functional on
general distributions versus entropy as a state function on macrostates, and
on the Lyapunov role of entropy in the 2nd law Boltzmann:second_law:86 ;
Fermi:TD:56  versus its large-deviation role in fluctuations Ellis:ELDSM:85 ,
through which the state-function entropy is most generally definable
Touchette:large_dev:09 . Both Shannon’s information entropy Shannon:MTC:49
and the older Hartley function Hartley:information:28 , will appear here as
they do in stochastic thermodynamics generally Seifert:stoch_thermo_rev:12 .
The different meanings these functions carry, which sound contradictory when
described and can be difficult to compare in derivations with different aims,
become clear as each appears in the course of a single calculation. Most
important, they have unambiguous roles in relation to either large deviations
or system decomposition without need of a reference to heat.

#### Main results and order of the derivation

Sec. II explains what is meant by multi-level systems with respect to a robust
separation of timescales, and introduces a family of models constructed
recursively from nested population processes. Macrostates are related to
microstates _within_ levels, and if timescale separations arise that create
new levels, they come from properties of a subset of long-lived, metastable
macrostates. Such states within a level, and transitions between them that are
much shorter than their characteristic lifetimes, map by coarse-graining to
the elementary states and events at the next level.

The environment of a system that arises in some level of a hierarchical model
is understood to include both the thermalized substrate at the level below,
and other subsystems with explicit dynamics within the same level. Sec. II.3
introduces the problem of system/environment partitioning of changes in
entropy, and states555The technical requirements are either evident or known
from fluctuation theorems Speck:HS_heat_FT:05 ; Harris:fluct_thms:07 ; variant
proofs are given in later sections. Their significance for system
decomposition is the result of interest here. the first main claim: the
natural partition employs a relative entropy Cover:EIT:91  within the system
and a housekeeping entropy rate Hatano:NESS_Langevin:01  to embed the system
in the environment. It results in two independently non-negative entropy
changes (plus a third entirely within the environment that is usually
ignored), which define the intensive thermodynamic of the focal system and the
extensive, or embedding thermodynamics of the system into the
$\mbox{system}\otimes\mbox{environment}$ whole.

Sec. III expresses the relation of macrostates to microstates in terms of the
large-deviations concept of separation of scale from structure
Touchette:large_dev:09 . Large-deviations scaling, defined as the convergence
of distributions for aggregate statistics toward exponential families, creates
a formal concept of a macroworld having definite structure, yet separated by
an indefinite or even infinite range of scale from the specifications of
micro-worlds.

Large-deviations for population processes are handled with the Hamilton-Jacobi
theory Bertini:macro_NEThermo:09  following from the time dependence of the
cumulant-generating function (CGF), and its Legendre duality Amari:inf_geom:01
to a large-deviation function called the effective action Smith:LDP_SEA:11 .
Macrostates are identified with the distributions that can be assigned large-
deviation probabilities from a system’s stationary distribution. Freedom to
study any CGF makes this definition, although concrete, flexible enough to
require choosing what Gell-Mann and Lloyd GellMann:EC:96 ;
GellMann:eff_complx:04  term “a judge”. The resulting definition, however,
does not depend on whether the system has any underlying mechanics or
conservation laws, explicit or implicit.

To make contact with multi-level dynamics and problems of interest in
stochastic thermodynamics Seifert:stoch_thermo_rev:12 ;
Polettini:open_CNs_I:14 ; Polettini:stoch_macro_thermo:16 , while retaining a
definite notation, Sec. IV assumes generators of the form used for Chemical
Reaction Networks (CRNs) Horn:mass_action:72 ; Feinberg:notes:79 ;
Krishnamurthy:CRN_moments:17 ; Smith:CRN_moments:17 . The finite-to-infinite
mapping that relates macrostates to microstates has counterparts for CRNs in
maps from the generator matrix to the transition matrix, and from mass-action
macroscopic currents to probability currents between microstates.

This section formally distinguishes the Lyapunov Fermi:TD:56  and large-
deviation Ellis:ELDSM:85  roles of entropy, and shows how the definition of
macrostates from Sec. III first introduces scale-dependence in the state
function entropy that was not present in the Lyapunov function for arbitrary
multi-scale models of Sec. II. In the Hamilton-Jacobi representation, Lyapunov
and large-deviation entropy changes occur within different manifolds and have
different interpretations. These differences are subordinate to the more
fundamental difference between the entropy state function and entropy as a
functional on arbitrary distributions over microstates. A properly formulated
2nd law, which is never violated, is computed from both deterministic and
fluctuation macrostate-entropies. The roles of Hartley informations
Hartley:information:28  for stationary states in stochastic thermodynamics
Speck:HS_heat_FT:05 ; Seifert:FDT:10 ; Seifert:stoch_thermo_rev:12  enter
naturally as macrostate relative entropies.

Sec. V derives the proofs of monotonicity of the intrinsic and extrinsic
entropy changes from Sec. II, using a cycle decomposition of currents in the
stationary distribution.666The decomposition is related to that used by
Schnakenberg Schnakenberg:ME_graphs:76  to compute dissipation in the
stationary state, but more cycles are required in order to compute dynamical
quantities. Whether only cycles or more complex hyperflows
Andersen:generic_strat:14  are required in a basis for macroscopic currents
distinguishes complexity classes for CRNs. Because cycles are _always_ a
sufficient basis in the microstate space, the breakdown of a structural
equivalence between micro- and macro-states for complex CRNs occurs in just
the terms responsible for the complex relations between state-function and
global entropies derived in Sec. IV.

Sec. VI illustrates two uses of the intensive/extensive decomposition of
entropy changes. It studies a simple model of polymerization and hydrolysis
under competing spontaneous and driven reactions in an environment that can be
in various states of disequilibrium. First a simple linear model of the kind
treated by Schnakenberg Schnakenberg:ME_graphs:76  is considered, in a limit
where one elementary reaction becomes strictly irreversible. Although the
housekeeping entropy rate becomes uninformative because it is referenced to a
diverging chemical potential, this is a harmless divergence analogous to a
scale divergence of an extensive potential. A Legendre dual to the entropy
rate remains regular, reflecting the existence of a thermodynamics-of-events
about which energy conservation, although present on the path to the limit,
asymptotically is not a source of any relevant constraints.

A second example introduces autocatalysis into the coupling to the
disequilibrium environment, so that the system can become bistable. The
components of entropy rate usually attributed to the “environment” omit
information about the interaction of reactions responsible for the
bistability, and attribute too much of the loss of large-deviation
accessibility to the environment. Housekeeping entropy rate gives the correct
accounting, recognizing that part of the system’s irreversibility depends on
the measure for _system_ relative entropy created by the interaction.

Sec.VII offers an alternative characterization of thermodynamic descriptions
when conservation and reversibility are not central. In 20th century physics
the problem of the nature and source of macroworlds has taken on a clear
formulation, and provides an alternative conceptual center for thermodynamics
to relations between work and heat. System decomposition and the conditional
independence structures within the loss of large-deviation accessibility
replace adiabatic transformation and heat flow as the central abstractions to
describe irreversibility, and define the entropy interpretation of the Hartley
information. Such a shift in view will be needed if path ensembles are to be
put on an equal footing with state ensembles to create a fully non-equilibrium
thermodynamics. The alternative formulation is meant to support the
development of a native thermodynamics of stochastic rule-based systems and a
richer phenomenology of living states.

##  II Multi-level systems

To provide a concrete class of examples for the large-deviation relations that
can arise in multi-scale systems, we consider systems with natural levels,
such that within each level a given system may be represented by a discrete
population process. The population process in turn admits leading-exponential
approximations of its large-deviation behavior in the form of a continuous
dynamical system. These dynamical systems are variously termed momentum-space
WKB approximations Assaf:mom_WKB:17 , Hamilton-Jacobi representations
Bertini:macro_NEThermo:09 , or eikonal expansions Freidlin:RPDS:98 . They
exist for many other classes of stochastic process besides the one assumed
here Martin:MSR:73 ; Kamenev:DP:02 , and may be obtained either directly from
the Gärtner-Ellis theorem Ellis:ELDSM:85  for cumulant-generating functions
(the approach taken here), by direct WKB asymptotics, or through saddle-point
methods in 2-field functional integrals.

A connection between levels is made by supposing that the large-deviation
behavior possesses one or more fixed points of the dynamical system. These are
coarse-grained to become the elementary states in a similar discrete
population process one level up in the scaling hierarchy. Nonlinearities in
the dynamical system that produce multiple isolated, metastable fixed points
are of particular interest, as transitions between these occur on timescales
that are exponentially stretched, in the large-deviation scale factor,
relative to the relaxation times. The resulting robust criterion for
separation of timescales will be the basis for the thermodynamic limits that
distinguish levels and justify the coarse-graining.

####  II.0.1 Micro to macro, within and between levels

Models of this kind may be embedded recursively in scale through any number of
levels. Here we will focus on three adjacent levels, diagrammed in Fig. 1, and
on the separation of timescales within a level, and the coarse-grainings that
define the maps between levels. The middle level, termed the mesoscale, will
be represented explicitly as a stochastic process, and all results that come
from large-deviations scaling will be derived within this level. The
microscale one level below, and the macroscale one level above, are described
only as needed to define the coarse-graining maps of variables between
adjacent levels. Important properties such as bidirectionality of escapes from
metastable fixed point in the stationary distribution, which the large-
deviation analysis in the mesoscale supplies as properties of elementary
transitions in the macroscale, will be assumed self-consistently as inputs
from the microscale to the mesoscale.

Figure 1:  A multiscale population process with a lattice of fixed points.
Middle layer is the grid of microstates and transitions in the mesoscale.
Surface $-\log\underline{\rho}$ for a stationary distribution is indicated by
the color gradient. Fixed points (asterisks) are at well bottoms. Modes of
$\underline{\rho}$ concentrated near fixed points are indicated in greyscale
beneath the surface. A classical deterministic trajectory (black, with arrow)
starts at the open circle and relaxes to the fixed point for that basin.
First-passage trajectories (colors) then mediate well switching by rare large
deviations. Fixed points of trajectory equations in the mesoscale are coarse-
grained to become elementary microstates (black dots) in the macroscale (top
layer), and first passages become elementary transition events (arrows).

Within the description at a single level, the discrete population states will
be termed microstates, as is standard in (both classical and stochastic)
thermodynamics. Microstates are different in kind from macrostates, which
correspond to a sub-class of distributions (defined below), and from fixed
points of the dynamical system. There is no unique recipe for coarse-graining
descriptions to reduce dimensionality in a multi-scale system, because the
diversity of stochastic processes is vast. Here, to obtain a manageable
terminology and class of models, we limit to cases in which the dynamical-
system fixed points at one level can be put in correspondence with the
microstates at the next higher level. Table 1 shows the terms that arise
within, and correspondences between, levels.

macroscale | mesoscale | microscale
---|---|---
|  microstate |  H-J fixed point
|  thermalization of the microscale |  H-J relaxation trajectories
|  elementary transitions between microstates |  H-J first-passages
|  arbitrary distribution on microstates |  coarse-grained distribution on fixed points
|  macrostate (distribution) $\leftrightarrow$ H-J variables |
microstate |  H-J fixed point |
… |  … |
Table 1:  Structures from micro to macro within and between levels. H-J refers
to the Hamilton-Jacobi dynamical system that arises in the large-deviation
approximation. Note that the distinction between microstate and macrostate is
a distinction of kind, which occurs within a level, whereas microscale,
mesoscale, and macroscale are designations of relative levels with respect to
large-deviation timescale separations. Terms in the same row are
representations of the same quantity in different levels, related through
coarse-graining.

###  II.1 Models based on population processes

The following elements furnish a description of a system:

##### The multiscale distribution:

The central object of study is any probability distribution $\rho$, defined
down to the smallest scale in the model, and the natural coarse-grainings of
$\rho$ produced by the dynamics. To simplify notation, we will write $\rho$
for the general distribution, across all levels, and let the indexing of
$\rho$ indicate which level is being used in a given computation.

##### Fast relaxation to fixed points at the microscale:

The counterpart, in our analysis of the mesoscale, to Prigogine’s
Glansdorff:structure:71 ; Prigogine:MT:98  assumption of local equilibrium in
a bath, is fast relaxation of the distribution on the microscale, to a
distribution with modes around the fixed points. In the mesoscale these fixed
points become the elementary population states indexed ${\rm n}$, and the
coarse-grained probability distribution is denoted ${\rho}_{{\rm n}}$. If a
population consists of individuals of types indexed
$p\in\left\\{1,\ldots,P\right\\}$, then ${\rm n}\equiv\left[{{\rm
n}}_{p}\right]$ is a vector in which the non-negative integer coefficient
${{\rm n}}_{p}$ counts the number of individuals of type $p$.

##### Elementary transitions in the mesoscale:

First-passages occurring in the (implicit) large-deviations theory, between
fixed points ${{\rm n}}^{\prime}$ and ${\rm n}$ at the microscale, appear in
the mesoscale as elementary transitions ${{\rm n}}^{\prime}\rightarrow{\rm
n}$. In Fig. 1, the elementary states are grid points and elementary
transitions occur along lines in the grid in the middle layer. We assume as
input to the mesoscale the usual condition of weak reversibility
Polettini:open_CNs_I:14 , meaning that if an elementary transition ${{\rm
n}}^{\prime}\rightarrow{\rm n}$ occurs with nonzero rate, then the transition
${\rm n}\rightarrow{{\rm n}}^{\prime}$ also occurs with nonzero rate. Weak
reversibility is a property we will derive for first passages within the
mesoscale, motivating its adoption at the lower level.

##### Thermalization in the microscale:

The elementary transitions ${{\rm n}}^{\prime}\leftrightarrow{\rm n}$ are
separated by typical intervals exponentially longer in some scale factor than
the typical time in which a single transition completes. (In the large-
deviation theory, they are instantons.) That timescale separation defines
thermalization at the microscale, and makes microscale fluctuations
conditionally independent of each other, given the index ${\rm n}$ of the
basin in which they occur. Thermalization also decouples components ${{\rm
n}}_{p}$ in the mesoscale at different $p$ except through the allowed
elementary state transitions.

##### A System/Environment partition within the mesoscale:

Generally, in addition to considering the thermalized microscale a part of the
“environment” in which mesoscale stochastic events take place, we will choose
some partition of the type-indices $p$ to distinguish one subset, called the
system ($s$), from one or more other subsets that also form part of the
environment ($e$). Unlike the thermalized microscale, the environment-part in
the mesoscale is slow and explicitly stochastic, like the system. The vector
${\rm n}$ indexes a tensor product space $s\otimes e$, so we write ${\rm
n}\equiv\left({{\rm n}}_{s},{{\rm n}}_{e}\right)$.

##### Marginal and conditional distributions in system and environment:

On $s\otimes e$, ${\rho}_{{\rm n}}$ is a joint distribution. A marginal
distribution ${\rho}^{s}$ for the system is defined by ${\rho}^{s}_{{{\rm
n}}_{s}}\equiv\sum_{{\rm n}\mid{{\rm n}}_{s}}{\rho}_{{\rm n}}$, where ${\rm
n}\mid{{\rm n}}_{s}$ fixes the $s$ component of ${\rm n}$ in the sum. From the
joint and the marginal a conditional distribution ${\rho}^{e\mid s}$ at each
${{\rm n}}_{s}$ is given by ${\rho}_{{\rm n}=\left({{\rm n}}_{s},{{\rm
n}}_{e}\right)}\equiv{\rho}^{s}_{{{\rm n}}_{s}}{\rho}^{e\mid s}_{{\rm n}}$.
The components of ${\rho}^{e\mid s}$ can fill the role often given to
chemostats in open-system models of CRNs. Here we keep them as explicit
distributions, potentially having dynamics that can respond to changes in $s$.

##### Notations involving pairs of indices:

Several different sums over the pairs of indices associated with state
transitions appear in the following derivations. To make equations easier to
read, the following notations are used throughout:

* $\left<{\rm n},{{\rm n}}^{\prime}\right>$ is an unordered pair of indices.

* $\sum_{{\rm n}}\sum_{{{\rm n}}^{\prime}\neq{\rm n}}$ counts every pair in both orders.
$\sum_{\left<{\rm n},{{\rm n}}^{\prime}\right>}$ counts every unordered pair
once.

* Therefore for any function $f\\!\left({\rm n},{{\rm n}}^{\prime}\right)$, $\sum_{{\rm n}}\sum_{{{\rm n}}^{\prime}\neq{\rm n}}f\\!\left({\rm n},{{\rm n}}^{\prime}\right)=\sum_{\left<{\rm n},{{\rm n}}^{\prime}\right>}\left[f\\!\left({\rm n},{{\rm n}}^{\prime}\right)+f\\!\left({{\rm n}}^{\prime},{\rm n}\right)\right]$.

* $\sum_{{\rm n}\mid{{\rm n}}_{s}}$ is a sum on the ${{\rm n}}_{e}$ component of ${\rm n}=\left({{\rm n}}_{s},{{\rm n}}_{e}\right)$.

* $\sum_{\left<{\rm n},{{\rm n}}^{\prime}\right>\mid{{\rm n}}_{s}}$ counts all unordered pairs $\left<{\rm n},{{\rm n}}^{\prime}\right>$ with common $s$-component ${{\rm n}}_{s}$.

###  II.2 Stochastic description within the mesoscale

The coarse-grained distribution ${\rho}_{{\rm n}}$ at the mesoscale evolves in
time under a master equation

$\displaystyle{\dot{\rho}}_{{\rm n}}=\sum_{{{\rm
n}}^{\prime}}{\mathbb{T}}_{{\rm n}{{\rm n}}^{\prime}}{\rho}_{{{\rm
n}}^{\prime}}$ $\displaystyle=\sum_{{{\rm n}}^{\prime}\neq{\rm
n}}\left(w_{{\rm n}{{\rm n}}^{\prime}}{\rho}_{{{\rm n}}^{\prime}}-w_{{{\rm
n}}^{\prime}{\rm n}}{\rho}_{{\rm n}}\right).$ (1)

Here and below, $(\dot{\mbox{ }})$ indicates the time derivative. The
generator $\mathbb{T}$ is a stochastic matrix on the left, which we write
$1^{T}\mathbb{T}=0$, where $1^{T}$ is the row-vector on index ${\rm n}$
corresponding to the uniform (unnormalized) measure. $w_{{\rm n}{{\rm
n}}^{\prime}}$ (following a standard notation Polettini:open_CNs_I:14 ) is the
component of $\mathbb{T}$ giving the transition rate from state ${{\rm
n}}^{\prime}$ to state ${\rm n}$.

For all of what follows, it will be necessary to restrict to systems that
possess a stationary, normalizable distribution denoted $\underline{\rho}$,
satisfying $\mathbb{T}\underline{\rho}=0$. The stationary distribution will
take the place of conservation laws as the basis for the definition of
macrostates. Moreover, if $\underline{\rho}$ is everywhere continuous, and the
number of distinct events generating transitions in the mesoscale (in a sense
made precise below) is finite, first passages between fixed points
corresponding to modes of $\underline{\rho}$ will occur at rates satisfying a
condition of detailed balance. The joint, marginal, and conditional stationary
distributions are denoted ${\underline{\rho}}_{{\rm n}=\left({{\rm
n}}_{s},{{\rm n}}_{e}\right)}\equiv{\underline{\rho}}^{s}_{{{\rm
n}}_{s}}{\underline{\rho}}^{e\mid s}_{{\rm n}}$.

##### The marginal stochastic process on $s$:

The system-marginal distribution ${\rho}^{s}$ evolves under a master equation
${\dot{\rho}}^{s}={\mathbb{T}}^{s}\\!\left({\rho}^{e\mid s}\right){\rho}^{s}$,
for which the transition matrix ${\mathbb{T}}^{s}$ has components that are
functions of the instantaneous environmental distribution, given by

$\displaystyle w^{s}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}\\!\left({\rho}^{e\mid s}\right)$
$\displaystyle\equiv\sum_{{\rm n}\mid{{\rm n}}_{s}}\sum_{{{\rm
n}}^{\prime}\mid{{\rm n}}^{\prime}_{s}}w_{{\rm n}{{\rm
n}}^{\prime}}{\rho}^{e\mid s}_{{{\rm n}}^{\prime}}.$ (2)

##### Time-independent overall stationary distribution as a reference:

Now we make an assumption that is crucial to being able to define a Lyapunov
function for the whole multi-level system distribution $\rho$: namely, that
the parameters in the mesoscale transition matrix $\mathbb{T}$ and hence the
stationary distribution $\underline{\rho}$ are time-independent. All dynamics
is made explicit as dynamics of distributions ${\rho}^{s}$ and ${\rho}^{e\mid
s}$, but the explicitly-written distributions are the _only_ source of
dynamics: once the microscale has thermalized, there are no further un-written
sources of time-dependence in the system.

##### Detailed balance propagating up from the microscale:

Finally we assume, propagating up from the microscale, a condition of detailed
balance that we will prove as a property of first-passages in the mesoscale,
and then apply recursively:

$\displaystyle w_{{\rm n}{{\rm n}}^{\prime}}{\underline{\rho}}_{{{\rm
n}}^{\prime}}$ $\displaystyle=w_{{{\rm n}}^{\prime}{\rm
n}}{\underline{\rho}}_{{\rm n}}.$ (3)

Note that condition (3) is _not_ an assumption of microscopic reversibility in
whatever faster stochastic process is operating below in the microscale. To
understand why, using constructions that will be carried out explicitly within
the mesoscale, see that even with rate constants satisfying Eq. (3), the
system-marginal transition rates (2) need not satisfy a condition of detailed
balance. Indeed we will want to be able to work in limits for the
environment’s conditional distributions ${\rho}^{e\mid s}$ in which some
transitions can be made completely irreversible: that is $w^{s}_{{{\rm
n}}_{s}{{\rm n}}^{\prime}_{s}}\neq 0$ but $w^{s}_{{{\rm n}}^{\prime}_{s}{{\rm
n}}_{s}}=0$. Even from such irreversible dynamics among microstates, first-
passage rates with detailed balance in the stationary distribution will
result, and it is that property that is assumed in Eq. (3).

From these assumptions on $\mathbb{T}$ and $\underline{\rho}$, it follows that
the relative entropy of any distribution $\rho$ from the stationary
distribution $\underline{\rho}$ is non-decreasing,

$\displaystyle-\dot{D}\\!\left(\rho\parallel\underline{\rho}\right)$
$\displaystyle=\sum_{{\rm n}}\sum_{{{\rm n}}^{\prime}\neq{\rm
n}}\log\left(\frac{{\left(\rho/\underline{\rho}\right)}_{{{\rm
n}}^{\prime}}}{{\left(\rho/\underline{\rho}\right)}_{{\rm n}}}\right)w_{{\rm
n}{{\rm n}}^{\prime}}{\rho}_{{{\rm n}}^{\prime}}$ $\displaystyle=\sum_{{\rm
n}}\sum_{{{\rm n}}^{\prime}\neq{\rm n}}\log\left(\frac{w_{{\rm n}{{\rm
n}}^{\prime}}{\rho}_{{{\rm n}}^{\prime}}}{w_{{{\rm n}}^{\prime}{\rm
n}}{\rho}_{{\rm n}}}\right)w_{{\rm n}{{\rm n}}^{\prime}}{\rho}_{{{\rm
n}}^{\prime}}$ $\displaystyle=\sum_{\left<{\rm n},{{\rm
n}}^{\prime}\right>}\log\left(\frac{w_{{\rm n}{{\rm n}}^{\prime}}{\rho}_{{{\rm
n}}^{\prime}}}{w_{{{\rm n}}^{\prime}{\rm n}}{\rho}_{{\rm
n}}}\right)\left(w_{{\rm n}{{\rm n}}^{\prime}}{\rho}_{{{\rm
n}}^{\prime}}-w_{{{\rm n}}^{\prime}{\rm n}}{\rho}_{{\rm n}}\right)\geq 0.$ (4)

The result is elementary for systems with detailed balance, because each term
in the third line of Eq. (4) is individually non-negative. We use relative
entropy to refer to minus the Kullback-Leibler divergence of $\rho$ from
$\underline{\rho}$ Cover:EIT:91 , to follow the usual sign convention for a
non-decreasing entropy.

###  II.3 System-environment decompositions of the entropy change

The exchange of heat for work is central to classical thermodynamics because
energy conservation is a constraint on joint configurations across sub-
systems, either multiple thermal systems in contact or a mechanical subsystem
having only deterministic variables, some of which set the boundary conditions
on thermal subsystems that also host fluctuations. Only the state-function
entropy, however, is a “function” of energy in any sense, so the only notion
of a limiting partition of irreversible effects between subsystems derivable
from energy conservation is the one defined by adiabatic transformations
passing through sequences of macrostates.

In more general cases, with or without conservation laws, the boundary
conditions on a system are imposed only through the elements of the marginal
transition matrix ${\mathbb{T}}^{s}$. The problem remains, of understanding
how one subsystem can limit entropy change in another through a boundary, but
it is no longer organized with reference to adiabatic transformations.

We wish to understand what constitutes a thermodynamically natural
decomposition of the mesoscale process into a system and an environment. A
widely-adopted decomposition Seifert:stoch_thermo_rev:12 ;
Polettini:open_CNs_I:14 ; Polettini:stoch_macro_thermo:16  for systems with
energy conservation777The decomposition is the same one used to define an
energy cost of computation Landauer:IHGCP:61 ; Bennett:TC:82  by constructing
logical states through the analogues to heat engines Szilard:MD:29 ;
Smith:NS_thermo_II:08 . separates a Shannon entropy of ${\rho}_{{\rm n}}^{s}$
from heat generation associated with terms $-\log w^{s}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}$ by the local equilibrium assumption for the bath. We begin
by writing down this information/heat decomposition, and arguing that it is
not the natural partition with respect to irreversibility.

##### Entropy relative to the stationary state rather than Shannon entropy:

The following will differ from the usual construction in replacing Shannon
entropy with a suitable relative entropy, without changing the essence of the
decomposition. Two arguments can be given for preferring the relative entropy:
It would be clear, for a system with a continuous state space in which
${\rho}_{{\rm n}}$ would become a density, that the logarithm of a dimensional
quantity is undefined. Hence some reference measure is always implicitly
assumed. A uniform measure is not a coordinate-invariant concept, and a
measure that is uniform in one coordinate system makes those coordinates part
of the system specification. Since discrete processes are often used as
approximations to continuum limits, the same concerns apply. The more general
lesson is that a logarithmic entropy unit is always given meaning with respect
to some measure. Only for systems such as symbol strings, for which a
combinatorial measure on integers is the natural measure, is Shannon entropy
the corresponding natural entropy. For other cases, such as CRNs, the natural
entropy is relative entropy referenced to the Gibbs equilibrium
Smith:NS_thermo_I:08 , and its change gives the dissipation of chemical work.
For the processes described here, the counterpart to Shannon entropy that
solves these consistency requirements, but does not yet address the question
of naturalness, is the relative entropy referenced to the steady state
marginal ${\underline{\rho}}^{s}$. Its time derivative is given by

$\displaystyle-\dot{D}\\!\left({\rho}^{s}\parallel{\underline{\rho}}^{s}\right)$
$\displaystyle=\sum_{\left<{{\rm n}}_{s},{{\rm
n}}^{\prime}_{s}\right>}\log\left(\frac{{\left(\rho/\underline{\rho}\right)}^{s}_{{{\rm
n}}^{\prime}_{s}}}{{\left(\rho/\underline{\rho}\right)}^{s}_{{{\rm
n}}_{s}}}\right)\left(w^{s}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}{\rho}^{s}_{{{\rm n}}^{\prime}_{s}}-w^{s}_{{{\rm
n}}^{\prime}_{s}{{\rm n}}_{s}}{\rho}^{s}_{{{\rm n}}_{s}}\right).$ (5)

The quantity (5) need not be either positive or negative in general.

A second term that separates out of the change in total relative entropy (4)
comes from changes in environmental states through events that do not result
in net change of the system state.888Note $w_{{\rm n}{{\rm n}}^{\prime}}$ may
depend on the system index-component ${{\rm n}}_{s}$ shared by both ${\rm n}$
and ${{\rm n}}^{\prime}$, so these rates can depend on system state. Catalysis
acts through such dependencies. The relative entropy of the conditional
distribution ${\rho}^{e\mid s}$ at a particular index ${{\rm n}}_{s}$ from its
stationary reference ${\underline{\rho}}^{e\mid s}$ has time derivative

$\displaystyle-{\dot{D}}^{e\mid s}_{{{\rm n}}_{s}}\\!\left({\rho}^{e\mid
s}\parallel{\underline{\rho}}^{e\mid s}\right)=$
$\displaystyle\sum_{\left<{{\rm n}}_{s},{{\rm n}}^{\prime}_{s}\right>\mid{{\rm
n}}_{s}}\log\left(\frac{{\left(\rho/\underline{\rho}\right)}^{e\mid s}_{{{\rm
n}}^{\prime}}}{{\left(\rho/\underline{\rho}\right)}^{e\mid s}_{{\rm
n}}}\right)\left(w_{{\rm n}{{\rm n}}^{\prime}}{\rho}^{e\mid s}_{{{\rm
n}}^{\prime}}-w_{{{\rm n}}^{\prime}{\rm n}}{\rho}^{e\mid s}_{{\rm n}}\right)$
$\displaystyle\sum_{\left<{{\rm n}}_{s},{{\rm n}}^{\prime}_{s}\right>\mid{{\rm
n}}_{s}}\log\left(\frac{w_{{\rm n}{{\rm n}}^{\prime}}{\rho}^{e\mid s}_{{{\rm
n}}^{\prime}}}{w_{{{\rm n}}^{\prime}{\rm n}}{\rho}^{e\mid s}_{{\rm
n}}}\right)\left(w_{{\rm n}{{\rm n}}^{\prime}}{\rho}^{e\mid s}_{{{\rm
n}}^{\prime}}-w_{{{\rm n}}^{\prime}{\rm n}}{\rho}^{e\mid s}_{{\rm n}}\right).$
(6)

Unlike the change of system relative entropy (5), Eq. (6) is non-negative
term-by-term, in the same way as Eq. (4).

The remaining terms to complete the entropy change (4) come from joint
transformations in system and environment indices ${{\rm n}}_{s}$ and ${{\rm
n}}_{e}$, and in usual treatments have the interpretation of dissipated
heats.999When more than one environment transition couples to the same system
transition, there can be reasons to further partition these terms; an example
is given in Sec. VI.1. They are functions of the pair of indices $\left({{\rm
n}}^{\prime}_{s},{{\rm n}}_{s}\right)$. An average change in relative entropy
of the environment, over all processes that couple to a given system state-
change, is

$\displaystyle{\sigma}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}\\!\left({\rho}^{e\mid s}\right)$
$\displaystyle\equiv\frac{1}{w^{s}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}}\sum_{{\rm n}\mid{{\rm n}}_{s}}\sum_{{{\rm
n}}^{\prime}\mid{{\rm
n}}^{\prime}_{s}}\log\left(\frac{{\left(\rho/\underline{\rho}\right)}^{e\mid
s}_{{{\rm n}}^{\prime}}}{{\left(\rho/\underline{\rho}\right)}^{e\mid s}_{{\rm
n}}}\right)w_{{\rm n}{{\rm n}}^{\prime}}{\rho}^{e\mid s}_{{{\rm
n}}^{\prime}}.$ (7)

(Note that if we had wished to use the un-referenced Shannon entropy
$-\sum_{{{\rm n}}_{s}}{\rho}^{s}_{{{\rm n}}_{s}}\log{\rho}^{s}_{{{\rm
n}}_{s}}$ in place of the relative entropy (5) – for instance, in an
application to digital computing – we could shift the measures
${\underline{\rho}}^{s}$ to the dissipation term to produce what is normally
considered the “environmental” heat dissipation, given by

$\displaystyle{\sigma}^{\rm env}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}w^{s}_{{{\rm n}}_{s}{{\rm n}}^{\prime}_{s}}$
$\displaystyle\equiv\left[{\sigma}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}-\log\left(\frac{{\underline{\rho}}^{s}_{{{\rm
n}}^{\prime}_{s}}}{{\underline{\rho}}^{s}_{{{\rm
n}}_{s}}}\right)\right]w^{s}_{{{\rm n}}_{s}{{\rm n}}^{\prime}_{s}}$
$\displaystyle=\sum_{{\rm n}\mid{{\rm n}}_{s}}\sum_{{{\rm
n}}^{\prime}\mid{{\rm n}}^{\prime}_{s}}\log\left(\frac{w_{{\rm n}{{\rm
n}}^{\prime}}{\rho}^{e\mid s}_{{{\rm n}}^{\prime}}}{w_{{{\rm n}}^{\prime}{\rm
n}}{\rho}^{e\mid s}_{{{\rm n}}}}\right)w_{{\rm n}{{\rm
n}}^{\prime}}{\rho}^{e\mid s}_{{{\rm n}}^{\prime}}.$ (8)

The quantity (8) is regarded as a property of the environment (both slow
variables and the thermal bath) because it is a function only of the
transition rates $w_{{{\rm n}}^{\prime}{\rm n}}$ and of the marginal
distributions ${\rho}^{e\mid s}$.)

####  II.3.1 The information/heat decomposition of total relative-entropy
change

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

$\displaystyle-\dot{D}\\!\left(\rho\parallel\underline{\rho}\right)$
$\displaystyle=-\dot{D}\\!\left({\rho}^{s}\parallel{\underline{\rho}}^{s}\right)$
$\displaystyle\mbox{}+\sum_{{{\rm n}}_{s}}\sum_{{{\rm
n}}^{\prime}_{s}\neq{{\rm n}}_{s}}{\sigma}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}\\!w^{s}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}{\rho}^{s}_{{{\rm n}}^{\prime}_{s}}$
$\displaystyle\mbox{}-\sum_{{{\rm n}}_{s}}{\rho}^{s}_{{{\rm
n}}_{s}}{\dot{D}}^{e\mid s}_{{{\rm n}}_{s}}\\!\left({\rho}^{e\mid
s}\parallel{\underline{\rho}}^{e\mid s}\right).$ (9)

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 .101010More 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.

####  II.3.2 Relative entropy referencing the system steady state at
instantaneous parameters

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.111111This 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 heat121212These 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”.

$\displaystyle{\sigma}^{\rm HK}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}w^{s}_{{{\rm n}}_{s}{{\rm n}}^{\prime}_{s}}$
$\displaystyle\equiv\left[{\sigma}^{\rm env}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}+\log\left(\frac{{\bar{\rho}}^{s}_{{{\rm
n}}^{\prime}_{s}}}{{\bar{\rho}}^{s}_{{{\rm n}}_{s}}}\right)\right]w^{s}_{{{\rm
n}}_{s}{{\rm n}}^{\prime}_{s}}$ $\displaystyle=\sum_{{\rm n}\mid{{\rm
n}}_{s}}\sum_{{{\rm n}}^{\prime}\mid{{\rm
n}}^{\prime}_{s}}\log\left(\frac{{\bar{\rho}}^{s}_{{{\rm
n}}^{\prime}_{s}}{\rho}^{e\mid s}_{{{\rm n}}^{\prime}}/\underline{\rho}_{{{\rm
n}}^{\prime}}}{{\bar{\rho}}^{s}_{{{\rm n}}_{s}}{\rho}^{e\mid s}_{{{\rm
n}}}/\underline{\rho}_{{\rm n}}}\right)w_{{\rm n}{{\rm
n}}^{\prime}}{\rho}^{e\mid s}_{{{\rm n}}^{\prime}}.$ (10)

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

$\displaystyle-\dot{D}\\!\left(\rho\parallel\underline{\rho}\right)$
$\displaystyle=-\dot{D}\\!\left({\rho}^{s}\parallel{\bar{\rho}}^{s}\right)$
$\displaystyle\mbox{}+\sum_{{{\rm n}}_{s}}\sum_{{{\rm
n}}^{\prime}_{s}\neq{{\rm n}}_{s}}{\sigma}^{\rm HK}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}\\!w^{s}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}{\rho}^{s}_{{{\rm n}}^{\prime}_{s}}$
$\displaystyle\mbox{}-\sum_{{{\rm n}}_{s}}{\rho}^{s}_{{{\rm
n}}_{s}}{\dot{D}}^{e\mid s}_{{{\rm n}}_{s}}\\!\left({\rho}^{e\mid
s}\parallel{\underline{\rho}}^{e\mid s}\right).$ (11)

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.

####  II.3.3 Intrinsic and extrinsic thermodynamics

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.

####  II.3.4 System Hartley information as a temporal connection

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.

##  III Hamilton-Jacobi theory for large deviations

The central concepts in thermodynamics are that of the macrostate, and of the
entropy as a state function from which properties of macrostates and
constraints on their transformations are derived. In classical thermodynamics
Fermi:TD:56 ; Kittel:TP:80 , macrostates are introduced in association with
average values of conserved quantities (e.g. energy, particle numbers),
because conservation laws naturally generate conditional independence between
subsystems in contact, given a simple boundary condition (the partition of the
conserved quantity between them).

Here we wish to separate the construction that _defines_ a macrostate from the
properties that make one or another class of macrostates dynamically _robust_
in a given system (conservation laws are important for the latter). The
defining construction can be quite general, but it must in all cases create a
dimensional reduction by an indefinite (or in the limit infinite) factor, from
the dimensionality of the microstate space that is definite but arbitrarily
large, to the dimensionality of a space of macrostate variables that is fixed
and independent of the dimensionality of microstates. Only dimensional
reductions of this kind are compatible with the large-deviation definition of
macroworlds as worlds in which structure can be characterized asymptotically
separate from scale Touchette:large_dev:09 . Robustness can then be
characterized separately within the large-deviation analysis in terms of
closure approximations or spectra of relaxation times for various classes of
macrostates.

Dimensional reduction by an indefinite degree is achieved by associating
macrostates with particular classes of distributions over microstates: namely,
those distributions produced in exponential families to define generating
functions. The coordinates in the tilting weights that define the family are
independent of the dimension of the microstate space for a given family and
become the intensive state variables (see Amari:inf_geom:01 ). They are
related by Legendre transform to deviations that are the dual extensive state
variables. Relative entropies such as
$-D\\!\left(\rho\parallel\underline{\rho}\right)$, defined as functionals on
arbitrary distributions, and dual under Legendre transform to suitable
cumulant-generating functions, become state functions when restricted to the
distributions for macrostates. The extensive state variables are their
arguments and the intensive state variables their gradients. Legendre duality
leads to a system of Hamiltonian equations Bertini:macro_NEThermo:09  for time
evolution of macrostate variables, and from these the large-deviation scaling
behavior, timescale structure, and moment closure or other properties of the
chosen system of macrostates are derived.

In the treatment of Sec. II, the relative entropy increased deterministically
without reference to any particular level of system timescales, or the size-
scale factors associated with various levels. As such it fulfilled the
Lyapunov role of (minus) the entropy, but not the large-deviation role that is
the other defining characteristic of entropy Ellis:ELDSM:85 ;
Touchette:large_dev:09 . The selection of a subclass of distributions as
macrostates introduces level-dependence, scale-dependence, and the large-
deviation role of entropy, and lets us construct the relation between the
Lyapunov and large-deviation roles of entropy for macroworlds, which are
generally distinct. As a quantity capable of fluctuations, the macrostate
entropy can decrease along subsets of classical trajectories; these
fluctuations are the objects of study in stochastic thermodynamics. The
Hamiltonian dynamical system is a particularly clarifying representation for
the way it separates relaxation and fluctuation trajectories for macrostates
into distinct sub-manifolds. In distinguishing the unconditional versus
conditional nature of the two kinds of histories, it shows how these macro-
fluctuations are not “violations” of the 2nd law, but rather a partitioning of
the elementary events through which the only properly formulated 2nd law is
realized.131313See a brief discussion making essentially this point in Sec.
1.2 of Seifert:stoch_thermo_rev:12 . Seifert refers to the 2nd law as
characterizing “mean entropy production”, in keeping with other
interpretations of entropies such as the Hartley information in terms of heat.
The characterization adopted here is more categorical: the Hartley function
and its mean, Shannon information, are not quantities with the same
interpretation; likewise, the entropy change (11) is the only entropy relative
to the boundary conditions in $\mathbb{T}$ that is the object of a well-
formulated 2nd law. The “entropy productions” resulting from the Prigogine
local-equilibrium assumption are conditional entropies for macrostates defined
through various large-deviation functions, shown explicitly below.

###  III.1 Generating functions, Liouville equation, and the Hamilton-Jacobi
construction for saddle points

####  III.1.1 Relation of the Liouville operator to the cumulant-generating
function

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}$:

$\displaystyle a^{\dagger}_{p}a_{p}\left|{\rm n}\right)$ $\displaystyle={{\rm
n}}_{p}\left|{\rm n}\right).$ (12)

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,

$\displaystyle\frac{d}{dt}\left|\rho\right)$
$\displaystyle=-\mathcal{L}\\!\left(a^{\dagger},a\right)\left|\rho\right),$
(13)

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)$,141414We 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,

$\displaystyle e^{-\Gamma\left(\log z\right)}$
$\displaystyle=\left(0\right|e^{za}\left|\rho\right).$ (14)

$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

$\displaystyle{n_{\rho}}_{p}\\!\left(z\right)$
$\displaystyle\equiv\frac{\left(0\right|e^{za}a^{\dagger}_{p}a_{p}\left|\rho\right)}{\left(0\right|e^{za}\left|\rho\right)}.$
(15)

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

$\displaystyle\frac{\partial}{\partial t}e^{-\Gamma\left(\log z\right)}$
$\displaystyle=-\left(0\right|e^{za}\mathcal{L}\\!\left(a^{\dagger},a\right)\left|\rho\right).$
(16)

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

$\displaystyle-\frac{\partial}{\partial t}\Gamma\\!\left(\log z\right)$
$\displaystyle=-\mathcal{L}\\!\left(z,n_{\rho}\\!\left(z\right)/z\right).$
(17)

##### From coherent-state to number-potential coordinates:

With a suitable ordering of operators in $\mathcal{L}$ (with all lowering
operators to the right of any raising operator), $z$ is the exact value
assigned to $a^{\dagger}$ in the expectation (16), and only the value
$n_{\rho}\\!\left(z\right)/z$ for $a$ depends on saddle-point approximations.
In special cases, where $\left|\rho\right)$ are eigenstates of $a$ known as
coherent states, the assignment to $a$ is also exact. Therefore the arguments
of $\mathcal{L}$ in Eq. (17) are called coherent-state coordinates.

However, $\log z$, which we henceforth denote by
$\theta\equiv\left[{\theta}_{p}\right]$, is the affine coordinate system in
which the CGF is locally convex,151515$\theta$ also provides the affine
coordinate system in the exponential family, which defines contravariant
coordinates in information geometry Amari:inf_geom:01 . and it will be
preferable to work in coordinates $\left(\theta,n_{\rho}\right)$, which we
term number-potential coordinates, because for applications in chemistry
$\theta$ has the dimensions of a chemical potential. We abbreviate
exponentials and other functions acting component-wise on vectors as $z\equiv
e^{\theta}$, and simply assign $n_{\rho}\\!\left(\theta\right)\equiv
n_{\rho}\\!\left(z\right)$. Likewise, $\mathcal{L}$ is defined in either
coordinate system by mapping its arguments:
${\mathcal{L}}^{\rm(N-P)}\\!\left(\theta,n\right)\equiv{\mathcal{L}}^{\rm(CS)}\\!\left(z,n\\!\left(z\right)/z\right)$.

####  III.1.2 Legendre transform of the CGF

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

$\displaystyle-\frac{\partial\Gamma}{\partial\theta}=n_{\rho}\\!\left(\theta\right)$
(18)

gives the mean (15) in the tilted distribution.

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

$\displaystyle S_{\rm eff}\\!\left(n\right)$
$\displaystyle=\max_{\theta}\left\\{\theta
n+\Gamma\\!\left(\theta\right)\right\\}.$ (19)

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$,

$\displaystyle\frac{\partial S_{\rm eff}}{\partial
n}={\theta}_{\rho}\\!\left(n\right),$ (20)

and the CGF is obtained by inverse Legendre transform.

$\displaystyle-\Gamma\\!\left(\theta\right)$
$\displaystyle=\max_{n}\left\\{\theta n-S_{\rm
eff}\\!\left(n\right)\right\\}.$ (21)

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

$\displaystyle{\left.\frac{\partial S_{\rm eff}}{\partial
t}\right|}_{n}=\mathcal{L}\\!\left({\theta}_{\rho}\\!\left(n\right),n\right).$
(22)

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.)

##### Multiple modes, Legendre-Fenchel transform, and locally-defined extrema

Systems with interesting multi-level structure do not have $S_{\rm eff}$ or
$\rho$ globally convex, but rather only locally convex. For these the
Legendre-Fenchel transform takes the place of the Legendre transform in Eq.
(19), and if constructed with a single coordinate $\theta$, $-\Gamma$ may have
discontinuous derivative.

For these one begins, rather than with the CGF, with $S_{\rm eff}$, evaluated
as a line integral of Eq. (20) in basins around the stationary points
$n_{\rho}\\!\left(0\right)$ at $\theta=0$. Each such basin defines invertible
pair of functions ${\theta}_{\rho}\\!\left(n\mid
n_{\rho}\\!\left(0\right)\right)$ and $n_{\rho}\\!\left(\theta\mid
n_{\rho}\\!\left(0\right)\right)$. We will not be concerned with the large-
deviation construction for general distributions in this paper, which is
better carried out using a path integral. We return in a later section to the
special case of multiple metastable fixed points in the stationary
distribution, and the modes associated with the stationary distribution
$\underline{\rho}$, and provide a more complete treatment.

####  III.1.3 Hamiltonian equations of motion and the action

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

$\displaystyle\frac{\partial{\theta}_{\rho}\\!\left(n\right)}{\partial t}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial n}.$ (23)

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

$\displaystyle\frac{\partial n_{\rho}\\!\left(\theta\right)}{\partial t}$
$\displaystyle=-\frac{\partial\mathcal{L}}{\partial\theta}.$ (24)

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, as161616The same action
functionals are arrived at somewhat more indirectly via 2-field functional
integral constructions such as the Doi-Peliti method.

$\displaystyle S$ $\displaystyle=\int
dt\left\\{-\left(d_{t}z\right)\left(n/z\right)+{\mathcal{L}}^{\rm(CS)}\\!\left(z,n/z\right)\right\\}$
$\displaystyle=\int
dt\left\\{-\left(d_{t}\theta\right)n+{\mathcal{L}}^{\rm(N-P)}\\!\left(\theta,n\right)\right\\}.$
(25)

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 .

##### Circulation-free vector field of $\theta$ for the stationary
distribution

In order for $S_{\rm eff}$ to be a continuum approximation to $-\log\rho$, if
$\rho$ exists and is smooth everywhere, the vector field $\theta$ obtained
from stationary trajectories of the action (25) must have zero circulation in
order to be a well-defined gradient through Eq. (20). To check that this is
the case, consider the increment of $\theta$ under a small interval $dt$ under
Eq. (23):

$\displaystyle d{\theta}_{p}$ $\displaystyle\equiv
dt\frac{\partial\mathcal{L}}{\partial n_{p}}.$ (26)

The gradient in $n$ of $\theta$ therefore increments in time as

$\displaystyle\frac{\partial\left({\theta}_{p}+d{\theta}_{p}\right)}{\partial
n_{q}}$ $\displaystyle=\frac{\partial{\theta}_{p}}{\partial
n_{q}}+dt\frac{\partial{\mathcal{L}}^{2}}{\partial n_{p}\partial n_{q}}.$ (27)

Contraction of Eq. (27) with the antisymmetric symbol in $p$ and $q$ vanishes

$\displaystyle\frac{d}{dt}{\epsilon}_{pq}\frac{\partial{\theta}_{p}}{\partial
n_{q}}$ $\displaystyle=0,$ (28)

so the circulation of $\theta$ is the same everywhere as at the fixed points.

From Eq. (20), it is required to be the case that
$\partial{\theta}_{p}/\partial n_{q}={\partial}^{2}S_{\rm eff}/\partial
n_{p}\partial n_{q}\equiv g^{-1}_{pq}$, the inverse of the Fisher metric
Amari:inf_geom:01 , symmetric by construction and thus giving
${\epsilon}_{pq}g^{-1}_{pq}=0$. The only difference between the fixed point
and any other point is that for distant points we are relying on Hamiltonian
trajectories to evaluate $\theta$, whereas the the fixed point the Fisher
metric may be calculated by means not relying on the large-deviation saddle-
point approximation. Therefore Eq. (28) may be read as a check that symmetry
of the Fisher metric is preserved by Hamiltonian trajectories directly from
the symmetric partial derivative of $\mathcal{L}$ in Eq. (27).

###  III.2 The stationary distribution and macrostates

Up to this point only the Lyapunov role (4) of the relative entropy has been
developed. While the increase of
$-D\\!\left(\rho\parallel\underline{\rho}\right)$ has the appearance of the
classical 2nd law, we can understand from three observations that this
relative entropy is not the desired generalization of the entropy state
function of classical thermodynamics to express the phenomenology of multi-
level systems:

  1. 1.

The relative entropy $-D\\!\left(\rho\parallel\underline{\rho}\right)$ is a
functional on arbitrary distributions, like the Shannon entropy that is a
special case. It identifies no concept of macrostate, and has no dependence on
state variables.

  2. 2.

In a multilevel system that may have arbitrarily fine-grained descriptions,
there is no upper limit to $D\\!\left(\rho\parallel\underline{\rho}\right)$,
and no appearance of the system scale at any particular level, which
characterizes state-function entropies.

  3. 3.

Eq. (4) describes a deterministic increase of relative entropy; the large-
deviation role of entropy as a log-probability for macrostate fluctuations
Ellis:ELDSM:85  does not appear.

The step that has not yet been taken in our construction is, of course, the
identification of a macrostate concept. Here we depart from the usual
development based on conservation laws, and follow Gell-Mann and Lloyd
GellMann:EC:96 ; GellMann:eff_complx:04  in claiming that the concept of
macrostate is not inherent in features of a system’s dynamics, but requires
one to explicitly choose a procedure for aggregation or coarse-graining – what
they call a “judge” – as part of the commitment to which phenomenology is
being described.

We will put forth the definition of macrostates as _the tilted distributions
arising in generating functions for the stationary distribution
$\underline{\rho}$_. In the case of generating functions for number, these are
the distributions appearing in Eq. (14), which we will denote by
${\rho}^{\left(\bar{n}\right)}$. They are the least-improbable distributions
with a given non-stationary mean to arise through aggregate microscopic
fluctuations, and therefore dominate the construction of the large-deviation
probability.

The extensive state variable associated with this definition of macrostate is
the tilted mean from Eq. (18), which we will denote
$\bar{n}=n_{\underline{\rho}}\\!\left(\theta\right)$. If $\underline{\rho}$
and ${\rho}^{\left(\bar{n}\right)}$ are sharply peaked – the limit in which
the large-deviation approximation is informative – the relative entropy of the
macrostate is dominated at the saddle point of
${\rho}^{\left(\bar{n}\right)}$, where $\log{\rho}^{\left(\bar{n}\right)}\sim
0$, and thus

$\displaystyle
D\\!\left({\rho}^{\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
$\displaystyle\sim{\left.-\log{\underline{\rho}}_{{\rm n}}\right|}_{{\rm
n}\sim\bar{n}}={\underline{S}}_{\rm eff}\\!\left(\bar{n}\right).$ (29)

The general relative entropy functional, applied to the macrostate, becomes
the entropy state function ${\underline{S}}_{\rm eff}$, which takes as its
argument the extensive state variable $\bar{n}$. Moreover, because probability
under $\underline{\rho}$ is concentrated on configurations with the scale that
characterizes the system, tilted means $\bar{n}$ that are not suppressed by
very large exponential probabilities will have comparable scale. If $\bar{n}$
is also the scale factor in the large-deviations function (a property that may
or may not hold, depending on the system studied), then ${\underline{S}}_{\rm
eff}\sim\bar{n}$ in scale, and the entropy state function now has the
characteristic scale of the mesoscale level of the process description.

The three classes of distributions that enter a thermodynamic description are
summarized in Table 2.

notation | definition | comment
---|---|---
$\underline{\rho}$ , ${\underline{\rho}}^{s}$ , ${\underline{\rho}}^{e\mid s}$ |  whole-system, $s$-marginal, $\left(e\mid s\right)$-conditional distributions in the global stationary state |
${\bar{\rho}}^{s}$ |  marginal system steady-state distribution |  function of instantaneous environment conditional distribution ${\rho}^{e\mid s}$
${\rho}^{\left(\bar{n}\right)}$ |  macrostate tilted to saddle-point value $\bar{n}$ |  defined relative to global stationary distributions $\underline{\rho}$; may be defined for whole-system, $s$, or $e\mid s$
Table 2:  Notation conventions adopted for three classes of distributions
arising in large-deviation systems with multi-scale relaxation structure.

####  III.2.1 Coherent states, dimensional reduction, and the $f$-divergence

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

$\displaystyle-\log{\rho}^{\left(\bar{n}\right)}_{{\rm n}}$
$\displaystyle\approx{\rm n}\cdot\log\left(\frac{{\rm
n}}{\bar{n}}\right)-\left({\rm n}-\bar{n}\right)\cdot 1\equiv
D_{f}\\!\left({\rm n}\parallel\bar{n}\right).$ (30)

$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

$\displaystyle
D\\!\left({\rho}^{\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
$\displaystyle\approx D_{f}\\!\left(\bar{n}\parallel\underline{n}\right),$
(31)

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.

####  III.2.2 Multiple fixed points and instantons

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

$\displaystyle{\underline{S}}_{\rm eff}\\!\left(\bar{n}\right)$
$\displaystyle=\max_{\theta}\left\\{\theta\,\bar{n}+{\Gamma}_{\underline{n}}\\!\left(\theta\right)\right\\}.$
(32)

$\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

$\displaystyle{\left.-\frac{\partial{\Gamma}_{\underline{n}}}{\partial\theta}\right|}_{\theta\left(\bar{n}\right)}=\bar{n}.$
(33)

Next, let the function $\underline{n}\\!\left(\bar{n}\right)$171717Note 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

$\displaystyle
D\\!\left({\rho}^{\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
$\displaystyle\sim{\underline{S}}_{\rm
eff}\\!\left(\bar{n}\right)-{\underline{S}}_{\rm
eff}\\!\left(\underline{n}\\!\left(\bar{n}\right)\right)+D\\!\left(1_{\underline{n}\left(\bar{n}\right)}\parallel\underline{\rho}\right).$
(34)

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.

##  IV Population processes with CRN-form generators

The results up to this point apply to general processes with discrete state
spaces, normalizable stationary distributions, and some kind of
$\mbox{system}\otimes\mbox{environment}$ tensor-product structure on states.
For the next steps we restrict to stochastic population processes that can be
written in a form equivalent to Chemical Reaction Networks Horn:mass_action:72
; Feinberg:notes:79 ; Gunawardena:CRN_for_bio:03 . CRNs are an expressive
enough class to include many non-mechanical systems such as evolving Darwinian
populations Smith:evo_games:15 , and to implement algorithmically complex
processes Andersen:NP_autocat:12 . Yet they possess generators with a compact
and simple structure Krishnamurthy:CRN_moments:17 ; Smith:CRN_moments:17 , in
which similarities of microstate and macrostate phenomena are simply reflected
in the formalism.

###  IV.1 Hypergraph generator of state-space transitions

Population states ${\rm n}$ give counts of individuals, grouped according to
their types $p$ which are termed species. Stochastic Chemical Reaction
Networks assume independent elementary events grouped into types termed
reactions. Each reaction event removes a multiset181818A multiset is a
collection of distinct individuals, which may contain more than one individual
from the same species. of individuals termed a complex from the population,
and places some (generally different) multiset of individuals back into it.
The map from species to complexes is called the stoichiometry of the CRN.
Reactions occur with probabilities proportional per unit time to rates. The
simplest rate model, used here, multiplies a half-reaction rate constant by a
combinatorial factor for proportional sampling without replacement from the
population to form the complex, which in the mean-field limit leads to mass-
action kinetics.

Complexes will be indexed with subscripts $i,j,\ldots$, and an ordered pair
such as $ji$ labels the reaction removing $i$ and creating $j$. With these
conventions the Liouville-operator representation of the generator from Eq.
(13) takes the form Krishnamurthy:CRN_moments:17 ; Smith:CRN_moments:17

$\displaystyle-\mathcal{L}$
$\displaystyle=\sum_{ji}k_{ji}\left[{\psi}_{j}\\!\left(a^{\dagger}\right)-{\psi}_{i}\\!\left(a^{\dagger}\right)\right]{\psi}_{i}\\!\left(a\right)$
$\displaystyle\equiv{\psi}^{T}\\!\left(a^{\dagger}\right)\mathbb{A}\,\psi\\!\left(a\right).$
(35)

$\left\\{k_{ji}\right\\}$ are half-reaction rate constants, organized in the
second line into an adjacency matrix $\mathbb{A}$ on complexes.
$\psi\equiv\left[{\psi}_{i}\right]$ is a column vector with components
${\psi}_{i}\\!\left(a\right)\equiv\prod_{p}a_{p}^{y_{ip}}\equiv a^{Y_{i}}$
that, in the Doi algebra, produce the combinatorial factors and state shifts
reflecting proportional sampling without replacement. Here
$Y\equiv\left[y_{ip}\right]$ is the matrix of stoichiometric coefficients,
with entry $y_{ip}$ giving the number of individuals of species $p$ that make
up complex $i$. Further condensing notation, $Y_{i}\equiv\left[y_{ip}\right]$
are column vectors in index $p$, and $Y_{p}\equiv\left[y_{ip}\right]$ are row
vectors on index $i$. $a^{Y_{i}}$ is understood as the component-wise product
of powers $a_{p}^{y_{ip}}$ over the species index $p$. ${\psi}^{T}$ is a
transpose (row) vector, and vector and matrix products are used in the second
line of Eq. (35).

The generator (35) defines a hypergraph Berge:hypergraphs:73  in which
$\mathbb{A}$ is the adjacency matrix on an ordinary graph over complexes known
as the complex graph Gunawardena:CRN_for_bio:03 ; Baez:QTRN_eq:14 . The
stoichiometric vectors $\left\\{Y_{i}\right\\}$, defining complexes as
multisets, make the edges in $\mathbb{A}$ directed hyper-edges relative to the
population states that are the Markov states for the process.191919Properly,
the generator should be called a “directed multi-hypergraph” because the
complexes are multisets rather than sets. The concurrent removal or addition
of complexes is the source of both expressive power and analytic difficulty
provided by hypergraph-generators.

A master equation (1) acting in the state space rather than on the generating
function may be written in terms of the same operators, as

$\displaystyle{\dot{\rho}}_{{\rm n}}$
$\displaystyle={\psi}^{T}\\!\left(e^{-\partial/\partial{\rm
n}}\right)\mathbb{A}\operatorname{diag}\left[\psi\\!\left(e^{\partial/\partial{\rm
n}}\right)\right]\Psi\\!\left({\rm n}\right){\rho}_{{\rm n}}.$ (36)

Here a formal shift operator $e^{\partial/\partial{\rm
n}}\equiv\left[e^{\partial/\partial{{\rm n}}_{p}}\right]$ is used in place of
an explicit sum over shifted indices, and
${\psi}_{i}\\!\left(e^{\pm\partial/\partial{\rm n}}\right)=e^{\pm
Y_{i}^{T}\partial/\partial{\rm n}}$ creates shifts by the stoichiometric
vector $Y_{i}$. $\operatorname{diag}\left[\psi\right]$ refers to the matrix
with diagonal entries given by the components ${\psi}_{i}$.

In the master equation the combinatorial factors must be given explicitly.
These are written as a vector $\Psi\equiv\left[{\Psi}_{i}\right]$ with
components ${\Psi}_{i}\\!\left({\rm n}\right)\equiv\prod_{p}\left[{{\rm
n}}_{p}!/\left({{\rm n}}_{p}-y_{ip}\right)!\right]\equiv{{\rm
n}}^{\underline{Y_{i}}}$ that are falling factorials from ${\rm n}$, denoted
with the underscore as ${{\rm n}}^{\underline{Y_{i}}}$.

The matrix elements of Eq. (1) may be read off from Eq. (36) in terms of
elements in the hypergraph, as

$\displaystyle w^{\left(ji\right)}_{{{\rm n}}^{\prime}{\rm n}}$
$\displaystyle\equiv k_{ji}{\Psi}_{i}\\!\left({\rm n}\right)$
$\displaystyle\mbox{where}\qquad{{\rm n}}^{\prime}$
$\displaystyle=Y_{j}-Y_{i}+{\rm n},$ (37)

between all pairs $\left({{\rm n}}^{\prime},{\rm n}\right)$ separated by the
stoichiometric difference vector $Y_{j}-Y_{i}$. If multiple reactions produce
transitions between the same pairs of states, the aggregate rates become

$\displaystyle w_{{{\rm n}}^{\prime}{\rm n}}$
$\displaystyle\equiv\sum_{ji\;\mid\;Y_{j}-Y_{i}={{\rm n}}^{\prime}-{\rm
n}}w^{\left(ji\right)}_{{{\rm n}}^{\prime}{\rm n}}.$ (38)

In this way a finite adjacency matrix $\mathbb{A}$ on complexes may generate
an infinite-rank transition matrix $\mathbb{T}$, which is the adjacency matrix
for an ordinary graph over states. We will see that for CRNs, the hypergraph
furnishes a representation for macrostates similar to the representation given
by the simple graph for microstates.

From Eq. (38), marginal transition rates $w^{s}_{{{\rm n}}_{s}{{\rm
n}}^{\prime}_{s}}$ for the system may be defined using the average (2) over
${\rho}^{e\mid s}$. Note that the dependence of the activity products
${\Psi}_{i}$ on species $p$ within the system remains that of a falling
factorial, even if the average over activities of species in the environment
is complicated. Denote by ${\Psi}^{s}\\!\left({{\rm n}}_{s}\right)$ and
${\psi}^{s}\\!\left(e^{\pm\partial/\partial{{\rm n}}_{s}}\right)$ the
restrictions of the activity and shift functions to the species in the system
$s$, and by $k^{s}_{ji}\\!\left({\rho}^{e\mid s}\right)$ the rate constants
after computing the sum in Eq. (2) over the index for species in the
environment.

####  IV.1.1 Descaling of transition matrices for microstates

Proofs of monotonic change, whether of total relative entropy (1) or for the
components of entropy change partitioned out in Sec. II.3, take a particularly
simple form for CRNs generated by finitely many reactions. They make use of
finite cycle decompositions of the current through any microstate or complex,
which are derived from descaled adjacency matrices.

The descaling that produces a transition matrix that annihilates the uniform
measure on both the left and the right is

$\displaystyle\hat{\mathbb{T}}$
$\displaystyle\equiv{\psi}^{T}\\!\left(e^{-\partial/\partial{\rm
n}}\right)\mathbb{A}\operatorname{diag}\left[\psi\\!\left(e^{\partial/\partial{\rm
n}}\right)\right]\Psi\\!\left({\rm
n}\right)\operatorname{diag}\left[{\underline{\rho}}_{{\rm n}}\right]$ (39)

for the whole mesoscale, or

$\displaystyle{\hat{\mathbb{T}}}^{s}$
$\displaystyle\equiv{\psi}^{sT}\\!\left(e^{-\partial/\partial{{\rm
n}}_{s}}\right){\mathbb{A}}^{s}\operatorname{diag}\left[{\psi}^{s}\\!\left(e^{\partial/\partial{{\rm
n}}_{s}}\right)\right]{\Psi}^{s}\\!\left({{\rm
n}}_{s}\right)\operatorname{diag}\left[{\bar{\rho}}^{s}_{{{\rm
n}}_{s}}\right]$ (40)

for the subsystem $s$, by definition of the stationary states.
${\mathbb{A}}^{s}$ in Eq. (40) is the adjacency matrix on complexes that gives
rate constants $w^{s}_{{{\rm n}}_{s}{{\rm n}}^{\prime}_{s}}$ from Eq. (2).
These descalings are familiar as the ones leading to the dual time-reversed
generators in the fluctuation theorem for housekeeping heat
Speck:HS_heat_FT:05 .

We return to use the descaled microstate transition matrices (39,40) in
monotonicity proofs for general CRNs in Sec. V, but before doing that, we use
the descaling in the state space to motivate an analogous and simpler
descaling for macrostates at the level of the hypergraph. That descaling
illustrates the cycle decomposition on finitely many states, though it only
yields the $f$-Divergence of Eq. (31) as a Lyapunov function for the complex-
balanced CRNs.

####  IV.1.2 Descaling of transition matrices for macrostates

The coherent states, which are the moment-generating functions of Poisson
distributions and their duals in the Doi Hilbert space, are defined as

$\displaystyle\left({\phi}^{\dagger}\right|$ $\displaystyle\equiv
e^{-{\phi}^{\dagger}\phi}\left(0\right|e^{{\phi}^{\dagger}a}$
$\displaystyle\left|\phi\right)$ $\displaystyle\equiv
e^{a^{\dagger}\phi}\left|0\right),$ (41)

where $\phi\equiv\left[{\phi}_{p}\right]$ is a vector of (generally complex)
numbers and ${\phi}^{\dagger}$ is its Hermitian conjugate. They are
eigenstates respectively of the raising and lowering operators with
eigenvalues
$\left({\phi}^{\dagger}\right|a^{\dagger}_{p}=\left({\phi}^{\dagger}\right|{\phi}^{\ast}_{p}$
and $a_{p}\left|\phi\right)={\phi}_{p}\left|\phi\right)$.

On the constant trajectories corresponding to a fixed point of the Hamiltonian
equations (23, 24), ${\phi}^{\dagger}\equiv 1$ and
$\phi\equiv\underline{\phi}=\underline{n}$, the fixed-point number. We may
descale the coherent-state parameters $\phi$ and ${\phi}^{\dagger}$, which
correspond to classical state variables, by defining

$\displaystyle{\phi}_{p}$
$\displaystyle\equiv{\underline{\phi}}_{p}{\varphi}_{p}$
$\displaystyle{\phi}^{\ast}_{p}{\underline{\phi}}_{p}$
$\displaystyle\equiv{\varphi}^{\ast}_{p}.$ (42)

In vector notation,
$\phi\equiv\operatorname{diag}\left[\underline{\phi}\right]\varphi$,
${\phi}^{\dagger}\operatorname{diag}\left[\underline{\phi}\right]\equiv{\varphi}^{\dagger}$.
The scaling (42) was introduced by Baish Baish:DP_duality:15  to study
dualities of correlation functions in Doi-Peliti functional integrals.

The compensating descaling of the adjacency matrix

$\displaystyle\hat{\mathbb{A}}$
$\displaystyle\equiv\mathbb{A}\operatorname{diag}\left[\psi\\!\left(\underline{\phi}\right)\right]$
(43)

may be compared with Eq. (39) for $\mathbb{T}$. A similar descaling may be
done for ${\mathbb{A}}^{s}$ using ${\bar{\rho}}^{s}$. Henceforth we omit the
duplicate notation, and carry out proofs with respect to whatever is the
stationary distribution for a given adjacency matrix and hypergraph.

The coherent-state action (25) with Liouville operator (35) has a symmetric
form in descaled variables that is particularly useful for understanding the
relaxation and escape trajectories in the Hamilton-Jacobi system:

$\displaystyle S$ $\displaystyle=\int
dt\left\\{{\phi}^{\dagger}\operatorname{diag}\left[\underline{\phi}\right]d_{t}\varphi-{\psi}^{T}\\!\left({\phi}^{\ast}\right)\hat{\mathbb{A}}\,\psi\\!\left(\varphi\right)\right\\}.$
(44)

####  IV.1.3 Equations of motion and the $\mathcal{L}=0$ manifold

The Hamiltonian equations derived by variation of the action (44) can be
written

$\displaystyle\operatorname{diag}\left[\underline{\phi}\right]\dot{\varphi}$
$\displaystyle=\frac{\partial{\psi}^{T}}{\partial{\phi}^{\dagger}}\hat{\mathbb{A}}\psi\\!\left(\varphi\right)$
$\displaystyle\underset{{\phi}^{\dagger}=1}{\rightarrow}Y\hat{\mathbb{A}}\psi\\!\left(\varphi\right)$
$\displaystyle{\dot{\phi}}^{\dagger}\operatorname{diag}\left[\underline{\phi}\right]$
$\displaystyle=-{\psi}^{T}\\!\left({\phi}^{\ast}\right)\hat{\mathbb{A}}\frac{\partial{\psi}}{\partial\varphi}$
$\displaystyle\underset{\varphi=1}{\rightarrow}-{\psi}^{T}\\!\left({\phi}^{\ast}\right)\hat{\mathbb{A}}Y^{T}.$
(45)

The general solution is shown first in each line and then particular limiting
forms are shown.

The sub-manifold ${\phi}^{\dagger}\equiv 1$ contains all relaxation
trajectories for any CRN. It is dynamically stable because $\mathbb{A}$ is a
stochastic matrix on the left, and thus
${\psi}^{T}\\!\left(1\right)\hat{\mathbb{A}}=0$ in the second line of Eq.
(45). The image of $Y\mathbb{A}$ is called the stoichiometric subspace, and
its dimension is denoted
$s\equiv\dim\left[\operatorname{im}\left(Y\mathbb{A}\right)\right]$.

An important simplifying property of some CRNs is known as complex balance of
the stationary distributions.202020Complex balance can be ensured at all
parameters if a topological character of the CRN known as deficiency equals
zero, but may also be true for nonzero-deficiency networks for suitably tuned
parameters. In this work nothing requires us to distinguish these reasons for
complex balance. It is the condition that the fixed point
$\psi\\!\left(\underline{n}\right)\in\ker\mathbb{A}$ and not only $\in\ker
Y\mathbb{A}$. Since $\underline{n}=\underline{\phi}$ corresponds to
$\varphi=1$ and thus $\psi\\!\left(\varphi\right)=1$,
$1\in\ker\hat{\mathbb{A}}$ and $\dot{\varphi}=0$ in the first line of Eq. (45)
at any ${\phi}^{\dagger}$. Complex-balanced CRNs (with suitable conditions on
$\mathbb{A}$ to ensure ergodicity on the complex network
Gunawardena:CRN_for_bio:03 ) always possess unique interior fixed points
Feinberg:notes:79 ; Feinberg:def_01:87  and simple product-form stationary
distributions at these fixed points Anderson:product_dist:10 .

For non-complex-balanced stationary solutions, although escapes may have
$\varphi=1$ as initial conditions, that value is not dynamically maintained.
Recalling the definition (12) of the number operator, the field
$n_{p}={\phi}^{\dagger}_{p}{\phi}_{p}$, so non-constant $\phi$ is required for
instantons to escape from stable fixed points and terminate in saddle fixed
points, in both of which limits ${\phi}^{\dagger}\rightarrow 1$.

All relaxations and also the escape trajectories from fixed points (the
instantons) share the property that $\mathcal{L}\equiv 0$ for a CRN with time-
independent parameters.212121We can see that this must be true because
$-\mathcal{L}$ is a Hamiltonian conserved under the equations of motion, and
because instantons trace the values of $n$ and $\theta$ in the stationary
distribution $\underline{\rho}$, which must then have $\mathcal{L}=0$ to
satisfy the Hamilton-Jacobi equation (22). This submanifold separates into two
branches, with ${\phi}^{\dagger}\equiv 1$ for relaxations and
${\phi}^{\dagger}\neq 1$ for escapes.

####  IV.1.4 The Schlögl cubic model to illustrate

The features of CRNs with multiple metastable fixed points are exhibited in a
cubic 1-species model introduced by Schlögl Schlogl:near_SS_thermo:71 , which
has been extensively studied Dykman:chem_paths:94 ;
Krishnamurthy:CRN_moments:17 ; Smith:CRN_moments:17  as an example of
dynamical bistability.

The reaction schema in simplified form is

$\displaystyle\varnothing$
$\displaystyle\overset{k_{d}}{\underset{{\bar{k}}_{d}}{\rightleftharpoons}}A$
$\displaystyle 2A$
$\displaystyle\overset{k_{c}}{\underset{{\bar{k}}_{c}}{\rightleftharpoons}}3A.$
(46)

We choose rate constants so that the mass-action equation for the number of
$A$ particles, given by $n$, is

$\dot{n}=\left(n_{3}-n\right)\left(n_{2}-n\right)\left(n_{1}-n\right).$ (47)

The three fixed points are
$\underline{n}\in\left\\{n_{1},n_{2},n_{3}\right\\}$, of which $n_{1}$ and
$n_{3}$ are stable, and $n_{2}$ is a saddle.

The relaxation and escape branches of the $\mathcal{L}=0$ manifold are shown
in Fig. 2. Because the Schlögl model is 1-dimensional, the condition
$\mathcal{L}=0$ fixes $\theta\\!\left(n\right)$ along escape trajectories.
Because the stochastic model is a birth-death process, it is also exactly
solvable vanKampen:Stoch_Proc:07 . We will return to the example in Sec. VI.2
to study properties of the intensive and extensive thermodynamic potentials
for it.

Figure 2:  The two branches of the $\mathcal{L}=0$ manifold, corresponding to
the rate equation (47). ${\phi}^{\dagger}\equiv 1$ (blue);
${\phi}^{\dagger}\neq 1$ (red). Trajectories at common $n$ have the same value
on the vertical axis. Time is shown moving toward the front.

###  IV.2 Large-deviation and Lyapunov roles of the effective action

The role of entropy in the understanding of Boltzmann and Gibbs was that of a
Lyapunov function Fermi:TD:56 ; Kittel:TP:80 , accounting for
unidirectionality from microscopically reversible mechanics. The much later
understanding of entropy in relation to fluctuations Ellis:ELDSM:85 ;
Touchette:large_dev:09  – precisely the opposite of deterministic evolution –
is that of a large-deviation function.

The chain rule in Eq. (34) of Sec. III.2.2 relates the conditional relative
entropy associated with quasi-deterministic relaxation to an additional
unconditioned entropy of metastable states that are stable fixed points of the
Hamiltonian dynamical system. Here we complete the description of the relation
between the Lyapunov and large-deviation roles of the macrostate entropy (29),
and show how the sum over instantons rather than a single Hamiltonian
trajectory results in deterministic increase of the unconditioned relative
entropy $D\\!\left(1_{\underline{n}}\parallel\underline{\rho}\right)$. In the
Hamilton-Jacobi representation this means constructing relations between the
${\phi}^{\dagger}=1$ and the ${\phi}^{\dagger}\neq 1$ branches of the
$\mathcal{L}=0$ manifold. For special cases, such as CRNs with detailed
balance or one-dimensional systems, the mapping is one of simple time reversal
of the $n$ fields in trajectories. More generally, even for complex-balanced
CRNs where the Lyapunov and large-deviation functions are the same, the
relations between relaxation and escape trajectories become more variable.

####  IV.2.1 Convexity proof of the Lyapunov property of macrostate entropy in
Hamilton-Jacobi variables

Begin with the dynamics of the relative entropy ${\underline{S}}_{\rm
eff}\\!\left(\bar{n}\right)$ from Eq. (29), as the state variable $\bar{n}$
evolves along a relaxation solution to the Hamiltonian equations of motion.
The properties of ${\underline{S}}_{\rm eff}\\!\left(\bar{n}\right)$ follow
from its construction as a large-deviation function along escape trajectories.
From Eq. (20), the time derivative of ${\underline{S}}_{\rm eff}$ along an
escape trajectories is given by

$\displaystyle\frac{{\underline{S}}_{\rm eff}\\!\left(n\right)}{\partial
n}{\dot{n}}_{\rm esc}$
$\displaystyle={\theta}_{\underline{\rho}}\\!\left(n\right)\cdot{\dot{n}}_{\rm
esc}\equiv{\left.{\dot{\underline{S}}}_{\rm eff}\right|}_{\rm esc}.$ (48)

Hamiltonian trajectories are least-improbable paths of fluctuation, so escapes
are conditionally dependent along trajectories. The conditional probability to
extend a path, having reached any position along the path, is always positive,
giving ${\left.{\dot{\underline{S}}}_{\rm eff}\right|}_{\rm esc}\geq 0$ in Eq.
(48).

Next compute ${\dot{\underline{S}}}_{\rm eff}$ along a relaxation trajectory,
for simplicity considering $\mathbb{A}$ (an equivalent construction exists for
${\mathbb{A}}^{s}$ in terms of ${\bar{\rho}}^{s}$). The continuum limit of the
relative entropy from Eq. (4) replaces $\sum_{{\rm n}}\rightarrow\int dn$, and
continuously-indexed $-\log{\rho}^{\left(\bar{n}\right)}_{n}$ and
$-\log{\underline{\rho}}_{n}$ are defined through the large-deviation
functions.

Writing the CRN Liouville operator (35) in coherent-state arguments, the time
dependence is evaluated as

$\displaystyle\dot{D}\\!\left({\rho}^{\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
$\displaystyle\rightarrow\int
d^{s}\\!n\,\log\left(\frac{{\rho}^{\left(\bar{n}\right)}_{n}}{{\underline{\rho}}_{n}}\right){\dot{\rho}}^{\left(\bar{n}\right)}_{n}$
$\displaystyle\sim\int
d^{s}\\!n\,\log\left(\frac{{\rho}^{\left(\bar{n}\right)}_{n}}{{\underline{\rho}}_{n}}\right)\left[{\psi}^{T}\\!\left(z\right)\mathbb{A}\psi\\!\left(\frac{n}{z}\right)\right]{\rho}^{\left(\bar{n}\right)}_{n}$
$\displaystyle=\sum_{ji}k_{ji}\int
d^{s}\\!n\,\log\left(\frac{{\rho}^{\left(\bar{n}\right)}_{n}}{{\underline{\rho}}_{n}}\right)e^{\theta\left(Y_{j}-Y_{i}\right)}{\psi}_{i}\\!\left(n\right){\rho}^{\left(\bar{n}\right)}_{n}$
$\displaystyle\approx\sum_{ji}k_{ji}\int
d^{s}\\!n\,{\left.\frac{\partial}{\partial
n_{p}}\log\left(\frac{{\rho}^{\left(\bar{n}\right)}_{n}}{{\underline{\rho}}_{n}}\right)\right|}_{\bar{n}}{\left(n-\bar{n}\right)}_{p}$
$\displaystyle\mbox{}\times{\left(n-\bar{n}\right)}_{q}{\left.\frac{-{\partial}^{2}\log{\rho}^{\left(\bar{n}\right)}_{n}}{\partial
n_{q}\partial
n_{r}}\right|}_{\bar{n}}{\left(Y_{j}-Y_{i}\right)}_{r}{\psi}_{i}\\!\left(n\right){\rho}^{\left(\bar{n}\right)}_{n}$
$\displaystyle=\sum_{ji}k_{ji}{\left.\frac{-\partial\log{\underline{\rho}}_{n}}{\partial
n_{p}}\right|}_{\bar{n}}{\left(Y_{j}-Y_{i}\right)}_{p}{\psi}_{i}\\!\left(\bar{n}\right)$
$\displaystyle={\left.\frac{-\partial\log{\underline{\rho}}_{n}}{\partial
n}\right|}_{\bar{n}}\cdot Y\hat{\mathbb{A}}{\psi}_{i}\\!\left(\bar{n}\right)$
$\displaystyle={\theta}_{\underline{\rho}}\\!\left(\bar{n}\right)\cdot\dot{n}_{{\rho}^{\left(\bar{n}\right)}}\\!\left(0\right)$
$\displaystyle={\left.{\dot{\underline{S}}}_{\rm
eff}\\!\left(\bar{n}\right)\right|}_{\rm rel}.$ (49)

In Eq. (49) $z$ abbreviates
$z_{{\rho}^{\left(\bar{n}\right)}}\\!\left(n\right)\equiv\exp{\theta}_{{\rho}^{\left(\bar{n}\right)}}\\!\left(n\right)$
from Eq. (22), for the $S_{\rm eff}$ that is the continuum approximation of
$-\log{\rho}^{\left(\bar{n}\right)}$. The third through fifth lines expand
$\mathbb{A}$ explicitly in terms of rate constants $k_{ji}$ following Eq.
(35), to collocate all terms in $\theta\equiv\log z$ in the Liouville
operator. The fourth line expands
$\log\left({\rho}^{\left(\bar{n}\right)}/\underline{\rho}\right)$ and $\theta$
to linear order in $n-\bar{n}$ in neighborhoods of the saddle point $\bar{n}$
of ${\rho}^{\left(\bar{n}\right)}$. The matrix
$\left[-{\partial}^{2}\log\rho/\partial n_{q}\partial n_{r}\right]$ is the
inverse of the Fisher metric that is the variance
$\left[\left<{\left(n-\bar{n}\right)}_{p}{\left(n-\bar{n}\right)}_{q}\right>\right]$
Amari:inf_geom:01 , so the product of the two is just the identity
$\left[{\delta}_{pr}\right]$.

In the penultimate line of Eq. (49),
${\theta}_{\underline{\rho}}\\!\left(\bar{n}\right)$ is the value of
$\partial{\underline{S}}_{\rm eff}/\partial n$ for the _escape_ trajectory
passing through $\bar{n}$, and now
$\dot{n}_{{\rho}^{\left(\bar{n}\right)}}\\!\left(0\right)$ is the velocity
along the _relaxation_ trajectory rather than the Hamilton-Jacobi escape
solution at $\bar{n}$. So the net effect of the large-deviation approximation
on relative entropy has been to replace escape with relaxation velocity
vectors at a fixed value of $\theta$ Legendre dual to $\bar{n}$.

Lemma:
${\theta}_{\underline{\rho}}\\!\left(\bar{n}\right)\cdot\dot{n}_{{\rho}^{\left(\bar{n}\right)}}\\!\left(0\right)\leq
0$.

Proof: The proof follows from four observations:

  1. 1.

$\mathcal{L}=0$: As noted, for both escapes and relaxations,
$\mathcal{L}\\!\left(\theta,n\right)=0$.

  2. 2.

Convexity: Both the potential value ${\theta}_{\underline{\rho}}$ for the
escape trajectory, and the velocity $\dot{n}$ of the relaxation trajectory,
are evaluated at the same location
$\bar{n}=n_{{\rho}^{\left(\bar{n}\right)}}\\!\left(0\right)$. The Liouville
function
$-\mathcal{L}=\sum_{ji}k_{ji}\left(e^{\theta\left(Y_{j}-Y_{i}\right)}-1\right){\psi}_{i}\\!\left(n\right)$,
with all $k_{ji}>0$, is convex on the $s$-dimensional sub-manifold of fixed
$n$. $\mathcal{L}$ is bounded above at fixed $n$, and in order for cycles to
be possible, shift vectors $\left(Y_{j}-Y_{i}\right)$ giving positive
exponentials must exist for all directions of $\theta$ in the stoichiometric
subspace. Therefore $\mathcal{L}\rightarrow-\infty$ at large
$\left|\theta\right|$ in every direction, and the region $\mathcal{L}>0$ at
fixed $n$ is bounded. The boundary $\mathcal{L}\\!\left(\theta,n\right)=0$ at
fixed $n$ is likewise convex with respect to $\theta$ as affine coordinates,
and $\mathcal{L}>0$ is its interior.

  3. 3.

Chord: The vector $\left(\theta-0\right)$ is thus a chord spanning the
$\mathcal{L}=0$ submanifold of co-dimension 1 within the $s$-dimensional
manifold of fixed $n$.

  4. 4.

Outward-directedness: The equation of motion
$\dot{n}=-\partial\mathcal{L}/\partial\theta$ gives
$\dot{n}\\!\left(\theta,n\right)$ as the _outward_ normal function to the
surface $\mathcal{L}\\!\left(\theta,n\right)=0$. The outward normal at
$\theta=0$ is the classical relaxation trajectory. Every chord
$\left(\theta-0\right)$ of the surface lies in its interior, implying that
${\theta}_{\underline{\rho}}\\!\left(n\right)\cdot\dot{n}\\!\left(0,n\right)<0$
for any $n$, and thus
${\theta}_{\underline{\rho}}\\!\left(\bar{n}\right)\cdot\dot{n}_{{\rho}^{\left(\bar{n}\right)}}\\!\left(0\right)\leq
0$. $\blacksquare$

The conditional part of the relative entropy,
$-D\\!\left({\rho}^{\left(\bar{n}\right)}\parallel\underline{\rho}\right)$, is
thus shown to be monotone increasing along relaxation trajectories, which is
the Lyapunov role for the entropy state function familiar from classical
thermodynamics. That increase ends when $\bar{n}$ terminates in the trajectory
fixed point $\underline{n}\\!\left(\bar{n}\right)$.

####  IV.2.2 Instantons and the loss of Large-deviation accessibility from
first passages

The deterministic analysis of Eq. (49) is refined by the inclusion of
instanton trajectories through the following sequence of observations,
completing the discussion begun in Sec. III.2.2. Relevant trajectories are
shown in Fig. 1.

  1. 1.

The 2nd law as formulated in Eq. (4) is approximated in the large-deviation
limit not by a _single_ Hamiltonian trajectory, but by the _sum of all_
Hamiltonian trajectories, from an initial condition. Along a single
trajectory,
$-D\\!\left({\rho}^{\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
could increase or decrease.

  2. 2.

$-D\\!\left({\rho}^{\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
_increases_ everywhere in the submanifold $\theta=0$ of the manifold
$\mathcal{L}=0$, by Eq. (49). This is the classical increase of (relative)
entropy of Boltzmann and Gibbs.
$-D\\!\left({\rho}^{\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
_decreases_ everywhere in the submanifold $\theta\neq 0$ of the manifold
$\mathcal{L}=0$, by Eq. (48). This is the construction of the log-probability
for large deviations. These escape paths, however, simply lead to the
evaluations of ${\underline{S}}_{\rm
eff}\\!\left(n\right)\sim{\left.-\log{\underline{\rho}}_{{\rm
n}}\right|}_{{\rm n}\approx n}$, the stationary distribution.

  3. 3.

If a CRN has a single fixed point $\underline{n}$, there is a unique
$\theta=0$ trajectory from any starting $\bar{n}$ to it, and
$-D\\!\left({\rho}^{\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
increases deterministically by Eq. (49) along that single path. The black
trajectory with arrow converging exactly in a domain fixed point is such a
path in Fig. 1.

  4. 4.

If a CRN has multiple fixed points and instantons, all trajectories are
exponentially close to the exact $\theta=0$ trajectory before they enter a
small neighborhood around the terminus $\underline{n}\\!\left(\bar{n}\right)$
of the exact trajectory; that is: they give the appearance of being the black
deterministic trajectory in Fig. 1. The trajectory sum is correspondingly
close to the deterministic relaxation that increases the conditional entropy
${\underline{S}}_{\rm
eff}\\!\left(\underline{n}\\!\left(\bar{n}\right)\right)-{\underline{S}}_{\rm
eff}\\!\left(\bar{n}\right)$ in Eq. (34).

  5. 5.

On longer times, however, the infinite sum of formally distinct Hamiltonian
trajectories disperses into a sum over series of instantons making a random
walk among fixed points, with an integral for each passage over the possible
times at which the escape occurs. (See Coleman:AoS:85  Ch.7.) Such a sum is
shown as a tree of colored first-passage trajectories in Fig. 1. The “cross-
sectional” sum at a single observation time _over_ instantons distinguished by
their escaping times gives the same result as a longitudinal line integral
_along_ a single instanton between the start time and the observation. That
integral of $-\left(\partial\log{\rho}^{\left(\bar{n}\right)}_{n}/\partial
n\right)\dot{n}$ through a full passage (escape instanton + succeeding
relaxation) gives ${\underline{S}}_{\rm
eff}\\!\left(\underline{n}\right)-{\underline{S}}_{\rm
eff}\\!\left({\underline{n}}^{\prime}\right)=\log\left(w^{\rm
macro}_{{\underline{n}}^{\prime}\underline{n}}/w^{\rm
macro}_{\underline{n}{\underline{n}}^{\prime}}\right)$. The escape from fixed
point $\underline{n}$ to a saddle between $\underline{n}$ and a fixed point
${\underline{n}}^{\prime}$ in an adjacent basin, which we denote
$\bar{n}={\underline{n}}^{\prime}\\!\ddagger\\!\underline{n}$, is an integral
over Eq. (48), while the relaxation from the saddle $\bar{n}$ to the fixed
point ${\underline{n}}^{\prime}$ is an integral over Eq. (49). These are the
classical “entropy fluctuations” of stochastic thermodynamics.

  6. 6.

The contribution to the probability of a trajectory from each instanton comes
only from the $\theta\neq 0$ sub-manifold, and is given by $w^{\rm
macro}_{{\underline{n}}^{\prime}\underline{n}}\sim
e^{-\left[{\underline{S}}_{\rm
eff}\\!\left(\bar{n}\right)-{\underline{S}}_{\rm
eff}\\!\left(\underline{n}\right)\right]}\equiv e^{-\Delta{\underline{S}}_{\rm
eff}}$, just the leaving rate from the macrostate $1_{\underline{n}}$. The
result, upon coarse-graining to the macroscale (see Table 1 and the top
diagram in Fig. 1) where first-passages become instantaneous elementary
events, is a new stochastic process on discrete states corresponding to the
mesoscale Hamiltonian fixed points
$\left\\{\underline{n},{\underline{n}}^{\prime},\ldots\right\\}$. The coarse-
grained counterpart to
$D\\!\left(1_{\underline{n}\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
from Eq. (34) is the Lyapunov function reduced by a transition matrix
${\mathbb{T}}^{\rm macro}$ with matrix elements $w^{\rm
macro}_{{\underline{n}}^{\prime}\underline{n}}$. The coarse-graining and the
reduction of
$D\\!\left(1_{\underline{n}\left(\bar{n}\right)}\parallel\underline{\rho}\right)$
are described in detail for a 2-basin example in Smith:LDP_SEA:11 .

  7. 7.

The properties of ${\mathbb{T}}^{\rm macro}$ are exactly those we have assumed
for $\mathbb{T}$ as inputs to the mesoscale, completing the self-consistency
of our recursively-defined multi-level model universe.

##  V Cycle decompositions and non-decrease of intrinsic and extrinsic
relative entropies

Schnakenberg Schnakenberg:ME_graphs:76  introduced a method to solve for the
dissipation of a driven CRN in terms of steady cyclic currents within the
network, coupled to flows into or out of environmental buffers modeled as
chemostats Polettini:open_CNs_I:14 . A similar decomposition can be used here
to derive the positive-semidefiniteness of the relative entropy and
housekeeping terms in Eq. (11). The form of the CRN generator leads to a
decomposition of entropy changes into sums of densities convected around a
basis of cycles by the stationary-state currents. Positivity is proved by
convexity arguments similar to those for the standard fluctuation theorems
Speck:HS_heat_FT:05 ; Harris:fluct_thms:07 ; Esposito:fluct_theorems:10 , and
the cycle decomposition expresses the duality relations of those theorems in
terms of shifts forward or backward around cycles.

Sec. IV constructed the relation between the hypergraph adjacency matrix
$\mathbb{A}$ and stoichiometry $Y$, and the microstate transition matrix
$\mathbb{T}$. Each reaction in the hypergraph acts regularly as a rule on
indefinitely many microstates, making CRNs an example of rule-based systems
Danos:rule_based_modeling:08 . This section extends that mapping to cycle
decompositions. For finitely-generated CRNs, there is a finite basis of either
cyclic flows through complexes or stoichiometrically coupled cycles called
hyperflows Andersen:generic_strat:14 , with the interpretation of mass-action
currents. From these, a (generally infinite) basis of finite-length cycles can
be constructed for flows in the microstate space, which projects onto the
basis in the hypergraph. We first use the finite basis of macrostate currents
to prove monotonicity results for $f$-divergences where those exist, showing
how the projection of cycles can result in a projection of the Lyapunov
property between levels.

The parallelism between microstates and macrostates only holds, however, when
the bases for currents at both levels contain only _cycles_. More generally,
when macrostate currents include non-cyclic hyperflows, the analogy between
the two levels is broken, leading to loss of the $f$-divergence as a Lyapunov
function on macrostate variables and generally to much greater difficulty of
solving for the stationary distribution on microstates. The conditions for
different classes of flows on the hypergraph are the basis for a complexity
classification of CRNs, and both the form of the large-deviation function and
the symmetry between its Lyapunov and large-deviation roles differ across
classes.

###  V.1 Complexity classes and cycle decompositions of stationary currents on
the hypergraph

Remarkably, the properties of classical fixed points in the Hamilton-Jacobi
representation are sufficient to classify CRNs in a hierarchy with three
nested levels of complexity, according to the mass-action currents at the
fixed points. Let $\underline{n}$ again denote the fixed point.

  1. 1.

The CRNs with detailed balance are those in which
$k_{ji}{\psi}_{i}\\!\left(\underline{n}\right)=k_{ij}{\psi}_{j}\\!\left(\underline{n}\right)$,
for all pairs of complexes $\left<j,i\right>$ connected by a reaction. Under
the descaling (43), this condition is that
${\hat{\mathbb{A}}}^{T}=\hat{\mathbb{A}}$.

  2. 2.

The CRNs with complex balance only require
$\sum_{i}k_{ji}{\psi}_{i}\\!\left(\underline{n}\right)=\sum_{j}k_{ij}{\psi}_{j}\\!\left(\underline{n}\right)$,
for each complex $j$, or under descaling,
$\hat{\mathbb{A}}\psi\\!\left(1\right)=0$.

  3. 3.

The general case requires only $Y\mathbb{A}\psi\\!\left(\underline{n}\right)$,
the condition of current balance at each species $p$, or under descaling, only
$Y\hat{\mathbb{A}}\psi\\!\left(1\right)=0$.

The classes are summarized in Table. 3,

case | complexity class
---|---
$Y\hat{\mathbb{A}}1=0$ | general case
$\hat{\mathbb{A}}1=0$ | complex-balanced
${\hat{\mathbb{A}}}^{T}=\hat{\mathbb{A}}$ | detailed-balanced
Table 3:  Hierarchical categories of stationary distributions of CRNs.

####  V.1.1 Complex balance and relations of the Lyapunov and large-deviation
roles of $S_{\rm eff}$

It is known that complex-balanced CRNs (with technical conditions to ensure
ergodicity on the complex network under $\mathbb{A}$
Gunawardena:CRN_for_bio:03 ; Baez:QTRN:13 ) possess unique interior fixed
points Feinberg:notes:79 , and moreover that their _exact_ steady states
$\underline{\rho}$ are products of Poisson distributions or sections through
such products Anderson:product_dist:10  defined by conserved quantities under
the stoichiometry. Their generating functions are the coherent states from
Sec. III.2.1. It is immediate that all classical states obtained by
exponential tilts with parameter $z=e^{\theta}$ are likewise coherent states
with state variables

$\displaystyle{\bar{n}}_{p}$ $\displaystyle\equiv
e^{{\theta}_{p}}{\underline{n}}_{p},$ (50)

and that the $f$-divergence of Eq. (30) is the macrostate relative entropy
function.

From the equations of motion (45), within either the ${\phi}^{\dagger}\equiv
1$ or ${\phi}^{\dagger}\neq 1$ sub-manifolds, the time derivatives of
${\underline{S}}_{\rm eff}$ along escape and relaxation trajectories, Eq. (49)
and Eq. (48) respectively, take the forms

$\displaystyle{\left.{\dot{\underline{S}}}_{\rm eff}\right|}_{\rm rel}$
$\displaystyle=\log\left(\frac{n}{\underline{n}}\right)\cdot{\dot{n}}_{\rm
rel}=\log\left(\frac{n}{\underline{n}}\right)\cdot
Y\hat{\mathbb{A}}\psi\\!\left(\frac{n}{\underline{n}}\right)$
$\displaystyle{\left.{\dot{\underline{S}}}_{\rm eff}\right|}_{\rm esc}$
$\displaystyle=\log\left(\frac{n}{\underline{n}}\right)\cdot{\dot{n}}_{\rm
esc}=-\log\left(\frac{n}{\underline{n}}\right)\cdot
Y{\hat{\mathbb{A}}}^{T}\psi\\!\left(\frac{n}{\underline{n}}\right).$ (51)

By Eq. (50), ${\psi}_{i}$ appearing in Eq. (51) evaluate to

$\displaystyle{\psi}_{i}\\!\left(\frac{n}{\underline{n}}\right)$
$\displaystyle=e^{\theta Y_{i}}.$ (52)

##### Finite cycle decomposition of the steady state current:

By definition of complex balance, $\hat{\mathbb{A}}1=0$, in addition to the
general condition that $1^{T}\hat{\mathbb{A}}=0^{T}$. Any such matrix acting
over finitely many complexes may be written as a sum of adjacency matrices for
cyclic flows, which we index $\alpha$, with positive coefficients
${\bar{\jmath}}_{\alpha}$ equal to the currents around those cycles in the
stationary state. For the subclass with detailed balance, a decomposition in
cycles of length 2 is always possible.

Letting $\sum_{ji\mid\alpha}$ denote the sum over directed links $ji$ in order
around the cycle $\alpha$, the trajectory derivatives (51) may be decomposed
as

$\displaystyle{\left.{\dot{\underline{S}}}_{\rm eff}\right|}_{\rm
rel}=\theta\cdot{\dot{n}}_{\rm rel}$
$\displaystyle=-\sum_{\alpha}{\bar{\jmath}}_{\alpha}\sum_{ji\mid\alpha}e^{\theta
Y_{i}}\cdot\log\left(\frac{e^{\theta Y_{i}}}{e^{\theta Y_{j}}}\right)$
$\displaystyle{\left.{\dot{\underline{S}}}_{\rm eff}\right|}_{\rm
esc}=\theta\cdot{\dot{n}}_{\rm esc}$
$\displaystyle=\sum_{\alpha}{\bar{\jmath}}_{\alpha}\sum_{ji\mid\alpha}e^{\theta
Y_{j}}\cdot\log\left(\frac{e^{\theta Y_{j}}}{e^{\theta Y_{i}}}\right).$ (53)

($\cdot$ indicates the vector inner product over the species index $p$.)

Letting $\sum_{i\mid\alpha}$ denote the sum over all complexes in the cycle
$\alpha$, the complex activities (52) may be normalized to vectors
$p_{\alpha}\equiv\left[p_{i\mid\alpha}\right]$ with unit measure, as

$p_{i\mid\alpha}\equiv\frac{e^{\theta Y_{i}}}{\sum_{j\mid\alpha}e^{\theta
Y_{j}}}.$ (54)

Then the trajectory derivatives (53), themselves the time derivatives of the
$f$-divergence (34), may be written in terms of positive semidefinite KL
divergences of the measures $p_{\alpha}$ from their own images advanced or
retarded by one complex around each cycle:

$\displaystyle\theta\cdot{\left.\dot{n}\right|}_{\rm rel}$
$\displaystyle=-\sum_{\alpha}\left({\bar{\jmath}}_{\alpha}\sum_{i\mid\alpha}e^{\theta
Y_{i}}\right)\sum_{ji\mid\alpha}p_{i\mid\alpha}\cdot\log\left(\frac{p_{i\mid\alpha}}{p_{j\mid\alpha}}\right)$
$\displaystyle\theta\cdot{\left.\dot{n}\right|}_{\rm esc}$
$\displaystyle=\sum_{\alpha}\left({\bar{\jmath}}_{\alpha}\sum_{i\mid\alpha}e^{\theta
Y_{i}}\right)\sum_{ji\mid\alpha}p_{j\mid\alpha}\cdot\log\left(\frac{p_{j\mid\alpha}}{p_{i\mid\alpha}}\right).$
(55)

Non-negativity of the KL divergences in every term of Eq. (55) recovers the
monotonicities of ${\underline{S}}_{\rm eff}$ along relaxation and escape
trajectories from Sec. IV.2.1. The total time derivatives are decomposed into
terms over finitely many cycles, each independently having the same sign, a
stronger decomposition than can be obtained from the convexity proofs for
increase of relative entropy Schnakenberg:ME_graphs:76  for general CRNs.

##### Locally linear coordinates for classical entropy change in complex-
balanced CRNs

Note that the cycle-KL divergences, multiplied by the density factors
$\sum_{i\mid\alpha}e^{\theta Y_{i}}$, define a coordinate system that is
locally invertible for the coordinates $\theta$ except at isolated points:

$\displaystyle x_{\alpha-}$ $\displaystyle\equiv\sum_{ji\mid\alpha}e^{\theta
Y_{i}}\log\left(\frac{e^{\theta Y_{i}}}{e^{\theta Y_{j}}}\right)$
$\displaystyle x_{\alpha+}$ $\displaystyle\equiv\sum_{ji\mid\alpha}e^{\theta
Y_{j}}\log\left(\frac{e^{\theta Y_{j}}}{e^{\theta Y_{i}}}\right).$ (56)

The time derivatives of ${\underline{S}}_{\rm eff}$ from Eq. (53) may be
written in these coordinates as a linear system,

$\displaystyle{\left.{\dot{\underline{S}}}_{\rm eff}\right|}_{\rm
rel}={\dot{D}\\!\left(e^{\theta Y}\parallel 1\right)}_{\rm
rel}=\theta\cdot{\dot{n}}_{\rm rel}$
$\displaystyle=-\sum_{\alpha}{\bar{\jmath}}_{\alpha}x_{\alpha-}$
$\displaystyle{\left.{\dot{\underline{S}}}_{\rm eff}\right|}_{\rm
esc}={\dot{D}\\!\left(e^{\theta Y}\parallel 1\right)}_{\rm
esc}=\theta\cdot{\dot{n}}_{\rm esc}$
$\displaystyle=\sum_{\alpha}{\bar{\jmath}}_{\alpha}x_{\alpha+}.$ (57)

Fig. 3 shows examples of the coordinates (56) for a 2-cycle and a 3-cycle in
the simplex of normalized activities ${\psi}_{i}/\sum_{j}{\psi}_{j}$ for three
species $i\in\left\\{A,B,C\right\\}$.

Figure 3:  Simplex for a probability distribution
$p=\left(p_{A},p_{B},p_{C}\right)$, with two coordinates $x_{AB}$ and
$x_{ABC}$ of the form in Eq. (56) on a 2-cycle and a 3-cycle.

####  V.1.2 Vorticity in the flowfield of stationary trajectories

Recall (Table 3) that detailed balance means
${\hat{\mathbb{A}}}^{T}=\hat{\mathbb{A}}$. In this case Eq. (51) implies equal
and opposite change of ${\underline{S}}_{\rm eff}$ along relaxations and
escapes, because the two trajectories are time-reverses of each other.
Detailed balance thus generalizes the time-reversal symmetry of 1-dimensional
systems to any number of dimensions, and is a correspondingly restrictive
condition. This is the assumption in classical thermodynamics that identifies
the Lyapunov and large-deviation roles of the entropy state function.

Already for complex-balanced CRNs, exact time reversal is generally broken.
The gradients of the tangent vectors to relaxation and escape trajectories,
with respect to the exponential-family coordinates $\theta$ in Eq. (50),
evaluate to

$\displaystyle\frac{\partial{\left({\dot{n}}_{\rm
rel}\right)}_{p}}{\partial{\theta}_{q}}$
$\displaystyle=Y_{p}\hat{\mathbb{A}}\operatorname{diag}\left[\psi\\!\left(e^{\theta}\right)\right]Y^{T}_{q}$
$\displaystyle\frac{\partial{\left({\dot{n}}_{\rm
esc}\right)}_{p}}{\partial{\theta}_{q}}$
$\displaystyle=-Y_{p}{\hat{\mathbb{A}}}^{T}\operatorname{diag}\left[\psi\\!\left(e^{\theta}\right)\right]Y^{T}_{q}.$
(58)

A measure of the time-asymmetry of the CRN is the vorticity, defined as the
antisymmetric part of the matrices (58). This vorticity equals the gradient of
the _sum_ of tangent vectors to relaxations and escapes at the same point,

$\displaystyle\frac{\partial{\left({\dot{n}}_{\rm esc}+{\dot{n}}_{\rm
rel}\right)}_{p}}{\partial{\theta}_{q}}$
$\displaystyle=Y_{p}\left(\hat{\mathbb{A}}-{\hat{\mathbb{A}}}^{T}\right)\operatorname{diag}\left[\psi\\!\left(e^{\theta}\right)\right]Y^{T}_{q}.$
(59)

At the fixed point where $\theta=0$, the vorticity is the antisymmetric part
of the matrix $YAY^{T}$.

####  V.1.3 Hyperflow decomposition for non-complex-balanced CRNs

For complex-balanced CRNs, $\hat{\mathbb{A}}$ describes flows on the complex
network, which is an ordinary directed graph. For non-complex balanced flows,
$Y\hat{\mathbb{A}}1=0$ but ${\mathbb{A}}1\neq 0$. $\hat{\mathbb{A}}$ therefore
cannot be written as a sum of adjacency matrices for currents on cycles in the
complex graph. However, because net current still equals zero for every
species, a basis for the currents at any fixed point still exists in balanced
hyperflows. If the elements of $Y$ are all integer values (as they are for
chemistry or for biological population processes), the basis elements are
balanced integer hyperflows Andersen:generic_strat:14 . Each such flow is an
assignment of a set of non-negative integers to every current with nonzero
$k_{ji}$, such that the net currents at each species vanish.

We may continue to refer to basis hyperflows with index $\alpha$ by extension
of the complex-balanced case, and for each basis element there is a non-
negative current ${\bar{\jmath}}_{\alpha}$ which is the coefficient of that
flow in the stationary-state current solution. Integer hyperflows cannot be
used to express ${\dot{\underline{S}}}_{\rm eff}$ as sums of KL divergences
with uniform sign, in keeping with the variety of complications that arise for
${\dot{\underline{S}}}_{\rm eff}$ in for non-complex-balanced CRNs.

###  V.2 Cycle decomposition in the microstate space, and non-decrease of
relative entropy components

Next we apply a cycle decomposition similar to the one in the previous
section, but in the microstate space rather than on the hypergraph, to prove
non-negativity of the first two terms in Eq. (11). The crucial distinction
between the two levels is that balanced integer flows can always be mapped to
cycles in the state space. For cycles in the hypergraph a natural map is
unique; for more general flows many mappings may be possible. Relative
entropies thus give Lyapunov functions more generally for distributions over
microstates than for macrostate variables. However, unlike the hypergraph
where basis decompositions are _finite_ , in the state space they generally
require solution for infinitely many cycle currents and are difficult except
in special cases.

####  V.2.1 The system-marginal relative entropy from ${\bar{\rho}}^{s}$

The system $s$ and environment $e$ must now be considered explicitly because
housekeeping heat is an embedding of $s$ in $s\otimes e$. We therefore work in
system indices $\left\\{{{\rm n}}_{s}\right\\}$ and descale with the steady
state ${\bar{\rho}}^{s}$ at the instantaneous ${\rho}^{e\mid s}$. In the
following derivations ${\rho}^{s}$ and ${\rho}^{e\mid s}$ can be any
normalizable distributions; macrostates are no longer singled out. The time-
change of $D\\!\left({\rho}^{s}\parallel{\bar{\rho}}^{s}\right)$, fixing