Evolution of complexity and the origins of biochemical life

Praful Gagrani

Wisconsin Institute for Discovery,

University of Wisconsin-Madison,

Madison, WI 53715, USA

While modern physics and biology satisfactorily explain the passage from the Big Bang to the formation of Earth and the first cells to present-day life, respectively, the origins of biochemical life still remain an open question. Since life, as we know it, requires extremely long genetic polymers, any answer to the question must explain how an evolving system of polymers of ever-increasing length could come about on a planet that otherwise consisted only of small molecular building blocks. In this work, we show that, under realistic constraints, an abstract polymer model can exhibit dynamics such that attractors in the polymer population space with a higher average polymer length are also more probable.
We generalize from the model and formalize the notions of complexity and evolution for chemical reaction networks with multiple attractors.
The complexity of a species is defined as the minimum number of reactions needed to produce it from a set of building blocks, which in turn is used to define a measure of complexity for an attractor. A transition between attractors is considered to be a progressive evolution if the attractor with the higher probability also has a higher complexity. In an environment where only monomers are readily available, the attractor with a higher average polymer length is more complex. Thus, by this criterion, our abstract polymer model can exhibit progressive evolution for a range of thermodynamically plausible rate constants. We also formalize criteria for open-ended and historically-contingent evolution and explain the role of autocatalysis in obtaining them. Our work provides a basis for searching for prebiotically plausible scenarios in which long polymers can emerge and yield populations with even longer polymers. Additionally, the existence of features like history-dependence and open-endedness support the view that the path of complexification from chemistry to biology, was one of gradual complexification rather than an instantaneous origin of life.

§ INTRODUCTION

Phase transitions are the changes of a system from one regime to another exhibiting distinctly dissimilar properties, and the phase transition paradigm finds applications in a wide variety of fields such as cosmology [1], condensed matter physics [2], and evolutionary games [3]. Theorists have often viewed the origin-of-life problem through the lens of phase transitions [4, 5, 6]. In 1982, Freeman Dyson proposed an abstract model [7] for the origin of life that allowed a large population of molecular units to make the transition from a disordered phase of no-life to an ordered phase of life with a calculable probability. He concludes with a discussion of how the model could be used to explain the emergence of the modern genetic apparatus [8] in life. Smith and Morowitz [9], in contrast, argue that a genetic apparatus could not be the result of a single phase transition, but required a cascade of non-equilibrium phase transitions in order to open up the required chemical energy flows on the early Earth.

Arguably the most remarkable feature of cellular life is the capacity for evolving lineages to repeatedly transition to new phases (genotypes), each of which is characterized by a set of long polymers (DNA, RNA and/or proteins) that differ in sequence and population distribution from previously occupied phases. The production of specific long polymers in itself is noteworthy, since it is assumed that the environment provides monomers (or even smaller chemical building blocks) and hydrolysis (fragmentation reactions) typically have higher rate constants than condensation (polymerization) reactions. This is a reasonable assumption since, otherwise, we would be in a regime in which run-away chemistry would tend to yield arbitrary mixtures of very high molecular mass compounds. Under these assumptions, the current regime in which biology uses specific, extremely long polymers must be the result of a set of transitions to attractors enriched for longer and longer polymers. Even more remarkably, in polymer-based (genetic) evolution there seems to be no limit to the number of phases and the phases accessible in the future are highly sensitive to the current phase. Explaining this open-ended and historically-contingent evolution is a major challenge in origin of life research. How could simple chemistry enter a regime where long polymers are likely, where there are indefinitely many alternative phases, where viable phases can be visited in an open-ended series, and where historical contingency results in independent lineages' diverging over time?

The non-equilibrium potential (NEP) against the average polymer size for an abstract one-monomer polymer model is shown. The minima of the NEP correspond to attractors of the system, and the depth of the potential well around the attractor signifies the stability of the attractor. Notice that, for this model, there are several attractors with various average polymer sizes, and the attractors with a higher average size are also more stable. A physically realizable model with this property will imply that a system starting close to the origin will keep ratcheting itself towards attractors with higher average sizes, explaining how the preconditions to modern genetic apparatus could arise.

In this work, we introduce a phase transition formalism for stochastic chemical reaction networks (CRNs) and mathematically define key terms of relevance to evolutionary biology. We introduce concepts from dynamical systems, algorithmic complexity, and probability theory to provide rigorous criteria for when a stochastic CRN is capable of progressive, historically-contingent, or open-ended evolution. Moreover, we construct thermodynamically realistic simple polymerization CRN models where a sequence of phases with progressively higher average polymer lengths also show monotonically increasing stability (see Fig. <ref>). This, and related models are developed in the examples throughout the text. Along with an exemplification of our formalism, they constitute a mathematical demonstration of conditions modelled as a stochastic CRN that permit the origins of biochemical life.

The organization of the paper is as follows. In Sec. <ref>, we introduce the mathematical framework of the phase transition paradigm for stochastic CRNs and provide examples of polymerization models relevant to the origins of biochemical life. In Sec. <ref>, we remark on the role of entropy production and autocatalysis in our framework and examples. Finally, in Sec. <ref>, we summarize our contributions and elaborate on avenues for future research.

§ MATHEMATICAL FORMALISM AND EXAMPLES

The $a.$ attractor-transition graph, $b.$ complexity graph, and $c.$ evolution graph for the 1-monomer polymerization model $\mcl{P}_{1,4}$ from e.g. <ref>. Since attractor $a_h$ is both more stable and more complex than $a_l$, the system will be said to be capable of progressive evolution from $a_l$ (Sec. <ref>).

In this section, we formalize the notion of evolution for stochastic CRNs. Starting from a thermodynamically consistent CRN, we explain how to decompose it into system-environment partitions, detect how many attractors the system can have, and calculate the relative probabilities of attractors in the stationary distribution. Then, introducing tools from algorithmic complexity, we explain how to assign a measure of complexity to the different attractors, and finally specify criteria for when a CRN can be said to be capable of exhibiting evolution. As an application of the formalism to the origins of biochemical life, we demonstrate a class of abstract polymer models that are capable of evolving to attractors with higher average lengths. A graphical representation of a model with two attractors, that is capable of progressively evolving from the attractor with the lower average length to the higher one, is shown in Fig. <ref>.

§.§ Internal-boundary-external partition

A tripartite decomposition of a complete CRN into internal-internal, boundary, and external reaction sets of a 1-monomer polymerization model from e.g. <ref> with $g=1$ and $N=4$, which we denote as $\mcl{P}'_{1,4}$. An analysis with rate constant assignments of $\mcl{P}_{1,4}$ can be found in e.g. <ref>.

A CRN $\mcl{G}$ taken under mass-action kinetics (MAK) will be given as a triple
$$\mcl{G} = \{\Spc,\R,\K\},$$

where $\Spc,\R$ and $\K$ are the sets of species, reactions, and rate constants, respectively (the set of complexes is omitted as it can be inferred from the reaction set [10]). A brief review of CRNs, stochastic CRNs, and thermodynamics of CRNs is given in appendices <ref>, <ref>, and <ref>, respectively.

We refer to the CRN of all the relevant species and reactions in the system and environment taken under MAK as the complete CRN and denote it as

$$\mcl{G}' = \{\Spc',\R',\K'\}.$$

For thermodynamic consistency, as explained in App. <ref>, the rate constants should be such that the complete CRN has a detailed-balanced equilibrium. This assumption implies that, under the absence of any external control, the concentrations of all species will relax to their equilibrium values such that the CRN produces no more entropy. As a corollary, it is implied that the CRN is reversible, and henceforth we use $\R$ to refer to the set of reactions unique up to reversal (see Eq. <ref>).

In order to permit evolution-like behavior, we will assume that system is open, meaning that a subset of species, called external species and denoted by $\mcl{E}$, are controlled by the environment. The set of remaining species, that follow the assumptions made by MAK, are called internal species and denoted as $\Spc$, yielding

\[ \Spc' = \Spc \cup \phantom{\cdot} \mcl{E}. \]

This partition on the species set induces an internal-external partition of the reactions set $\R'$ into the set of all reactions in which any internal species is consumed or produced $\R$ and its complement where no internal species participates

\[ \R' = \R \cup \notR.\]

Denoting the restriction of reactions in $\R$ to the species set $\Spc$ as $\R_{\Spc}$, we define:

\begin{align*}
\textbf{Internal CRN  } & \mcl{G}:= \{\Spc,\R_{\Spc},\K\},\\
\textbf{External CRN  } & \mcl{G}_e:= \{\mcl{E},\notR,\K'_e\},
\end{align*}

where $\K$ and $\K'_e$ are modified rate constants obtained from $\K'$ (see Eq. <ref>).
The stoichiometric matrix of the complete CRN, denoted as $\nabla$, then takes the form

\begin{equation}
\nabla = \qquad \bordermatrix{
~  & \tikzmark{harrowleft} \R &  & \notR & \tikzmark{harrowright}  \cr
\tikzmark{varrowtop}
\mbox{\normalfont\Large $\mcl{E}$}  & \mbox{\normalfont\Large $\nabla^{\R}_\mathcal{E}$} &  \begin{matrix}
\phantom{a} \\
\phantom{c}
\end{matrix} \rvline &
\mbox{\normalfont\Large $\nabla^{\notR}_\mathcal{E}$}\cr
\hline \cr
\mbox{\normalfont\Large $\Spc$}
& \mbox{\normalfont\Large $\St$} & \begin{matrix}
\phantom{a} \\
\phantom{c}
\end{matrix}\rvline & \bigzero \cr
\tikzmark{varrowbottom} &  &  &
\tikz[overlay,remember picture] {
\draw[->] ([yshift=3ex]harrowleft) -- ([yshift=3ex]harrowright)
node[midway,above] {\scriptsize Reactions};
\draw[->] ([yshift=1.5ex,xshift=-2ex]varrowtop) -- ([xshift=-2ex]varrowbottom)
node[sloped, midway,below] {\scriptsize Species};
\label{eq:full_stoich}
\end{equation}

where $\St$ is the stoichiometric matrix of the internal CRN $\mcl{G}$.

While, in the absence of control, the complete CRN would relax to the detailed-balanced equilibrium, in the presence of control by the environment, the internal CRN does not necessarily do the same. Typically, the system is driven towards and maintained at a non-detailed-balanced equilibrium, called non-equilibrium steady state (NESS), which is characterized by nonzero fluxes and nonzero chemical gradients [11]. Moreover, in contrast to a detailed-balanced equilibrium, the entropy production rate (EPR) (Eq. <ref>) at a NESS is non-zero and equals the amount of thermodynamic work done by the chemical potential energy of the external species to maintain the NESS (see App. <ref>).

To facilitate investigation of NESSs, we further partition the internal reaction set $\R$ into internal-internal and boundary reaction sets $\R_i$ and $\R_b$, respectively. Furthermore, we partition the internal species set $\Spc$ into the internal-boundary set $\Spc_b$ and internal-internal set $\Spc_i$. Similarly, we partition the external species $\mcl{E}$ into the the external boundary, $\mcl{E}_b$, and the external-external set, $\mcl{E}_e$, yielding

\begin{align*}
\Spc &= \Spc_b \cup \Spc_i,\\
\mcl{E} &= \mcl{E}_b \cup \mcl{E}_e.
\end{align*}

This partition decomposes the CRN $\mcl{G}'$ into three subnetworks (as illustrated in Fig. <ref>):

\begin{align*}
\textbf{Internal-internal CRN  } & \mcl{G}'_i = \{\Spc_i \cup \Spc_b, \R_i,\K'_i\},\\
\textbf{Boundary CRN  } & \mcl{G}'_b = \{\Spc_b \cup \phantom{\cdot} \mcl{E}_b, \R_b,\K'_b\},\\
\textbf{External-external CRN  } & \mcl{G}'_e = \{\mcl{E}_b \cup \mcl{E}_e, \notR,\K'_e\},
\end{align*}

where $\K'_x$ is the restriction of the rate constants of $\mcl{G}'$ to $\mcl{G}'_x$ (explained below), and the stoichiometric matrix takes the form:

\begin{equation}
\nabla = \quad \bordermatrix{
~  &  \R_i & \R_b  & & \notR   \cr            \begin{matrix}
\largeentry{\mcl{E}_e} \\
\largeentry{\mcl{E}_b}
\end{matrix}  & \bigzero  & \begin{matrix}
\mbf{0} \\
\nabla^{\R_b}_{\mcl{E}_b}
\end{matrix} &
\rvline &
\Largeentry{\nabla^{\notR}_{\mcl{E}}}
\cr
\hline \cr
\begin{matrix}
\largeentry{\mcl{S}_b} \\
\largeentry{\mcl{S}_i}
\end{matrix}  & \Largeentry{\St_i}  & \begin{matrix}
\St_b \\
\mbf{0}
\end{matrix} &
\rvline &
\bigzero \cr
}. \label{eq:stoich_w_boundary}
\end{equation}

Internal-internal CRN:

The internal-internal CRN $\mcl{G}'_i$ consists only of reactions that do not involve any external species. For relevance to origins of life models, $\mcl{G}'_i$ can be polymerization models [12, 13], or cluster CRNs (CCRNs) (Sec. 5, [14]). In CCRNs, the species are counts of the constituents, forgetting bond structure, and reactions are generated by rules that satisfy certain conservation laws on the constituents.

For example, consider polymers constituted of two monomer types $A$ and $B$ such as $AAB, ABA, AAAB,$ etc. Following notation from [14], a cluster species in a CCRN with $m$ constituents will be denoted as $\ovl{\mbf{n}} = \ovl{n_1,\ldots,n_m}$. Thus, continuing the example, both polymers $AAB$ and $ABA$ map to the cluster $\ovl{2,1}$, $AAAB$ maps to $\ovl{3,1}$, and so on (where we have arbitrarily chosen an ordering, $A,B$). The advantage of using CCRNs over polymerization models is that the size of the species set is polynomial in a CCRN rather than exponential, as it is for polymers. Finally, the reaction set in a CCRN is algorithmically computable and will only contain reactions that satisfy integral conservation laws, of the type

\begin{align*}
\ce{ $\ovl{i,j} + \ovl{k,l}$ &<=> $\ovl{i\text{+}j,k\text{+}l}$},\\
\ce{ $\ovl{i,j} + \ovl{k,l}$ &<=> $\ovl{m,n} + \ovl{i\text{+}k\text{-}m,j\text{+}l\text{-}n}$},
\end{align*}

and so on, where $i,j,k,l, i\text{+}k\text{-}m, j\text{+}l\text{-}n \in \mathbb{Z}_{\geq 0}$.

External-external CRN:
For a complete CRN $\mcl{G}'$, there are several ways to open the system. One could fix the rate of change in concentration of boundary species, or employ more complicated control mechanisms, such as where the external CRN controls the concentration of boundary species through feedback [15]. The kind of control systems we consider here are those where the external-boundary species $\mcl{E}_b$ are chemostatted or fixed at a constant concentration by the environment. Thus, henceforth, we assume that the external-external species set $\mcl{E}_e$ and reaction set $\notR$ are empty.

Boundary CRN:

Let us assume that the external-boundary species are chemostatted at $q^*$. The rate constants $k_r \in \K_b$ on the boundary reactions $r \in \R_b$ are then obtained using $k'_r \in \K'_b$ and the chemostatted concentration such that

\begin{equation}
k_r = k_r' \prod_{s \in \mcl{E}_b} (q^*_s)^{r^-_s}. \label{eq:rate_const_renorm}
\end{equation}

Notice that the re-assignment of rate constants through this procedure leads to regimes where the Wegscheider's conditions [16] are not satisfied leading to non-detailed balanced equilibrium. This consideration is crucial to mulstistability and we return to this point in Sec. <ref>.

The internal-boundary and external-boundary species sets play different roles in an open-system.
As derived in Eq. <ref>, at a NESS, the EPR is accounted for only using the chemical potential of the external-boundary species set $\mcl{E}_b$ [17]. On the other hand, through their interaction with the external-boundary species, the internal-boundary species may break conservation laws in the internal CRN that are otherwise present in the internal-internal CRN (see Eg. <ref>).

Using the above decomposition, the internal CRN is

\begin{align*}
\mcl{G} &= \{ \Spc, \R,\K\} \\
&= \{ \Spc_i \cup \Spc_b, \R_i \cup (\R_b)_i, \K'_i \cup \K_b\},
\end{align*}

and its stoichiometric matrix is given by

\[ \St =
\begin{pmatrix}
\Largeentry{\St_i}  & \begin{matrix}
\St_b \\
\mbf{0}
\end{matrix}
\end{pmatrix}.
\]

By construction, the internal CRN can exhibit several NESSs and their entropy production is informed by thermodynamic considerations. However, as explained in Sec. <ref>, simple thermodynamic considerations of entropy production and free-energy transduction are insufficient to formalize evolution.

We define a class of models with relevance to the origins of biochemical life that we will use to illustrate our formalism, namely 1-monomer polymerization models. The class is parameterized by two positive integers $g$ and $N$, and we denote a complete CRN and internal CRN in the class as $\mcl{P}'_{g,N}$ and $\mcl{P}_{g,N}$, respectively. $N$ and $g$ count the number of polymers and functional species (see Sec. <ref>), respectively, in the internal CRN.

It can be shown that the complete CRN has deficiency zero (see App. <ref> for definition and e.g. <ref> for a similar example). Since it is also reversible, by Feinberg's deficiency zero theorem (Thm. <ref>), the complete CRN has a unique equilibrium in each stoichiometric compatibility class for any choice of rate constants. On the other hand, the internal CRN has deficiency non-zero, and it will be shown in subsequent examples that there do exist rate constants such that the internal CRN exhibits multistability.

The internal-internal CRN, denoted by $(\mcl{P}'_i)_{N}$ consists of the set of reactions representing polymerization by monomer addition at an end (since it does not explicitly depend on $g$, we suppress that parameter). Since a 1-monomer polymer model is indistinguishable from a 1-constituent CCRN, we employ the cluster notation to define the CRN.

\begin{align*}
(\mathcal{P}'_i)_N &= \{ \Spc_N, \R_N, \K_N\}\\
\Spc_N &= \{ \ovl{1}, \ovl{2}, \ldots, \ovl{N}\},\\
\R_N &= \{ \ce{$\ovl{1} + \ovl{1}$ <=>[$k_{2}^+$][$k_{2}^-$] $\ovl{2}$},\\
& \phantom{=2} \ce{$\ovl{1} + \ovl{2}$ <=>[$k_{3}^+$][$k_{3}^-$] $\ovl{3}$},\\
& \phantom{=2} \ldots \\
& \phantom{=2} \ce{$\ovl{1} + \ovl{N\text{-1}}$ <=>[$k_{N}^+$][$k_{N}^-$] $\ovl{N}$}\}.
\end{align*}

Employing notation from Sec. 5 of [14], we can condense the reaction set as the following rule set

\[\R_N = \{ \ce{$\ovl{1} + \ovl{j\text{-}1}$ <=>[$k_{j}^+$][$k_{j}^+$] $\ovl{j}$} \phantom{\cdot}\big| \phantom{\cdot} 1 \leq j \leq N \in \mathbb{Z}\}.\]

The stoichiometric matrix $\St_i$ for the internal-internal CRN becomes

\begin{equation*}
\St_i = \bordermatrix{
~  &   & &   &     \cr
\ovl{1} &-1 & \ldots  &-1  & -2  \cr
\ovl{2} & \phantom{-}0 & \ldots  &-1  & \phantom{-}1  \cr
\ovl{3} & \phantom{-}0 & \ldots  & \phantom{-}1  & \phantom{-}0  \cr
\ldots & \phantom{-}0  & \ldots  &  \phantom{-}0 & \phantom{-}0  \cr
\ovl{N\text{-1}} &-1 & \ldots  & \phantom{-}0  & \phantom{-}0 \cr
\ovl{N} & \phantom{-}1 & \ldots  & \phantom{-}0  & \phantom{-}0
\end{equation*}

As discussed earlier, the external-external CRN for the complete CRN will be taken to be empty. To specify the complete CRN, we now specify the boundary CRN $(\mcl{P}'_b)_g$ (suppressing $N$ because the CRN does not explicitly depend on the parameter). The first reaction is the production of the monomer $\ovl{1}$ by the boundary species $Y_0$, and all subsequent reactions are catalyzed production of the monomer in presence of different boundary species. In this model, different external boundary species have been coupled to the boundary reactions to simplify the CRN's decomposition into processes and transduction analysis (see App. <ref>).

\begin{align*}
(\mcl{P}'_b)_g &= \{ \Spc_g \cup \phantom{\cdot} \mcl{E}_g , (\R_b)_g, (\K_b)_g\}\\
(\Spc'_b)_g &= \{ \ovl{1}, \ovl{2}, \ovl{4}, \ldots, \ovl{2g}\},\\
(\mcl{E}_b)_g &= \{ Y_0, Y_2, Y_3,  \ldots, Y_{2g+1}\},\\
\{\mathcal{R}'_g,\mathcal{K}'_g\} &= \{\ce{Y_0 <=>[${c'}_0^+$][${c'}_0^-$] \ovl{1}},\\
& \phantom{=0} \ce{Y_2 + Y_0 + \ovl{2} <=>[${c'}_1^+$][${c'}_1^-$] \ovl{2} + \ovl{1} + Y_3},\\
& \phantom{=0} \ldots \\
& \phantom{=0}  \ce{Y_{2g} + Y_0 + \ovl{2g} <=>[${c'}_g^+$][${c'}_g^-$] \ovl{2g} + \ovl{1} + Y_{2g+1}}\}.
\end{align*}

The stoichiometric matrix for the boundary CRN is

\begin{equation*}
\bordermatrix{
~  &   & &   &     \cr
Y_{2g+1} & \phantom{-} 1 & \ldots  & \phantom{-}0  &  \phantom{-}0  \cr
Y_{2g} & -1 & \ldots  & \phantom{-}0  & \phantom{-}0  \cr
\ldots & \phantom{-}0  & \ldots  &  \phantom{-}0 & \phantom{-}0  \cr
Y_3 & \phantom{-} 0 & \ldots  &\phantom{-}1  & \phantom{-}0  \cr
Y_2 & \phantom{-}0 & \ldots  &-1  & \phantom{-}0  \cr
Y_0 & -1 & \ldots  & -1  & -1  \cr
\hline \cr
\ovl{1} & \phantom{-}1  & \ldots  &  \phantom{-}1 & \phantom{-}1  \cr
\ovl{2} & \phantom{-} 0 & \ldots  & \phantom{-}(1)  & \phantom{-}0 \cr
\ldots & \phantom{-}0  & \ldots  &  \phantom{-}0 & \phantom{-}0  \cr
\ovl{2g} & \phantom{-} (1) & \ldots  & \phantom{-}0  & \phantom{-}0
\end{equation*}

yielding the stoichiometric matrix of the external-boundary motif

\begin{equation*}
\nabla^{\mcl{E}}_{\Spc} = \quad \bordermatrix{
~  &   & &   &     \cr
Y_{2g+1} & \phantom{-} 1 & \ldots  & \phantom{-}0  &  \phantom{-}0  \cr
Y_{2g} & -1 & \ldots  & \phantom{-}0  & \phantom{-}0  \cr
\ldots & \phantom{-}0  & \ldots  &  \phantom{-}0 & \phantom{-}0  \cr
Y_3 & \phantom{-} 0 & \ldots  &\phantom{-}1  & \phantom{-}0  \cr
Y_2 & \phantom{-}0 & \ldots  &-1  & \phantom{-}0  \cr
Y_0 & -1 & \ldots  & -1  & -1  \cr
\end{equation*}

and the stoichiometric matrix of the internal-boundary CRN

\begin{equation*}
\St_b = \quad \bordermatrix{
~  &   & &   &     \cr
\ovl{1} & \phantom{-}1  & \ldots  &  \phantom{-}1 & \phantom{-}1  \cr
\ovl{2} & \phantom{-} 0 & \ldots  & \phantom{-}(1)  & \phantom{-}0 \cr
\ldots & \phantom{-}0  & \ldots  &  \phantom{-}0 & \phantom{-}0  \cr
\ovl{2g} & \phantom{-} (1) & \ldots  & \phantom{-}0  & \phantom{-}0
\end{equation*}

Since the reactions involve direct catalysis (same species in reactant and product sets) but the stoichiometric matrix misses that information, we denote the stoichiometric coefficient of the catalyzed species in parentheses.

Finally, we chemostat the boundary species $\{Y_0,Y_2, \ldots,Y_{2g+1}\}$ to obtain the reactions and rate constants (Eq. <ref>) of the internal CRN $\mcl{P}_{g,N}$. The reaction set of the internal CRN consists of the polymerization reactions, a single monomer production reaction, and catalyzed monomer production reactions. Assuming $N\geq 2g$, the internal CRN is given as:

\begin{align*}
\mcl{P}_{g,N} &= \{ \Spc_{N,g}  , \R_{N,g}\}\\
\Spc_{N,g} &= \{ \ovl{1}, \ovl{2}, \ldots, \ovl{N}\},\\
(\R_b)_g &= \{\ce{$\vn$ <=>[$c_0^+$][$c_0^-$] \ovl{1}},\\
& \phantom{=0} \ce{\ovl{2} <=>[$c_1^+$][$c_1^-$] \ovl{2} + \ovl{1} },\\
& \phantom{=0} \ldots \\
& \phantom{=0}  \ce{\ovl{2g} <=>[$c_g^+$][$c_0^-$] \ovl{2g} + \ovl{1}}\}\\
\R_{N,g} &= \R_N \cup (\R_b)_g.
\end{align*}

§.§ Attractor-transition graph

[Aspects of the internal CRN from Example <ref>]
Fixed-points, NEP, stochastic simulation data, and escape times of the internal CRN from e.g. <ref>.
(a) Fixed-points along with their stability are shown. (b) An estimate for the NEP along escape paths parameterized by the average polymer length is obtained using the AFGD algorithm [18]. (c) Runs for stochastic simulation for different volumes overlaid with the average polymer length of the fixed points are shown. Simulations are only shown through the first time the higher attractor is visited. (d) The escape times for these simulations is overlaid with a numerical estimate using the value of the NEP at the attractor with the lower average polymer length.

Following terminology from dynamical systems theory, we identify NESSs in an internal CRN with stable fixed points of the mass-action kinetics of the mean values of populations or attractors of the system. In this subsection, given a stochastic CRN, we explain how to determine its attractors and the expected time required by the system to escape each of them. From these data, we obtain the Attractor-Transition (A-T) graph. The A-T graph is a directed graph and is obtained as follows:

* Any two attractors adjacent to the same transient point are connected in the A-T graph.

* The direction points towards the attractor with a higher probability in the stationary distribution.

We use the A-T graph to formalize a concept of evolution in Sec. <ref>.

The relationship between a Markov process specified by a local transition-rate matrix to conditional probabilities of events separated in state space and time is well-understood through path-integral methods (see Ch. 10.4 in [2] and Sec. VI 3 in [19]). A short introduction and reference to mathematical and physics literature can be found in App. <ref>. A brief summary of the results needed for the construction of A-T graph is as follows.

In a stochastic CRN, events are specified by the vector of species population and time

\[ (n,t) \text{ for } n \in \mathbb{Z}^{\Spc}_{\geq 0}, t \in \mathbb{R},\]

along with a volume (also called scale factor) $V$ of the system which yields the species concentration

\[ q = \frac{n}{V} \in \mathbb{Q}^{\Spc}_{\geq 0}.\]

For a reversible CRN, the Hamiltonian is given as

\begin{equation}
H(p,q) = \sum_{r \in \mathcal{R}}
(e^{r^+\cdot p}-e^{r^-\cdot p})\left( k_{r}^+(e^{-p}q)^{r^-}
-k_{r}^-(e^{-p}q)^{r^+} \right), \label{eq:reversible_Hamiltonian}
\end{equation}

where $p \in \mathbb{R}^{\Spc}$ is the momentum vector canonically conjugate to the concentrations.
As explained in App. <ref>, the Hamiltonian function can be used to obtain leading order approximations to conditional probabilities of events, and through that, other quantities for understanding the stochastic dynamics of the system.

Consider the Hamiltonian from Eq. <ref> for an internal CRN. As explained in App. <ref>, the equations of MAK are obtained as (Eq. <ref>)

\begin{equation}
\dot{q} = \pdv{H}{p}\bigg|_{p=0}. \label{eq:MAK_Ham_sec}
\end{equation}

The set of all concentrations $\eq$ such that

\[ \pdv{H}{p}\bigg|_{(q,p) = (\eq,0)} = 0\]

are the fixed points of the MAK. If we determine the eigenvalues of the Hessian

\[ J = \pdv{H}{p}{q}\bigg|_{(q,p) = (\eq,0)}\]

at the fixed points, an attractor or a transient point are, respectively, fixed points with zero or a single eigenvalue with a strictly positive real part. As mentioned above, we identify a NESS as an attractor of the internal CRN and use the two interchangeably.

The solutions of the mass-action equations in Eq. <ref> are called decay or relaxation trajectories and connect states $q$ to fixed points $\eq$. A transient point is adjacent to an attractor if there is a mass-action solution that starts in the vicinity of the transient point and terminates at the attractor.
In stochastic CRNs, one can also find fluctuation or escape trajectories that take the system out of a NESS $\eq$ to a state $q$. As explained around Eq. <ref>, the momentum assignment along an escape trajectory is used to obtain the difference in non-equilibrium potential (NEP), denoted as $\V(q)$, which in turn can be used to obtain the stationary distribution $\pi(n)$ (when it exists,  [20, 21]) of the stochastic process using

\[ \pi(n=Vq) \asymp e^{-V \V(q)}.\]

In [18], we give a computational algorithm to estimate the NEP for multistable CRNs by finding the most probable escape or fluctuation paths.

The proportion of time the system spends in state $q$ converges to its stationary distribution value $\pi(q)$ as time runs to infinity (Ch. 1, [22]). The NEP can also be used to calculate the mean escape time of an attractor out of its domain of attraction.
Let the set $t^a = \{t^a_1,\ldots,t^a_N\}$ be the set of transient points adjacent to an attractor $a$. Let the NEP difference between $t^a_i$ and $a$ be given be $\Delta \V_i(a) \geq 0$. Define $\Delta \V^*(a)$ to be the minimum of $\Delta \V_i(a)$

\[\Delta \V^*(a) = \min_i \Delta \V(a).\]

Then, from Eq. 240 [19], the mean escape time from the domain of attraction of $a$ is given by

\begin{equation}
\tau_\text{esc}(a) = e^{V \Delta \V^*(a)}. \label{eq:escape}
\end{equation}

We can see the time spent by a system in an attractor gives a measure of the stability of the attractor. Thus, the depth of the NEP well around an attractor is a direct measure of the stability of the attractor.

We conclude this section by describing a trick for determining how the depth of the NEP wells around an attractor can be adjusted by varying the rate constants. Our argument relies on the continuity of the map from rate constants to fixed points of the MAK.
It is well-known (rederived in Th. <ref>) that the NEP is a Lyapunov function for solutions of MAK. Suppose a transient point $t$ has two adjacent attractors $a_1^t$ and $a_2^t$. Thus, the NEP at attractors is lesser or equal than the value at the adjacent transient point

\[ \V(t) \geq \V(a_1^t), \hspace{0.2em} \V(t) \geq \V(a_2^t).\]

This simply means that in the stationary distribution, the stochastic system is less likely to be at the transient point than the attractors.
However, notice that if the transient point $t$ was identified with an attractor, say $a_1$, i.e. $t = a_1$, then

\[ \V(t) = \V(a_1)> \V(a_2).\]

This means that $a_2$ would have a higher probability in the stationary distribution (lower NEP). By the same argument, if the transient point was very close to $a_2$, the probability in the stationary distribution of $a_1$ will be more than that of $a_2$. Using this observation, at least in simple cases, we can mathematically engineer the rate constants of multistable CRNs with a desired ordering on the stationary distribution. (Also, see e.g. <ref>).

Consider the class of internal CRNs denoted as $\mcl{P}_{g,N}$ from e.g. <ref>.
While the Hamiltonian for this model can be explicitly written, it is not particularly enlightening. Instead, we directly write the equations of MAK. Denoting the concentration of species $\overline{i}$ by $q_i$, and its rate of change by $\dot{q}_i$, the equations of MAK are

\begin{align*}
\dot{q}_n &= - \rcm{k}{n}q_n + \rcp{k}{n}q_1q_{n-1}, \\
\dot{q}_{n-1} &=  \rcm{k}{n}q_n - \rcp{k}{n}q_1q_{n-1}- \rcm{k}{n-1}q_{n-1} + \rcp{k}{n-1}q_1q_{n-2},\\
&\ldots \\
\dot{q}_{2} &=  \rcm{k}{3}q_3 - \rcp{k}{3}q_1q_{2} - \rcm{k}{2}q_{2} + \rcp{k}{2}q_1^2,\\
\dot{q}_1 &= ( \rcm{k}{2}q_{2} - \rcp{k}{2}q_1^2) + (\rcp{c}{0} - \rcm{c}{0}q_1)\\
& \phantom{=} + \sum_{i=2}^{N} \left(\rcm{k}{i}q_i - \rcp{k}{i}q_1 q_{i-1}\right) \\
& \phantom{=} +\sum_{j=1}^g\left( \rcp{c}{j} - \rcm{c}{j}q_1\right)q_{2j}.
\end{align*}

The fixed points of the system are obtained by solving for $\eq$ such that $\dot{q}(\eq) = 0$. Elementary considerations yield that the fixed points satisfy

\begin{align*}
\underline{q_i} &= \prod_{j=2}^{i}\left[\frac{\rcp{k}{j}}{\rcm{k}{j}}\right] (\underline{q_1})^i \text{ for } i \in [2,N],\\
0 &= (\rcp{c}{0} - \rcm{c}{0}\underline{q_1})+\sum_{j=1}^g\left( \rcp{c}{j} - \rcm{c}{j}\underline{q_1}\right)\underline{q_{2j}}.
\end{align*}

Substituting $\underline{q_{2j}}$ from the first equation into the second equation yields a polynomial whose roots are the equilibrium concentrations of the monomer. A simple application of Descartes' rule of signs yields that the maximum number of positive real roots is $2g+1$. Thus, we conclude that $\mcl{P}_{g,N}$ can have up to $g+1$ attractors.

As a concrete example, consider an internal CRN with $g=1$ and $N=4$ denoted as $\mcl{P}_{1,4}$

\begin{align*}
\mcl{P}_{1,4} &= \{ \Spc_{4,1}  , \R_{4,1}\}\\
\Spc_{4,1} &= \{ \ovl{1}, \ovl{2}, \ovl{3}, \ovl{4}\},\\
\R_{4,1} &= \{\ce{$\vn$ <=>[$c_0^+$][$c_0^-$] \ovl{1}},\\
& \phantom{=0} \ce{2\ovl{1} <=>[$k_2^+$][$k_2^-$] \ovl{2} <=>[$c_1^+$][$c_1^-$] \ovl{2} + \ovl{1} <=>[$k_3^+$][$k_3^-$] \ovl{3} },\\
& \phantom{=0}  \ce{\ovl{1} + \ovl{3} <=>[$k_4^+$][$k_4^-$] \ovl{4}} \},
\end{align*}

with rate constant assignment

\begin{align*}
k^+ &= 1 = \rcp{k}{i} &\text{ for } i &\in [2,4],\\
k^- &= 2.06 = \rcp{k}{i} &\text{ for } i &\in [2,4],\\
\rcp{c}{0} &= 0.1476, & \rcm{c}{0} &= 1,\\
\rcp{c}{1} &= 2.9, & \rcm{c}{1} &= 1.
\end{align*}

The vectors of fixed points $\eq = [\underline{q_1},\underline{q_2},\underline{q_3},\underline{q_4}]$ for this model are

\begin{align*}
\eq_1 &= [ 0.20, 0.019, 0.009, 0.00018]\\
\eq_2 &= [ 0.80, 0.31, 0.12, 0.046]\\
\eq_3 &= [ 1.9, 1.7, 1.6, 1.5].
\end{align*}

The first and the third fixed points are attractors and the second fixed point is a transient point, as shown in Fig. <ref>, panel (a).
Defining the average polymer length at a concentration as

\begin{equation}
\avg{\ell}_q = \sum_{i=1}^N i \frac{q_i}{\sum q}, \label{eq:avg_ell}
\end{equation}

the average polymer lengths at the three fixed points are:

\[ \langle \ell \rangle_1 = 1.1072, \hspace{1em} \langle \ell \rangle_2 = 1.54,\hspace{1em}  \langle \ell \rangle_3 = 2.399. \]

In the remainder of the example, we will call the first and third attractor the lower and higher attractor, respectively.

The escape paths joining the attractors to the transient point and the NEP along them are estimated using the action-functional gradient descent (AFGD) algorithm [18] and shown Fig. <ref>, panel (b), with average polymer length on the x-axis. As will be justified in the succeeding paragraph, we do not descend to the escape path completely but only run the algorithm for a few iterations to obtain an estimate. Thus, what is shown is an upper bound to the NEP. It can be seen that the attractor with the higher average polymer length (higher attractor) has a lower NEP value than that of the lower attractor. This means that, when running a stochastic simulation for long enough times, the process spends more time in the vicinity of the higher attractor than the lower one. The A-T graph for this model is drawn in panel $a$ of Fig. <ref>.

Stochastic simulations for different volumes or scale factors, namely $30$, $40$, and $50$, using the Gillespie algorithm [23] are displayed in Fig. <ref>, panel (c). An escape time is calculated for each run based on the first time step when the system crosses the concentration of transient point in all species. The escape times are then plotted in the bottom-right panel of the same figure along with a numerical estimate using the difference NEP values between the lower attractor and the transient point and substituting it in Eq. <ref>. One can see that estimates of the NEP are in good agreement with stochastic simulations. It is worth noting that the difference in value of the NEP at the transient point and higher attractor is a factor of $10$ higher, which means that the time of escape from the higher attractor is $e^{10}$ ($\approx 22,000$) more than that at the lower attractor. This is why stochastic simulation data for escapes from the higher attractor is omitted.

§.§ Complexity graph

In the previous subsection, we formalized NESSs as attractors of the internal CRN and outlined how to obtain an A-T graph. In this subsection, we define a complexity function $\mathfrak{C}$ for attractors

\[\mathfrak{C}: \text{NESS} \to \mathbb{R}_{\geq 0}\]

and define the complexity graph as follows. For every pair of adjacent attractors $a_1$ and $a_2$ in the A-T graph, make a directed edge $a_1 \to a_2$ if $\Cpx(a_1) < \Cpx(a_2)$. Thus, in the directed complexity graph, the arrow points towards the attractor with the higher complexity.

In computer science and mathematics, the subfield of algorithmic complexity assigns complexity measures to strings of data based on the complexity of algorithms that generate them  [24]. A popular measure is the logical depth introduced by C. Bennett in [25] and is defined as follows. Denote some universal computer as $U$, a program as $p$, the output of the computer on the program as $U(p)$, and the number of computational steps needed to perform $p$ as $T(p)$. The logical depth $D$ of a string $x$ is defined as the least number of computational steps needed to obtain it, or formally

\begin{equation}
D(x) = \min \{ T(p) | U(p) = x\}. \label{eq:logical_depth}
\end{equation}

In contrast, in mathematical and computational chemistry, it is convenient to model molecules as undirected graphs and chemical reactions as graph rewriting operations [26, 27, 28]. The species and reaction sets of a CRN can then be seen as generated through repeated application of graph rewriting operations on a starting set of molecules. Analogous to the complexity measure of logical depth for strings explained above and assembly in assembly theory [29, 30], we can define a measure of complexity for a molecular species as the minimum number of rewriting operations needed to produce that molecule. We can further generalize a complexity measure to a partial order where only molecules that are sequentially created downstream from a particular molecule have a higher complexity than that molecule, leaving the relative complexity of certain pairs of molecules generated through completely different mechanisms as incommensurable. Alternatively, one can also introduce a complexity measure for CRN configurations as the minimum number of reactions needed to reach [31, 32] it from a starting set of configurations. Notice that the above mentioned measures are purely topological and do not make explicit reference to the rate constant assignments.

While a detailed investigation of the notion of complexity for chemistry is outside the scope of this work, we nonetheless illustrate our formalism with a naive measure.
For models with relevance to the origins of biochemical life, we assume that we are given a complete CRN generated by an underlying set of graph rewriting rules along with a set of external species readily available from the environment. Following notation from Sec. <ref>, denote the sets of external and internal species as $\mcl{E}$ and $\Spc$, respectively.  We define the complexity $C(s)$ of an internal species $s$ as the minimum number of distinct reactions needed to generate it from $\mcl{E}$.

An algorithm to assign a complexity measure to species is as follows. Define a sequence of sets $S_0, S_1, \ldots, S_N$ such that:

Step $\mbf{0}$: Introduce all external species in set $S_0$

\[ S_0 = \{ s | s \in \mcl{E}\}\]

Step $\mbf{I}$: Introduce all species from the previous set $S_{I\text{-1}}$ in $S_I$ along with those species which can only be produced by reactions which only consumes species from $S_{I\text{-1}}$

\[ S_I = \{ s | s \in S_{I\text{-1}} \cup \{\text{supp}(r^+)|\text{supp}(r^-) \in S_{I\text{-1}} \hspace{0.2em} \forall r \in \R \}\}.\]

Terminate at $N$ when

\[S_N=S_{N+1}.\]

We define the complexity of species $C(s)$ of species $s$ as (compare with Eq. <ref>)

\begin{equation}
C(s) = \begin{cases}
\min \{ i | s \in S_i \},\\
\infty \text{ if } s \not\in S_N.
\end{cases}. \label{eq:complexity}
\end{equation}

Finally, we use the complexity of species to induce a measure of complexity $\Cpx(q)$ for a concentration $q \in \mathbb{R}^{\Spc}_{\geq 0}$. Since a NESS is specified by a concentration, we use this measure to assign a measure of complexity to attractors and obtain the complexity-graph described in the introduction to this subsection. Through any choice of non-linear function $h$ which depends on the concentration and complexity of species, one can abstractly define a complexity of concentrations as

\[ \Cpx(q) = \sum_s h(q,C(s)). \]

For this work, we investigate two of these choices namely:

Average complexity: sum of species complexity weighted by their relative concentrations

\[\Cpx^A(q) = \sum_s \frac{q_s}{\sum_t q_t} C(s). \]

Total complexity: sum of species complexity weighted by their absolute concentrations

\[\Cpx^T(q) = \sum_s q_s C(s). \]

Assembly: Assembly in Eq. 1 of [29, 30] is given as

\[A = \sum_s  \frac{q_s-1}{\sum_{t} q_{t}} e^{C(s)}.\]

Following e.g. <ref> and <ref>, in this example we assign measures of complexity to species and concentrations of the complete CRN $\mcl{P}'_{g,N}$.

The complexity of all external species $\mcl{E} = \mcl{E}_b$ is zero

\[C(\mcl{E})=0.\]

Due to the reaction $Y_0 \to \ovl{1}$, the complexity of the monomer $\ovl{1}$ is $1$

\[C(\ovl{1}) = 1.\]

Notice that the graph rewriting rule which corresponds to adding a new monomer to the end of the chain recovers the internal CRN.
In particular, due to reaction $2 \ovl{1} \to \ovl{2}$, the complexity of the dimer $\ovl{2}$ is 2

\[ C(\ovl{2}) = 2.\]

By a similar argument, it can be seen that the complexity of $\ovl{n}$ is $n$

\[ C(\ovl{n}) = n \hspace{1em} \forall n \in [1,N].\]

This follows because, in this CRN, only one monomer can be added at a time. Thus the species complexity simply measures the total number of monomers that constitute the polymer (also corresponds with the conserved quantity assigned to the cluster species).

Due to the observation above, for a concentration vector $q$ of polymers, the average complexity equals the average polymer size (see Eq. <ref>)

\[ \Cpx^A(q) = \avg{\ell}_q.  \]

Similarly, the total complexity equals the total number of monomers in the system

\[ \Cpx^T(q) = \sum_{i = 1}^N i q_i.\]

Using $\mcl{P}_{1,4}$ and indexing the lower and higher attractors in e.g. <ref> as $l$ and $h$, respectively, their complexity measures are

\begin{align*}
\Cpx^A(l) &= 1.1, & \Cpx^A(h) &= 2.4, \\
\Cpx^T(l) &= 0.27, & \Cpx^T(h) &= 16.
\end{align*}

Notice that by either measure, $\Cpx(l)< \Cpx(h)$. The complexity graph for this model is displayed in panel $b$ of Fig. <ref>.

§.§ Criteria for evolution

For obtaining the evolution graph for a CRN, $(1)$ start from the internal CRN, $(2)$ calculate non-equilibrium potential to extract the attractor-transition (A-T) graph, $(3)$ assign complexity to species and extract the complexity graph, $(4)$ overlay the two to get the evolution graph. In the shown graphs, the vertices are attractors (NESSs) of the system, dashed arrows represent the complexity relationships, and the solid arrows represent the relative probabilities of the pair of adjacent attractors under the stationary distribution.

Visualizing different types of evolution. Panel $a$ and $b$ display progressive evolution and historically contingent evolution, respectively. Panels $c$ and $d$ are both examples of open-ended evolution, and, in particular, panel $d$ shows historically contingent open-ended evolution.

Equilibrium concentrations of monomers for a 2-monomer copolymerization model CCRN from e.g. <ref> is shown. As explained there, there exists a rate constant assignment such that the attractors further from the origin can be made both more complex and more probable in the stationary distribution than the attractor closer to the origin. Furthermore, since there is no transient point between the two, the found system would be capable of exhibiting historical contingent evolution from the attractor closer to the origin.

Earlier (Sec. <ref>) we saw how to use a complete CRN and an external species set to obtain an internal CRN. Then (Sec. <ref>) we explained how to obtain its NESSs or attractors, estimate their probability in the stationary distribution, and encode them in a directed A-T graph. Furthermore, in Sec. <ref>, we explained how to assign a complexity to each attractor and obtain a complexity graph. In this subsection we explain how, using the A-T and complexity graphs, we can determine whether a CRN is capable of evolution, and, if so, of which type.

Recall that an A-T graph is a directed graph with attractors as vertices and directed edges pointing towards an adjacent attractor with higher probability in the stationary distribution. We use solid lines to denote the edges in the A-T graph. Furthermore, a complexity graph is also a directed graph with attractors as vertices and directed edges pointing towards a more complex adjacent attractor. We use dashed lines to denote edges in the complexity graph. We overlay the A-T graph and complexity graph to obtain the evolution graph and use it to identify evolvability properties of the CRN. A graphical outline of this procedure is shown in Fig. <ref>.

In the remainder of this subsection, we specify our criteria for progressive, historically contingent, and open-ended evolution.

Progressive evolution:

A CRN will be said to allow progressive evolution from an attractor if there exists an adjacent attractor that is:

* More complex than the current attractor.

* More probable than the current attractor in the stationary distribution.

Graphically, as shown in panel $a$ of Fig. <ref>, this corresponds to a pair of attractors on the graph connected by a solid and a dashed arrow in the same orientation. In that case, the CRN will be said to be capable of progressive evolution from the attractor at the base of the arrows.

The evolution graph of the model $\mcl{P}_{1,4}$ from e.g. <ref> is shown in the $c$ panel of Fig. <ref>. Since the higher attractor $a_h$ is more complex and more probable in the stationary distribution than the lower attractor $a_l$, the system is capable of progressively evolving from $a_l$ to $a_h$. Observe that the evolution graph for this model is equivalent to panel $a$ of Fig. <ref>.

Historically contingent evolution:
A CRN will be said to be capable of historically contingent transition from an attractor if there exist at least two adjacent attractors, such that they are:

* More probable than the current attractor in the stationary distribution.

* Not adjacent to each other.

Moreover, if the two attractors are also more complex than the current attractor, then the CRN will be said to be capable of historically contingent evolution from the less complex attractor.

Graphically, as shown in panel $b$ of Fig. <ref>, historically contingent evolution corresponds to having at least three attractors such that one of them is connected to the other two with a solid and a dashed edge with the tail at the first attractor, and the two are not connected with a solid edge.

We justify the association of the above criterion with historical contingency [33, 34] because (as explained in App. <ref>) the most likely transition between attractors goes through a transient point between them. Since the more complex attractors are not directly linked by a transient point between them, it is highly unlikely that the system would transition directly from one of these attractors to the other without going through the less complex attractor. However, since the less complex attractor is also less probable in the stationary distribution, the transition from the more complex to a less complex attractor is also relatively unlikely. Thus, the stochastic choice made by the system in transitioning out of the attractor with lower complexity, seals the fate of the system in a probabilistic sense. This phenomenon is known as symmetry breaking or ergodicity breaking in physics and mathematics literature, respectively [35, 36].

Consider a 2-monomer copolymerization model [12] consisting of two different types of monomers, labelled $A$ and $B$. Denoting a polymer as $\omega$, the reactions where only end additions and removals are allowed with some (polymer independent) rates $s$ and $d$, respectively, are given as:

\begin{align*}
\ce{$\omega$ + A <=>[s][d] $\omega A$}, \hspace{1em}   \ce{$\omega$ + B <=>[s][d] $\omega B$}.
\end{align*}

Notice that the state space scales exponentially in the size of the polymer as there are $2^n$ polymers of size $n$ in two monomers. To simplify analysis, we map strings to clusters (for example, the polymer $AABA$ maps to the cluster $\ovl{3,1}$) and analyze the cluster CRN (CCRN) [14] associated with the copolymerization model, as described below.

The species set of a two constituent CCRN is given by

\[\mathcal{S} = \{ \overline{1,0},\overline{0,1},\overline{1,1},\ldots,\overline{i,j},\ldots\},\]

for $i,j \geq 0 \in \mathbb{Z}$. The internal-internal reaction set induced by end additions and removals is given by

\begin{align*}
\mathcal{R}_i =&  \bigg\{ \ce{\overline{i,j} + \overline{1,0} <=>[s][d] \overline{i+1,j}}\hspace{0.1em}, \hspace{0.2em} \ce{\overline{i,j} + \overline{0,1} <=>[s][d] \overline{i,j+1}}  \bigg|\\
& \phantom{=} \forall i,j \in \mathbb{Z}_{\geq 0},  N \geq i+j\geq 1\bigg\},
\end{align*}

where $s$, $d$, and $N$ are rates of addition, removal, and maximum polymer length, respectively.

Analogous to the internal CRN of the 1-monomer polymerization models from e.g. <ref>, we introduce reactions for monomer production and catalyzed monomer production. We assume that the clusters $\ovl{2,0}$ and $\ovl{3,2}$ catalyze the production of $\ovl{1,0}$ and $\ovl{0,2}$ and $\ovl{2,3}$ catalyze the production of $\ovl{0,1}$. A concrete realization of the model with rate constants is given by the reaction set:

\begin{align*}
\mathcal{R} = & \mathcal{R}_i|_{s=1, d= 5.01} \cup \\
& \bigg\{ \hspace{0.5em} \ce{ $\overline{0,1}$ <=>[17][10] $\varnothing$ <=>[10][17] $\overline{1,0}$} \hspace{0.3em}, \\
& \phantom{=} \ce{$\overline{2,0}$ <=>[40.08][5.01] $\overline{2,0} + \overline{1,0}$ } \hspace{0.3em}, \hspace{0.5em}
\ce{$\overline{3,2}$ <=>[31.5][13.54] $\overline{3,2} + \overline{1,0}$ },\\
& \phantom{=} \ce{$\overline{0,2}$ <=>[40.08][5.01] $\overline{0,2} + \overline{0,1}$ } \hspace{0.3em}, \hspace{0.5em}
\ce{$\overline{2,3}$ <=>[56.7][36.82] $\overline{2,3} + \overline{0,1}$ }\bigg\}.
\end{align*}

The attractors in the monomer concentrations are shown in Fig. <ref>. It can be seen that the CRN has three attractors, roughly around $(1,1),(4,1)$ and $(1,3)$, and we denote them as $a_1$, $a_2$, and $a_3$, respectively.

Notice that while there is a transient point between $a_1$ and $a_2$, and $a_1$ and $a_3$, there is no transient point between $a_2$ and $a_3$. This means that $a_1$ is adjacent to both $a_2$ and $a_3$, however $a_2$ and $a_3$ are not adjacent to each other. By making the transient points much closer to $a_1$ than $a_2$ or $a_3$, as explained in Sec. <ref> and illustrated in e.g. <ref> and <ref>, the NEP at $a_1$ can be made higher than both $a_2$ or $a_3$. Thus, there exists an assignment of rate constants such that $a_2$ and $a_3$ are more probable than $a_1$ in the stationary distribution. Furthermore (by the same arguments as in e.g. <ref>) the polymer concentration in the equilibrium distribution follows an exponential distribution and the average polymer length as well as the total monomer concentration in $a_2$ and $a_3$ is higher than in $a_1$. Thus, using our complexity measure for concentrations, $a_1$ is the least complex attractor.
This means that there exists a rate constant assignment such that $a_2$ and $a_3$ are more complex and probable in the stationary distribution than $a_1$, and also not adjacent to each other. The evolution graph for such a system is equivalent to panel $b$ in Fig. <ref> and the system can exhibit historically-contingent evolution from $a_1$.

Open ended evolution:
Stochastic CRNs can exhibit a wide variety of dynamics such as limit cycles, chaotic attractors [37], and transience [38]. Transience in CRNs has been used to explore polymer models that grow indefinitely in [12, 13, 39]. However, since this work is limited to point attractors, a generalization to other attractors and transience will be left to future work.

A CRN will be said to be capable of open-ended evolution if there is at least one subset of infinitely many attractors such that each one is capable of progressive evolution to another attractor in the subset.

Some graphs displaying open-ended evolution are shown in panels $c$ and $d$ of Fig. <ref>. Notice that there is only one open-ended direction for the system shown in panel $c$, however panel $d$ has two historically-contingent directions for exhibiting open-endedness.

As explained earlier (in e.g. <ref>), $\mcl{P}_{g,N}$, with $N \geq 2g$, can have up to $2g + 1$ fixed points and $g+1$ attractors. Furthermore, by the arguments in Sec. <ref>, a rate constant assignment can always be made such that the NEP of the more complex adjacent attractor is also lower, thus ensuring that the more complex attractor is more probable in the stationary distribution. A one dimensional analog of the model with $3$ attractors is shown in e.g. <ref>.

Note that the extension of this model to infinite positive attractors is technically challenging as the rate constants also increase, raising issues of summability and convergence. Thus, while we propose a criterion for open-ended evolution, a proper treatment of reaction networks with infinite reactions and fixed-points will be left as future work.

Applying the criteria explained in this subsection to the evolution graph displayed in panel $(4)$, Fig. <ref>, it can be seen that the system is:

* Capable of progressive evolution from attractor $a_1$, $a_2$, $a_4$, $a_5$.

* Capable of historically contingent evolution from attractor $a_2$.

* Incapable of open-ended evolution because of a finite number of attractors.

§ REMARKS ON FORMALISM AND MODEL

In this section, we remark on the roles of entropy production rate (EPR) and autocatalysis in our formalism and models. In Sec. <ref>, we explain the indirect role of EPR in our definition of evolution. In Sec. <ref>, we defined a class of polymerization models and demonstrated that a model with $g$ polymer catalysts can exhibit up to $g+1$ attractors. For these models, in Sec. <ref>, we explain the role of autocatalysis in making multistability possible.

§.§ Role of dissipation and efficiency in formalism

In [40], using fluctuation theorems concerning dissipation and irreversibility proposed in [41, 42], England suggests a relationship between biological organization and EPR (for fluctuation theorems in CRNs, see [43]). In particular, it is conjectured that biological systems evolve to maximize entropy production or dissipation. Similarly, in [44], complex systems are found where the system temporally evolves to states of higher EPR. In this subsection, we briefly summarize the role played by EPR and free-energy transduction efficiency in our formalism. (See App. <ref> for a review of related concepts.)

First, one can show that the EPR only informs about the stationary probability distribution of the system when the system admits a single detailed-balanced or complex-balanced equilibrium in each stoichiometric compatibility class (see discussion after Eq. <ref>). In particular, the expression for EPR and the integrand of the NEP with time as the variable of integration coincide (up to sign) only for detailed-balanced systems. Since we are primarily interested in multistable systems with many NESSs, the EPR cannot be used to determine relative probabilities of the NESSs or transitions between them. Thus, EPR cannot be used to obtain an attractor-transition graph (see Sec. <ref>).

Second, in Sec. <ref>, we explained how internal reactions can be obtained through the application of graph-rewriting operations, and used them to assign a complexity measure to species that induces a complexity measure on the NESSs. However, because EPR at a NESS only depends on the chemical potential and fluxes of the external boundary species (see Eq. <ref>), it does not depend on the internal-internal species that make up the equilibrium concentration profile of the NESS. Hence, EPR cannot be used to obtain a complexity measure on the NESSs either.

In summary, the EPR at NESSs need not always be correlated with either their probability in the stationary distribution or their complexity. Nonetheless, as reviewed in App. <ref>, one can decompose a CRN into processes and obtain a free-energy transduction efficiency using the EPRs of reactions at each of the NESSs. Thus, in principle, the efficiency (or EPR) could be used as a measure for ranking the NESSs and obtaining a dissipation analog of a complexity graph. The relevance of such a ranking to the origin of life problem is, however, unclear.

As shown in e.g. <ref> and <ref>, the steady state behavior of the 1-monomer polymerization models from e.g. <ref> is identical to the complete generalized Schlögl model introduced in e.g. <ref>. Since the efficiency at NESSs only depends on the steady state behavior, the efficiency analysis of the polymer model is identical to that of the Schlögl model.

An efficiency analysis of NESSs of the generalized Schlögl model, shown in e.g. <ref>, demonstrates that (Fig. <ref>) the efficiency increases monotonically in the distance of the attractor from the origin and does not correlate with the NEP or the probability of the attractor in the stationary distribution. Although this example is a case where the efficiency at the various attractors is positively correlated with their complexity, unless a theorem is proved otherwise, this cannot be assumed to always be the case.

§.§ Role of autocatalysis in models

In Sec. <ref>, we defined a class of polymerization models and demonstrated that a model with $g$ polymer catalysts can exhibit up to $g$ attractors.
In this section, we explain that these polymer catalysts are autocatalysts, and a model with $g$ polymer catalysts has exactly $g$ minimal autocatalytic subnetworks. We term the species like the catalysts as functional species, and elaborate on conditions satisfied by the models that yield multiple attractors.

§.§.§ Review: Autocatalysis

Autocatalysis is an important property of CRNs that generalizes the biologically relevant notions of self-replicability and ecology [45, 46, 47, 14]. In [47], it is shown that different autocatalytic subnetworks can be composed to achieve a variety of ecologically relevant dynamical behavior including competition, predator-prey, and mutualism. In [48, 49], it is shown that autocatalytic subnetworks on a spatial lattice exhibit Turing patterns and other spatial patterns relevant to the origins of life. The kinetic role of autocatalysis is explained in [50], and some thermodynamic results about autocatalysis are shown in [51].

Let $\mcl{G} = (\Spc,\R)$ denote a subnetwork of $\mcl{G}'= (\Spc',\R')$, where $\R \subset \R'$ and $\Spc$ consists of all the species of $\Spc'$ that participate in $\R$. For $\mcl{G} = (\Spc,\R)$, we denote the input, output, and stoichiometric matrix as $\St^-$, $\St^+$, and $\St$, respectively (see App. <ref>), and the restriction of a $\Spc \times \R$ matrix $\mathbb{M}$ to a set $\mcl{M} \subset \Spc$ as $\mathbb{M}|_{\mcl{M}}$. The reactions in the complement of $\mcl{G}$ are given by all reactions in $\R'$ that are not in $\R$ and denoted as $\R'/\R$.

A subnetwork is autocatalytic in the autocatalytic set if there exists a reaction flow such that every species in the autocatalytic set is consumed but also net-produced.
Formally, a subnetwork $\mcl{G} = (\Spc,\R)$ is autocatalytic in the species set $\mcl{M}\subseteq \Spc$ if (see [14] and references therein):

* Autonomy: Every species in $\mcl{M}$ is consumed and produced by at least one reaction in $\R$, and every reaction consumes and produces at least one species in $\mcl{M}$,

\begin{align*}
\forall s \in \mcl{M}, \exists r,r' \in \R : &\quad (\St^{-})^r_s > 0 \text{ and } (\St^{+})^{r'}_s > 0,\\
\forall r \in \R, \exists s,s' \in \mcl{M} : &\quad (\St^{-})^r_s > 0 \text{ and } (\St^{+})^r_{s'} > 0.
\end{align*}

* Productivity: There exists a semi-positive reaction flow vector such that all species in $\mcl{M}$ are net produced, (see Sec. <ref> for notation)

\[ \exists \mbf{x} \in \mathbb{R}^{\R}_{\geq 0} : \quad \St|_{\mcl{M}} \mbf{x} \gg \mbf{0}.\]

Under these conditions, $\mcl{M}$ is an autocatalytic set. The restriction of the stoichiometric matrix to the autocatalytic set  $\mathbb{S}\big|_\mathcal{M}$ is a semi-positive (SP) matrix (Ch. 3, [52]).
If each $\mbf{x}$ that maps in the productive region is strictly positive, or $\kappa \gg \mathbf{0}$, then $\mathbb{S}\big|_\mathcal{M}$ is minimally SP (MSP) and the reaction network $\mathcal{G}$ is minimally productive in $\mathcal{M}$. By definition, no reaction-reduced subnetwork of a minimal productive network can be productive in the autocatalytic set. One can further restrict the set of autonomous species to find a minimal autocatalytic motif, also called an autocatalytic core. It can be shown that the stoichiometric matrix of an (autocatalytic) core is an invertible matrix [46]. We refer to the species in the core as core species.

§.§.§ Heuristics for multistability

There are two aspects of multistability: kinetic and topological.

With regards to the kinetic requirement for multistability, recall from Sec. <ref> that for thermodynamic consistency, the rate constants on the complete network are assumed to be such that the system has a single detailed-balanced equilibrium in each stoichiometric compatibility class. However, as explained below Eq. <ref>, by introducing an environmental control in the form of chemostatted external-boundary species, the rate constants of the system can depart from a detailed-balanced equilibrium. Thus, the first heuristic for multistability is that there must be at least one boundary reaction in the complete CRN.

In regards to the topological requirement for multistability, non-zero deficiency (Eq. <ref>) is a necessary condition for multistability in a reversible CRN (see App. <ref>). The appendix also explains that the deficiency, $\delta$, of a CRN counts the number of null flows, termed $\mbf{\delta}$-flows, that are not supported in a single linkage class. From Sec. <ref>, recall that each autocatalytic subnetwork has a core whose stoichiometric matrix is invertible. This means that there exists a reaction-flow vector through the autocatalytic subnetwork that creates exactly one core species. If there exists a reaction-flow in the complement of the core that also produces any of the core species, and if the above mentioned two reaction-flows have a support in at least two linkage classes, then the deficiency of the network must be non-zero. Thus, the presence of autocatalytic cores is our second heuristic for multistability.

To summarize, the existence of at least one boundary reaction and at least one autocatalytic core are our two heuristics for multistability. We term an internal-boundary species that is also a core species a functional species. Such a species can perform the dual role of transducing free-energy from the environment and super-linear growth in some concentration regimes, satisfying a few conditions necessary for multistability.

Consider a subnetwork $\mcl{H} = (\Spc_\mcl{H},\R_\mcl{H})$ of the internal CRN $\mcl{P}_{1,2}$ from e.g. <ref>

\begin{align*}
\Spc_\mcl{H} &= \{ \ovl{1},\ovl{2}\},\\
\R_\mcl{H} &= \{ 2\ovl{1} \to \ovl{2}, \ovl{2} \to \ovl{1}+\ovl{2}\}.
\end{align*}

Observe that $\mcl{H}$ is autonomous in $\Spc_\mcl{H}$, as both species are produced and consumed and each reaction produces and consumes both species. Moreover, the network is also productive as all the species are net produced under the flow $\mbf{x} = [1,3]^T$,

$$\St \mbf{x} = \begin{bmatrix}
\end{bmatrix} \gg \mbf{0}.$$
Thus, $\mcl{H}$ is an autocatalytic subnetwork. Furthermore, due to the reaction $\ovl{1} \to \vn$ in the CRN $\mcl{P}_{1,2}$, the autocatalytic subnetwork produces a $\delta$-flow

$$ \ovl{2} \to \ovl{1} + \ovl{2}, \quad \ovl{1} \to \vn.$$

As explained in e.g. <ref>, due to the external species in the network, $\ovl{1}$ and $\ovl{2}$ are both internal-boundary species. Since they are both internal-boundary species and autocatalytic species, we call them functional species. We have shown that there exist rate constants such that $\mcl{P}_{1,2}$ exhibits multistability, and we can attribute it to the presence of functional species. By the same argument, one can see that for a model $\mcl{P}_{g,N}$, there are $g$ independent autocatalytic subnetworks and $g+1$ functional species. Moreover, each autocatalytic subnetwork creates an independent $\delta$-flow, and there exist rate constants such that the system exhibits $g+1$ attractors.

§ DISCUSSION

In this work we have formalized phases in CRNs under mass-action kinetics (MAK) as attractors (stable fixed points). Although closed and thermodynamically consistent systems can only exhibit a single detailed-balanced attractor, open systems can exhibit multistability (i.e. multiple attractors or non-equilibrium steady states). We propose that multistability is a precondition for evolution, and define four steps to ascertain whether a multistable CRN is capable of progressive evolution. First, start with a complete CRN and obtain its internal-boundary-external partition by specifying the species that are internal to the system under investigation and its environment. We restrict this work to scenarios where the environmental species are chemostatted and refer to the resulting network as the internal CRN (Sec. <ref>). Second, determine the attractors for the internal CRN and their adjacencies using MAK, and assign to each attractor a stability using its non-equilibrium potential (NEP) (App. <ref>). This information is graphically represented in an attractor-transition graph  (A-T graph) with vertices as attractors, and directed edges between adjacent attractors pointing towards the more stable attractor (Sec. <ref>). Third, superimpose a complexity graph on the A-T graph with the same vertices and edges, but where the edge direction is from less complex to more complex adjacent attractors. This step requires a (partial) ordering of the attractors by their complexity, which we define using measures from algorithmic complexity (Sec. <ref>). Finally, we superimpose the two graphs to obtain an evolution graph, and specify that a CRN is capable of progressive evolution from an attractor if there exists another attractor adjacent to it that is more stable and more complex (Sec. <ref>).

In Sec. <ref>, we explain that evolution does not always correspond with an increase in entropy production rate or increase in free-energy transduction efficiency of the attractors. A further technical contribution of this work lies in defining the concept of a functional species and exploring its connection to multistability. A functional species has two key roles namely: (1) interacting with the environment to transduce chemical potential energy into the system needed to maintain it at a non-detailed balanced attractor, and (2) autocatalysis in the internal species. As remarked in Sec. <ref>, within a range of rate constants for our models, these interactions introduce additional attractors at higher average polymer lengths, while also enforcing that the attractor with a higher average length is also more stable. By a similar construction, in Sec. <ref>, we show that adding a functional polymer species in a two monomer polymerization CRN can make the system capable of historical contingency, and that adding a countable number of functional species can, in principle, make the system open-ended. Thus, we demonstrate that, at least mathematically, one can find CRN models that exhibit behaviors that resemble aspects of biochemical life.

Our work has several direct implications for the study of the origins of biochemical life. Given a prebiotically plausible CRN with rate constants, the proposed framework can be used to numerically ascertain whether the CRN is capable of progressive, historically-contingent, or open-ended evolution under given environmental conditions. The distribution of attractors and their relative stabilities can be found using numerical solvers for finding attractors (such as Bertini [53]) and algorithms for estimating the NEP between them [18]. The concept of an evolutionary graph could readily be employed beyond mass-action to other kinetic models such as Michaelis-Menten [54], used to account for modes of environmental control from the environment other than chemostatting [55], for protocols that deviate from continuously-stirred tank reactors (such as serial-dilution experiments [56]), or to incorporate other types of dynamics such as limit cycles, chaotic attractors [57], and transience [38].

By combining advanced techniques in probability theory with foundational concepts in evolutionary biology, we have rigorously defined the concept of evolution for stochastic CRNs. Similar to the Gibb's free energy in statistical mechanics [58], our framework relies on the NEP for stochastic processes [59] for identifying attractors and their relative stabilities. A question for advancing rigorous phenomenological investigations into biology that our approach opens up is whether all biological phenomena have a counterpart in stochastic CRNs, or some generalization thereof? If there is biological phenomenon without a stochastic counterpart, then it sheds light on the limitations of any stochastic framework for the origins of life, such as the one proposed here. If, however, all biologically relevant phenomena can be formalized in a stochastic framework, then it would suggest that there is a natural path of complexification from prebiotic chemistry to modern life [60], with foundations in non-equilibrium statistical physics.

§ ACKNOWLEDGEMENTS

This work was funded by the National Science Foundation, Division of Environmental Biology (Grant No: DEB-2218817). I thank Christoph Flamm, David Anderson, Diego Rojas La Luz, Massimo Bilancioni, Massimiliano Esposito, and Nicholas Guttenberg for technical conversations and references, Vladimir Sotirov, Eric Smith, and David Baum for help in conceptualizing and clarifying several aspects of this work, and David Baum for a thorough proofreading and editing of the manuscript.

§ MATHEMATICAL PRELIMINARIES

§.§ Basic terminology

* The set of real numbers and its subset of integers are denoted by
$\mathbb R$ and $\mathbb Z$.

* For any subset $X\subseteq\mathbb R$ of real numbers, its subsets of
non-negative and positive elements are denoted by
$X_{\ge 0}=\{x\in X:x\ge 0\}$ and $X_{>0}=\{x\in X:x>0\}$.

* $\mbf{0}$ denotes a vector of zeros.

* Given any sets $X$ and $Y$, the set of $Y$-valued functions with domain
$X$ is denoted by $Y^X$.

* A real-valued function with domain a set $X$, i.e. an element
$v\in\mathbb R^X$, is called an $X$-vector.

* The value at $x\in X$ of $\mbf{v}\in\mathbb R^X$ is its
$x$-component and is denoted by $\mbf{v}_x$.

* Given two $X$ vectors $a$ and $b$, we use $a^b$ to denote $\prod_{i\in X} a_i^{b_i}$.

* A $\mathcal S\times\mathcal R$ matrix is a two-variable
function with domain a pair of finite sets $(\mathcal S ,\mathcal R)$,
i.e. an element $\mathbb{M}\in\mathbb R^{\mathcal S\times\mathcal R}$.

* The value at $(s,r)\in\mathcal S\times\mathcal R$ of
$\mathbb{M}\in\mathbb R^{\mathcal S\times \mathcal R}$ is the
$(s,r)$ entry $\mathbb{M}^r_s\in\mathbb R$ of the matrix $\mathbb{M}$.

* $\mathbb{M}^r$ is the $\mathcal S$ vector satisfying $(\mathbb{M}^r)_s=\mathbb{M}^r_s$ for all $r\in\mathcal R$.

* $\mathbb{M}_s$ is the $\mathcal R$ vector satisfying $(\mathbb{M}_s)_r=\mathbb{M}^r_s$ for all	$s\in\mathcal S$.

§.§ Chemical reaction networks

A chemical reaction network (CRN) $\mathcal{G}=\{\Spc,\C,\R\}$ is a hypergraph[A graph is defined as a set of vertices and a set of edges linking the vertices. A hypergraph is a generalization of a graph with a set of hyperedges which can link multisets of vertices, also called hypervertices. (A multiset is a generalization of a set that allows for multiple instances of each element.)]
consisting of three sets (see [61]):

* a set $\mathcal{S}$, elements of which are the species of the network (vertices of the hypergraph).

* a set $\mathcal{C}$ of $\mathcal{S}$ (column) vectors in $\mathbb{Z}_{\geq 0}^{\mathcal{S}}$, elements of which of multisets of species and called the complexes of the network (hypervertices of the hypergraph).

* a set $\mathcal{R} \subset \mathcal{C} \times \mathcal{C}$ such that (hyperedges of the hypergraph)

* for each $y \in \mathcal{C}$, $(y,y) \not\in \mathcal{R}$.

* for each $y \in \mathcal{C}$ there is a $y' \in \mathcal{C}$ such that $(y,y') \in \mathcal{R}$ or $(y',y) \in \mathcal{R}$.

Members of $\mathcal{R}$ are the reactions of the network. A reaction $r \in \mathcal{R}$ is of the form $(r^-,r^+)$, where $r^-$ and $r^+$ are the input and output complexes, respectively, in $\mathcal{C}$. In what follows, we use the notation $r^- \to r^+$ to denote $(r^-,r^+)$.

The input (output) matrix $\St^-$ ($\Sp$) is an $\Spc \times \R$ matrix whose columns are the input (output) complexes, $(\St^-)^r = r^-$ ($(\Sp)^r = r^+)$. The difference of the output and input matrices is called the stoichiometric matrix and denoted as $\St$, where $\St = \Sp-\St^-$. The columns of $\St$ are called the reaction vectors and the $r^\text{th}$ column is denoted as $\Delta r$, where $\Delta r = \St^r = (r^+-r^-).$

Define the $\mathcal{S} \times \mathcal{C}$ complex stoichiometry matrix $Y$, such that $Y^y = y \in \mathcal{C}.$ Since the CRN is a hypergraph, it is a graph over the complexes (hypervertices). Define the $\mathcal{C}\times \mathcal{R}$ vertex-edge matrix $\mathbb{M}$ such that

\[ \mathbb{M}_y^r = \begin{cases}
+1 \text{ if } y = r^+,\\
-1 \text{ if } y = r^-\\
\phantom{+}0 \text{ otherwise.}
\end{cases} \]

Then the stoichiometric matrix is given as the product of the matrix of complex stoichiometries and the vertex-edge matrix
$\mathbb{S} = Y \mathbb{M} .$

A concentration vector is an $\mathcal{S}$ vector $q \in \mathbb{R}^\mathcal{S}_{\geq 0}$ that consists of the concentrations of the species.
An $\mathcal{R}$ (column) vector in the positive-orthant of the domain of the stoichiometric matrix $j \in \mathbb{R}_{\geq 0}^\mathcal{R}$ will be called a reaction-flow vector.
The image of any reaction-flow vector under the stoichiometric matrix will be called a species-flow vector. Under a reaction-flow vector $j$, the rate of change of the concentration $q$ will be given by the species-flow vector $\dot{q}$ such that

\begin{align}
\dot{q} &= \mathbb{S}j. \label{eq:rate_equation}
\end{align}

Note that $\dot{q}$ and $j$ have the same units.

Any $\mathcal{S}$ (column) vector $c^* \in \text{Ker}(\mathbb{S}^T)$ is called a conservation law [61]. We denote the number of independent conservation laws or dimension of $\text{Ker}(\mathbb{S}^T)$ by $\ell$.

The image of the stoichiometric matrix $\text{Im}(\St)$ is called the stoichiometric subspace. By the rank-nullity theorerm, we know that the dimension of the stoichiometric subspace $\text{dim}(\text{Im}(\St)) = |\Spc|-\ell$. Two concentration vectors $q$ and $q'$ are stoichiometrically compatible if $q - q' \in \text{Im}(\St)$. Stoichiometric compatibility is an equivalence relation that induces a partition of $\mathbb{R}^{\Spc}_{\geq 0}$ into equivalence classes called the stoichiometric compatibility classes (Def. 3.4.6, [61]).

Any reaction-flow vector $j \in \text{Ker}(\mathbb{S})$ is a null flow. Null flows can arise because either $j \in \text{Ker}(\mathbb{M})$ or $\mathbb{M}j \in \text{Ker}(Y)$. Any null flow in $\text{Ker}(\mathbb{M})$ and $\text{Ker}(Y)$ will be called a null cycle and $\mbf{\delta-}$flow [62], respectively. We denote the dimension $\text{Ker}(\mathbb{M})$ and $\text{Ker}(Y)$ as $\iota$ and $\delta$, respectively. Elementary considerations yield that for a CRN $\{\mathcal{S},\mathcal{C},\mathcal{R}\}$ (App. A, [14]),

\begin{equation}
\delta = |\mathcal{R}| - |\mathcal{S}| - \iota + \ell. \label{eq:deficiency}
\end{equation}

$\delta$ is called the deficiency of the network.

Mass-action kinetics (MAK) is a kinetic model (Def. 1, [50]) for a CRN $\{\mathcal{S},\mathcal{C},\mathcal{R}\}$ where:

* each reaction $r \in \mathcal{R}$ is assigned a rate constant $k_r \in \mathbb{R}_{>0}$.

* for a given concentration vector $q$, a reaction-flow vector $j$ is assigned such that for each reaction $r: r^- \to r^+ \in \R$,
$j_r = k_r q^{r^-}$.

Thus, under MAK, the dynamics of the concentration vectors $q$ is given by the ODES

\begin{equation}
\dot{q}(q) = \sum_{r \in \R}\mathbb{S}^r j_r = \sum_{r \in \R}\mathbb{S}^r k_r q^{r^-}. \label{eq:MAK_eqs}
\end{equation}

We denote the set of all rate constants as $\mathcal{K}$, and henceforth specify a CRN under MAK with the quadruple $\{\mathcal{S},\mathcal{C},\mathcal{R},\mathcal{K}\}$. A concentration $\eq$ is an equilibrium concentration of the system if
$$\dot{q}(\eq) = 0.$$

Reversible CRNs:

A CRN is called reversible if for every forward reaction $r:r^- \to r^+ \in \mathcal{R}$, there exists its reverse reaction $r': r^+ \to r^- \in \mathcal{R}$.
For reversible CRNs, we employ the following notation. First, we obtain a unique set of reactions $\R_u \subset \R$ which contains either the forward or reverse reaction for each reaction in $\R$, but not both. We denote the $\Spc \times \R_u$ stoichiometric matrix $\St$ and reaction $r \in \R_u$. For each $r$, we denote the forward and reverse rate constants as $k_r^+$ and $k_r^-$, respectively, forward and reverse reaction-flow vectors as

\[j_r^+ = k_r^+ q^{r^-} \text{ and  } j_r^- = k_r^- q^{r^+},\]

respectively, the net reaction-flow vector through $r$ as

\[j_r = j_r^+ - j_r^-,\]

and the resulting concentration-flow vector as

\begin{equation}
\dot{q} = \sum_{r\in \R_u} (\St)^r j_r = \sum_{r\in \R_u} (\St)^r (j_r^+-j_r^-). \label{eq:reversible_MAK}
\end{equation}

In this work, for reversible CRNs, we use the following strong version of the theorem due to Horn, Jackson, and Feinberg (Theorem 7.1.1, [61]).

A reversible CRN of deficiency zero taken under MAK, regardless of the rate constants, has precisely one stable equilibrium within each stoichiometric compatibility class.

Furthermore, it can be shown that the equilibrium for deficiency zero reversible CRNs is either detailed-balanced (Ch. 14, [61]), where the forward and backward reaction-flow for each reaction are identical

\begin{align}
j_r(\eq) &= 0, \nonumber\\
j_r^+(\eq) &= j_r^-(\eq), \nonumber\\
k_r^+ \eq^{r^-} &= k_r^- \eq^{r^+},  \label{eq:detailed_balanced}
\end{align}

or complex-balanced (see Ch. 15, [61] for details).

We define a class of CRNs, parametrized by $g \in \mathbb{Z}_{> 0} $, that we refer to as the complete generalized-Schlogl model and denote as $\mathcal{M}'_g$. Taken under MAK, $\mathcal{M}'_g = \{\mathcal{S}'_g,\mathcal{C}'_g,\mathcal{R}'_g,\mathcal{K}'_g\}$ where

\begin{align*}
\mathcal{S}'_g &= \{ X, Y_0, Y_2, Y_3, \ldots, Y_{2g+1}\},\\
\mathcal{C}'_g &= \{ X, Y_0, Y_2 + Y_0 + 2X, 3X + Y_3,  \ldots,\\
& \phantom{=0} Y_{2g} + Y_0 + (2g)X, (2g+1)X + Y_{2g+1}\},\\
\{\mathcal{R}'_g,\mathcal{K}'_g\} &= \{\ce{Y_0 <=>[$k_0'$][$k_1'$] X},\\
& \phantom{=0} \ce{Y_2 + Y_0 + 2X <=>[$k_2'$][$k_3'$] 3X + Y_3},\\
& \phantom{=0} \ldots \\
& \phantom{=0}  \ce{Y_{2g} + Y_0 + (2g)X <=>[$k_{2g}'$][$k_{2g+1}'$] (2g + 1)X + Y_{2g+1}}\}.
\end{align*}

Let us denote the stoichiometric matrix for $\mathcal{M}'_g$ as $\mathbb{S}'_g$. Then for $g=2$,

\[  \mathbb{S}'_2 =
\begin{bmatrix}
\phantom{-}1 & -1 & \phantom{-}1 & -1 & \phantom{-}1 & -1 \\
-1 & \phantom{-}1 & -1 & \phantom{-}1 & -1 & \phantom{-}1\\
\phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}1 &\phantom{-}0 & \phantom{-}0\\
\phantom{-}0 & \phantom{-}0 & \phantom{-}1 & -1 &\phantom{-}0 & \phantom{-}0\\
\phantom{-}0 & \phantom{-}0 &\phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}1 \\
\phantom{-}0 & \phantom{-}0 &\phantom{-}0 & \phantom{-}0 & \phantom{-}1 & -1 \\
\end{bmatrix}
\]

For any $g$, we have

\begin{align*}
|\mathcal{S}'| & = 2g + 2,\\
|\mathcal{R}'| & = 2g +2,\\
\iota &= g+1,\\
\ell & = g+1.
\end{align*}

$\iota$ can be counted by observing that any null cycle is given by going forward and backward on the same reaction vector (and these are the only null cycles). $\ell$ can be counted by observing that the pairs of species with the same conserved quantity are $\{X, Y_0\}$, $\{Y_{2i},Y_{2i+1}\}_{i \in [1,g]}$ (and these are the only conservation laws).
This yields that for any $g$, the deficiency of $\mathcal{M}_g'$ is zero, $\delta (\mathcal{M}_g') = 0$. Thus, applying Theorem <ref>, every CRN in the complete generalized-Schlögl model class has exactly one asymptotically stable equilibrium irrespective of the choice of rate constants.

We define another class of CRNs, parametrized by $g \in \mathbb{Z}_{> 0} $, that we refer to as the generalized-Schlogl model [63] and denote as $\mathcal{M}_g$. Taken under MAK, $\mathcal{M}_g = \{\mathcal{S}_g,\mathcal{C}_g,\mathcal{R}_g,\mathcal{K}_g\}$ where

\begin{align*}
\mathcal{S}_g &= \{ X\},\\
\mathcal{C}_g &= \{\vn, X, 2X, \ldots,(2g+1)X\},\\
\{\mathcal{R}_g,\mathcal{K}_g\} &= \{\ce{$\vn$ <=>[$k_0$][$k_1$] X},\\
& \phantom{=0} \ce{2X <=>[$k_2$][$k_3$] 3X },\\
& \phantom{=0} \ldots \\
& \phantom{=0}  \ce{(2g)X <=>[$k_{2g}$][$k_{2g+1}$] (2g + 1)X}\}.
\end{align*}

Let us denote the stoichiometric matrix for $\mathcal{M}_g$ as $\mathbb{S}_g$. Then for $g=2$,

\[  \mathbb{S}_2 =
\begin{bmatrix}
\phantom{-}1 & -1 & \phantom{-}1 & -1 & \phantom{-}1 & -1
\end{bmatrix}
\]

For any $g$, we have

\begin{align*}
|\mathcal{S}| & = 1,\\
|\mathcal{R}'| & = 2g +2,\\
\iota &= g+1,\\
\ell & = 0.
\end{align*}

$\iota$ can be counted by observing that any null cycle is given by going forward and backward on the same reaction vector (and these are the only null cycles). $\ell = 0$ as there are no conservation laws.
This yields that for any $g$, the deficiency of $\mathcal{M}_g$ is $g$, $\delta (\mathcal{M}_g') = g$. Thus, we cannot apply Theorem <ref> and this CRN may exhibit multiple equilibria.

§.§ Stochastic chemical kinetics

The discrete nature of CRNs is modeled stochastically as a pure jump Markov process where a state is given by a species population vector ($n \in \mathbb{Z}^{\Spc}_{\geq 0}$) [64]. A reaction, when it occurs, replaces an input combination of species sampled without-replacement from the system into an output combination and the system jumps to a new state. An object of study for stochastic CRNs is the probability distribution $\mathbb{P}(n,t)$ of finding the system in state $n$ at time $t$. It is well-known that this probability distribution evolves under the chemical master equation [65] and the dynamics of its first moment (average) reduces to the MAK in systems with a very large number of particles [66].

Denoting the volume of the system by $V$, consider the stochastic dynamics of the process $q^V_t = n(t)/V$ describing the concentration of species at time $t$ starting from some concentration $q^V_0$.
From path-integral methods (Ch. 4, [3]; Ch. 4, [67]), it is known that the probability of observing $q^V$ at time $t$ conditioned on the system being in $q^V_0$ at time $0$ is asymptotically given by the large-deviation (LD) principle of the form (Eq. 226, [68])

\begin{equation}
\mathbb{P}(q_t^V|q^V_0) \asymp e^{-V S(q_t,q_0)},
\end{equation}

where the rate function $S(q_t,q_0)$ is evaluated along the optimal trajectory with appropriate boundary conditions

\[ S(q_t,q_0) = \inf_{q(t):q(0)=q_0,q(t)=q_t} J[q]\]

of the functional

\begin{equation}
J[q] = \sup_{p(t)} \int_0^t \,dt \left( p \cdot \dv{q}{t} - H(p,q) \right). \label{eq:J_rate}
\end{equation}

Here, $H(p,q)$ is the Hamiltonian of the process and for CRNs, it takes the form

\begin{align}
H(p,q) &= \sum_{r \in \mathcal{R}} (e^{(r^+-r^-)\cdot p}-1)k_{r}q^{r^-}, \nonumber\\
&= \sum_{r \in \mathcal{R}} (e^{\Delta r \cdot p}-1)k_{r}q^{r^-},\label{eq:Ham-kinetics}
\end{align}

where $p \in \mathbb{R}^\mathcal{S}$ is the conjugate momentum. For more on variational methods in LD theory and LD theory for CRNs, see App. C in [69] and [70], respectively.

The rate function can be re-expressed as value along the saddle-point solutions

\[S(q_t,q_0) = \inf_{q(t):q(0)=q_0,q(t)=q_t} \sup_{p(t)} \mathcal{A}[p,q] \]

of the action functional

\begin{equation}
\mathcal{A}[p,q] = \int \,dt \left( p \dv{q}{t} - H(p,q) \right). \label{eq:action_func}
\end{equation}

It is then straightforward to show that the rate function satisfies the Hamilton-Jacobi PDEs  [71]

\[ \pdv{S}{t} = -H\left(\pdv{S}{q_t},q_t\right).\]

The stationary distribution of a stochastic process, if it exists, is the probability distribution over states (or its corresponding rate function) that does not change with time. It can be obtained from stochastic simulations [72, 23] by simulating the system for a very long time and normalizing the residence times of the system in (a small neighbourhood of) each state. Since the probability of a state $q_t$ in the stationary distribution is independent of the initial condition and time, the stationary distribution is only a function of the state where evaluated and we denote it as $\pi(q)$. Denote the rate function of the stationary distribution, also called the non-equilibrium potential (NEP) [59] by $\V(q)$. Thus, the NEP must satisfy

\begin{equation}
H\left(\pdv{\V}{q},q\right) = 0. \label{eq:NEP}
\end{equation}

Thus, the NEP is calculated by finding solutions to the Hamilton-Jacobi equations in the $H=0$ submanifold (Sec. 4.10, [67]).

It is well-known that the integral curves to the Hamilton-Jacobi PDEs are the solutions of Hamilton's equations of motion (EoM)

\begin{align}
\dv{q}{t} &= \phantom{-}\pdv{H(p,q)}{p}, \nonumber\\
\dv{p}{t} &= -\pdv{H(p,q)}{q}.
\label{eq:HamEoM}
\end{align}

Since, in this work, we are only interested in understanding the stationary distribution, we only investigate the classes of solutions of Hamilton's EoM in the $H=0$ submanifold, namely: relaxation and escape trajectories [73].

The equations of MAK (compare with Eq. <ref>) are obtained as the $p=0$ solution

\begin{align}
\dot{q}_\text{rel} := \dv{q}{t}\bigg|_{p=0} &= \pdv{H}{p}\bigg|_{p=0}\nonumber\\
&= \sum_r (r^+-r^-)k_{r}q^{r^-}, \nonumber\\
&= \sum_r \St^r k_{r}q^{r^-} \label{eq:relax_Ham}
\end{align}

and their solutions are called relaxation trajectories of the system.

The fixed points (or equilibrium) of the relaxation trajectories, when they exist, inform about the long-term behavior dynamics of CRNs. Denoting an arbitrary fixed point as $\eq$ and their set as $\{\eq\}$, we have $\dot{q}_\text{rel}(\eq)=0.$
The stability of a fixed-point can be determined by evaluating the eigenvalues of the Jacobian matrix $J$, where

\[J = \pdv{\dot{q}_\text{rel}}{q} = \pdv{H}{q}{p}\bigg|_{p=0},\]

at that point. In particular, we term fixed-points where all non-zero eigenvalues have negative real parts as attractors, all but one non-zero eigenvalues have negative real parts as transient points, and all other fixed-points as unstable points. CRNs can also exhibit periodic or chaotic attractors, but in this work we only consider CRNs with (point) attractors as defined above. A system with multiple attractors is called multistable.

The solutions with $p \neq 0$ in the $H=0$ submanifold are termed as the escape trajectories. These are trajectories in the phase space (concentrations and their conjugate momenta) and we denote them as $(q_\text{esc},p_\text{esc},t)$. While any two points $q_1$ and $q_2$ may not be joined together by a single escape trajectory, it is possible to join them with a concatenation of escape trajectories with the time variable possibly running against the direction of the curve parametrization (Sec. 3.2.2, [73]), and we denote it as $(q^*,p^*,t^*)$ where $q^*(0) = q_1$ and $q^*(\max(t^*)) = q_2$. Using Eqs. <ref> and <ref>, the difference in NEP between the two points is given as

\begin{align}
\mathcal{V}(q_2) - \mathcal{V}(q_1) &= \int_{q^*} \pdv{\V}{q} \cdot \,dq = \int_{q^*} p^*(q) \cdot \,dq. \label{eq:NEP_diff}
\end{align}

For an information-geometric explanation of the role of escape trajectories in estimating the NEP, see Sec. IV A 5 [74].

A transient point is adjacent to an attractor if there is a relaxation trajectory that starts in the vicinity of the transient point and terminates at the attractor. By the mountain pass theorem [75, 76], it is known that any two nearby attractors must have at least one transient point in between them. Two attractors adjacent to the same transient point will be called adjacent to each other. A multistable stochastic CRN spends indefinite amount of time in the vicinity of an attractor before transitioning out of it due to fluctuations and entering another attractor adjacent to it. For reversible CRNs, wherein for every reaction, its reverse is also possible, there must exist trajectories, particularly an optimal or escape trajectory, that emanate out of the attractor and terminate at the transient point. Estimating these escape trajectories can be computationally challenging and can be carried out by using the shooting-method [77] or functional gradient descent algorithms as proposed in [18] or [78]. Once determined, they can then be used in Eq. <ref> to estimate the difference in the NEP at the attractors (Sec. 4.2.2, [73]).

Suppose a transient point $t$ has two adjacent attractors $a_1$ and $a_2$. The following theorem (Lemma 1, [73]) proves that, for a multistable CRN, the NEP at attractors is lesser or equal than that at the adjacent transient point

\[ \V(t) \geq \V(a_1), \hspace{0.2em} \V(t) \geq \V(a_2).\]

The NEP is a Lyapunov function along the relaxation trajectories given by MAK.

First, recall that the equations of MAK are the gradient of the Hamiltonian in the momentum variable evaluated at $p=0$,

\[\dot{q}_\text{rel} = \pdv{H}{p}\bigg|_{p=0}.\]

Moreover, the Hamiltonian is identically zero at $p=0$

\[H(0,q) = 0.\]

Secondly, notice that the Hessian of the Hamiltonian in the momentum variables is positive-definite,

$$\forall v \in \R^{\Spc}, \quad v^T\pdv{H}{p^T}{p}v \geq 0,$$

or equivalently the Hamiltonian is convex in the momentum variables. This means that for any concentration $q$,

\[H(p,q) \geq H(0,q) + p \cdot \pdv{H}{p}\bigg|_{p=0}.\]

Since the escape momentum $p_\text{esc}$ lies in the $H=0$ submanifold, we have

\[H(p_\text{esc},q) = 0 \geq 0 + p_\text{esc} \cdot \pdv{H}{p}\bigg|_{p=0}, \]

which yields

\[p_\text{esc} \cdot \pdv{H}{p}\bigg|_{p=0} \leq 0 \hspace{0.25em} \forall q.\]

Finally, using Eqs. <ref> and <ref>, we get

\[\pdv{\mathcal{V}}{q} \cdot \dot{q}_\text{rel} \leq 0,\]

which concludes the proof.

Stochastic simulation

Model 1

Model 2

Model 3

Stochastic simulation using Gillespie algorithm, non-equilibrium potential (NEP), and mass-action kinetic (MAK) polynomial is showed for the three models from Example <ref>. The MAK polynomial is arbitrarily scaled to overlay on the NEP to indicate the positions of the fixed-points and their stability.

For a CRN with a detailed balanced (Eq. <ref>) or complex-balanced equilibrium, it can be shown that (see App. B, [18]):

\begin{align}
\pdv{\V}{q} = p_\text{esc} &= \ln\left(\frac{q}{\eq}\right), \nonumber\\
\V(q) - \V(\eq)  &= \int_{\eq}^q  \ln\left(\frac{q}{\eq}\right) \cdot \,dq \nonumber\\
&= q\ln\left(\frac{q}{\eq}\right) - (q-\eq). \label{eq:NEP_comp_balanced}
\end{align}

This a well-known result for complex-balanced zero deficiency CRNs and the NEP $\V$ is also called the Horn-Jackson potential [79, 80]. Since the complete generalized Schlögl model from e.g. <ref> is a deficiency zero reversible CRN and has a unique detailed-balanced or complex-balanced equilibrium for any choice of rate constants, its NEP must be of the form <ref>.

Consider the generalized Schlögl model from e.g. <ref> with one species and non-zero deficiency. Observe that there are only two distinct reaction vectors in the CRN $\Delta r = \pm 1$. Every reaction with $\Delta r=+1$ or $-1$ can be seen as a birth or death process, respectively, and the CRN is a birth-death process [59].

Define polynomials $f_\pm(q)$ (suppressing the dependency on rate constants) with the rates of birth and death as follows:

\begin{align*}
f_+(q) &= \sum_{i=0}^{g} k_{2i} q^{2i},\\
f_-(q) &= \sum_{i=0}^{g} k_{2i+1} q^{2i+1}.
\end{align*}

Using Eq. <ref>, the Hamiltonian for $\mathcal{M}_g$ is

\[H(p,q) = (e^p-1)f_+(q) + (e^{-p}-1)f_-(q). \]

Using Eq. <ref>, the relaxation trajectory or MAK equation is

\[\dot{q}_\text{rel} = f_+(q)-f_-(q) =\sum_{i = 0}^{2g +1} (-1)^i k_i x^i. \]

The escape solution for the system is obtained as the $p \neq 0$ solution along the $H=0$ submanifold. Since the phase space is two-dimensional, one constraint uniquely identifies $p(q) \neq 0$ in the $H=0$ submanifold, yielding

\[p_\text{esc} = \log{\left(\frac{f_-(q)}{f_+(q)}\right)}.\]

One can verify that $H(p_\text{esc},q)=0$ for all $q$ and the this solution satisfies Eq. <ref>.
Since this is a one-dimensional system, the escape trajectories in concentration space are trivial, and the difference of the NEP between any two points is given by

\begin{align*}
\mathcal{V}(q_2) - \mathcal{V}(q_1) &= \int_{q_1}^{q_2} \log{\left(\frac{f_-(q)}{f_+(q)}\right)} \,dq.
\end{align*}

For this model, as an application of Viéte's formula [81], there is a bijection between rate constants and fixed points.
Suppose the rate constants were such that $\dot{q}_\text{rel}$ had $2g+1$ positive real fixed points,

\[ 0\leq \eq_1 \leq \ldots \leq \eq_{2g+1}.\]

Then, we can also express the relaxation trajectory as

\[\dot{q}_\text{rel} = - \prod_{i=0}^{2g+1} (q - \eq_i).\]

For concreteness, we consider three models $\mathcal{M}^{(i)}_2, i \in \{1,2,3\}$ with $g=2$. The set of rate constants $\mathcal{K} = \{k_0,k_1,k_2,k_3,k_4,k_5\}$ are given as:

\begin{align*}
\mathcal{K}^{(1)} &= \{124.2, 286.44, 236.72, 88.88, 15.4, 1\},\\
\mathcal{K}^{(2)} &= \{84.6, 218.22, 201.46, 81.74, 14.9, 1\},\\
\mathcal{K}^{(3)} &= \{59.4, 158.58, 154.14, 67.46, 13.5, 1\}.
\end{align*}

The fixed-points for these models are:

\begin{align*}
\{\eq\}^{(1)} &= \{ 1,1.8,3,4.6,5\},\\
\{\eq\}^{(2)} &= \{ 1,1.2,3,4.7,5\},\\
\{\eq\}^{(3)} &= \{ 1,1.8,3,3.3,5\}.
\end{align*}

A stochastic run for each of these models and their NEPs overlaid with (arbitrarily scaled) MAK polynomial can be found in Fig. <ref>. Notice that the relative positions of the roots determine the phase of the dynamics and affect of the adjacent attractors has a higher residence time.

§.§ Thermodynamics of CRNs

For reversible CRNs, the entropy production rate (EPR) $\sigma$ is defined as (Eq. 55, [82])

\begin{equation}
\sigma(j) =  \sum_{r\in \R_u} \ln\left(\frac{j^+_r}{j^-_r}\right)j_r, \label{eq:EPR}
\end{equation}

where $\R_u, j^+_r, j^-_r$ are defined above Thm. <ref> in App. <ref>. In this subsection, we will restrict to reversible CRNs and henceforth use $\R$ to mean $\R_u$. Observe that each term of the type $\log(a/b)(a-b)$ is non-negative, and thus $\sigma(j) \leq 0$ for all $j$. In fact, $\sigma(j)=0$ only at detailed-balanced equilibrium where $j^+_r = j^-_r$.

A CRN will be said to be thermodynamically consistent if each species can be assigned a chemical potential of formation $\mu_f \in \mathbb{R}^\mathcal{S}$ such that each reaction $r \in \R_u$ can be assigned a barrier height $E^b_r$ and the rate constants are of the form

\begin{align}
k_r^+ &= e^{-(E^b_r - \mu_f \cdot r^-)/(k_BT)}, \nonumber\\
k_r^- &= e^{-(E^b_r - \mu_f \cdot r^+)/(k_BT)}, \nonumber
\end{align}

where $k_B$ and $T$ are the Boltzmann constant and temperature, respectively. In what follows, we choose units where $k_B T = 1$. We can re-express the condition in terms of the ratio of $K = k_r^+/k_r^-$, also called the equilibrium constant, as

\begin{equation}
\ln K_r := \ln\left( \frac{k^+_r}{k^-_r}\right) = -\mu_f \cdot (r^+-r^-). \label{eq:equilibrium_constant}
\end{equation}

It is known from a generalization of Wegscheider's conditions [16] that $\ln K \in \text{Im}(\St^T)$, as in the above equation, if and only if there is a detailed balanced equilibrium $\eq$ that satisfies Eq. <ref>

\[  K_r = \frac{\eq^{r^+}}{\eq^{r^-}}. \]

On comparing the above two equations, we find that

\[ \St^T (\mu_f + \ln \eq) = 0.\]

Define the equilibrium chemical potential
$$\mu_0 := \mu_f + \ln \eq,$$
and the chemical potential

\[ \mu :=  \mu_f + \ln(q) = \mu_0 + \ln\left(\frac{q}{\eq}\right). \]

This yields the affinity

\begin{align}
\ln\left(\frac{j^+_r}{j^-_r}\right) &= -(\mu_f+\ln(q)) \cdot (r^+-r^-) \nonumber\\
&= - \mu \cdot (r^+-r^-), \label{eq:affinity}
\end{align}

and the EPR (Eq. <ref>) becomes

\begin{equation}
\sigma = -\mu^T \St j = -\ln\left(\frac{q}{\eq}\right) \cdot \dot{q}, \label{eq:balanced_EPR}
\end{equation}

where we use $\mu_0 \St = 0$ and $\dot{q}$ to refer to the MAK flow from Eq. <ref>. The chemical potential determines the gain in free energy to the system by adding a particle of the species at a concentration (Sec. III C, [82]). Thus, the affinity measures the loss in free energy due to one occurrence of a reaction at that concentration, and the EPR is rate of loss in free energy for the complete CRN. Define the entropy produced (EP) as the integral of the EPR in time

\begin{equation}
\Sigma = \int \,dt \sigma. \label{eq:EP}
\end{equation}

Observe that the EP is negative of the the Horn-Jackson potential or the NEP for detailed balanced CRNs (Eq. <ref>). Thus, from the Lyapunov property of NEP, the EP by a thermodynamically consistent CRN never decreases as the system evolves deterministically.

It is known that CRNs with a detailed-balanced or complex-balanced equilibrium have precisely one equilibrium in each stoichiometric compatibility class (Theorems 14.2.3 and 15.2.2, [61]). While Feinberg's characterization of quasi-thermodynamic kinetic systems (Sec. 13.4, [61]) admits both detailed-balanced and complex-balanced systems, only CRNs with a detailed-balanced equilibrium admit a thermodynamic description as described above. Eq. <ref> is used as a starting point  in Sec. III, [82] to define the Gibb's free energy of reaction

\begin{equation}
\Delta_r G := \mu^T \mathbb{S}^r =  -\ln\left(\frac{j^+_r}{j^-_r}\right), \label{eq:Gibbs_reaction}
\end{equation}

also called the thermodynamic force driving reaction $r$, to define the thermodynamics of CRNs. For reaction flow $j_r \in \mathbb{R}$, the EPR per reaction $r$ is

\begin{equation}
\sigma_r(j_r) = - \Delta_r G j_r, \label{eq:EPR_reaction}
\end{equation}

and the EPR of the CRN is given by their sum over all reactions

\begin{equation}
\sigma = \sum_r \sigma_r. \label{eq:EPR_CRN_2}
\end{equation}

For various thermodynamic decompositions of EPR for complex-balanced and multistable CRNs, see [83].

The considerations above demonstrate that in a (thermodynamically consistent) closed system, where the dynamics of all species is modeled by MAK, the concentration will always approach a detailed-balanced equilibrium. However, in open systems, where the dynamics of only a species subset is modeled by MAK whereas the others are controlled by the environment, the system can be driven towards and maintained at a non-equilibrium steady state (NESS).

Henceforth, we follow terminology and notation from Sec. <ref>.
Assuming that there are no reactions in the complete CRN not in $\R$, i.e. $\notR = \vn$, the complete stoichiometric matrix

\[ \nabla = \begin{bmatrix}
\nabla^{\R}_\mcl{E} \\
\St
\end{bmatrix}.\]

If the open system is allowed to evolve for a sufficiently long time, it will relax to a NESS assuming that the internal CRN has point attractors. The species-flow through the boundary species at the NESS is not necessarily zero and is the source of thermodynamic work done by the environment to maintain the NESS. Any non-vanishing boundary species-flow at the NESS is called an emergent cycle [84], denoted by $\epsilon \in \mathbb{R}^{\R}$, and satisfies

\begin{equation*}
\nabla \epsilon = \begin{bmatrix}
\nabla^{\R}_\mcl{E} \epsilon \\
\end{bmatrix}.
\end{equation*}

The basis of emergent cycles (Sec. V, [85])

$$\{\epsilon\} = \text{Ker}(\St)/\text{Ker}(\nabla)$$

consists of all reaction-flows in the kernel of the internal stoichiometric matrix that are not in the kernel of the embedding stoichiometric matrix. Observe that the EPR due to $\epsilon$ only depends on the chemical potential of the external-boundary species (Sec. II D, [74])

\begin{equation}
\sigma(\epsilon) = -\mu^T \nabla \epsilon = - \sum_{s \in \mcl{E}} \mu_s (\nabla^{\R}_\mcl{E} \epsilon)_s. \label{eq:EPR_emergent}
\end{equation}

From physical considerations, it is useful to decompose the complete CRN $\mcl{G}'$ through a basis of processes [17]
$$p\in \mathcal{P}: p^- \to p^+ ,$$

that respect the conservation laws and span the stoichiometric subspace.
Thus, analogously to the stoichiometric matrix, define the process matrix $\mathbb{P}$ such that

\begin{align*}
(\mathbb{P})^p &= p^+-p^-,\\
\text{Ker}(\mathbb{P}^T) &= \text{Ker}(\nabla^T),\\
\text{Im}(\mathbb{P}) &= \text{Im}(\nabla).
\end{align*}

To transform from the reaction basis to process basis, we define a $\mcl{P} \times \R$ process-reaction matrix $\mbf{C}$ such that

\[ \sum_p (\mathbb{P})^p_s \mbf{C}_p^r = \nabla^r_s \hspace{1em} \forall s \in \Spc'. \]

Using unprimed variables for reactions and primed variables for processes, analogous to Eq. <ref>, we define Gibb's free energy of process

\[\Delta'_p G := \mu^T (\mathbb{P})^p. \]

For reaction flow $j \in \mathbb{R}^{\R}$, the flux through process $p$ $j'p \in \mathbb{R}^\mathcal{P}$ is

\[j'_p = \sum_r \mbf{C}_p^r j_r.\]

The EPR for a process for a given reaction-flow $j$

\begin{align}
\sigma'_p(j) &=  - \mu^T (\mathbb{P})^p (\mbf{C})_p j, \nonumber \\
&= -\Delta'_p G  j'_p, \label{eq:process_EPR}
\end{align}

and the total EPR (analogous to Eq. <ref>) is

\begin{equation}
\sigma = \sum_{p \in \mathcal{P}} \sigma'_p. \label{eq:EPR_CRN_3}
\end{equation}

From the definition of EPR in Eq. <ref>, we know each term is individually nonnegative or a reaction always dissipates chemical potential energy. However, under a process decomposition of a reaction, the EPR of a process, $\sigma'_p$ in Eq. <ref>, can be negative in which case the process will be said to be doing chemical work or transducing free energy at the expense of other processes which dissipate chemical potential at a higher rate to ensure that the total EPR in Eq. <ref> is positive. The transduction efficiency of the processes is defined as (Sec. V, [17])

\begin{equation}
\eta = - \frac{\text{Output}}{\text{Input}}, \label{eq:efficiency}
\end{equation}

\begin{align*}
\text{Output} &= \left\{\sum_p \sigma'_p \bigg| \sigma'_p < 0 \right\}, \\
\text{Input} &= \left\{\sum_p \sigma'_p \bigg| \sigma'_p > 0 \right\}. \\
\end{align*}

Using Eqs. <ref> and <ref>, we can calculate the efficiency of the environment in maintaining a NESS at emergent cycle $\epsilon$ with EPR given by Eq. <ref>.

The transduction efficiency at the fixed points and EPR for all concentrations of the 3 models from e.g. <ref> are shown in the left and right panel, respectively. For these models, the attractors are fixed at {1,2,3} and the transient points in between vary. While the NEP profile for these models are different (Fig. <ref>), the transduction efficiency and EPR monotonically increase in the distance of the fixed points from the origin.

Consider the complete generalized-Schlogl model from e.g. <ref> as the embedding CRN. Since it is a reversible zero deficiency CRN and every linkage class has a single reaction, the equilibrium must be detailed-balanced and hence the embedding CRN is thermodynamically consistent. Let the boundary set
$$\Spc_B = \{Y_0, Y_2, \ldots , Y_{2g + 1}\}$$
chemostatted at concentrations $q^*$ by the environment.
The resulting internal CRN then is the generalized Schlogl model from e.g. <ref> with rate constants

\begin{align*}
k_0 &= k_0' q^*_{Y_0}, \\
k_1 &= k_1',\\
k_{2j} &= k_{2j}' q^*_{Y_0} q^*_{Y_{2j}}, \\
k_{2j+1} &= k_{2j+1}' q^*_{Y_{2j+1}},
\end{align*}

for $j \in [1,g].$

Recall that the unique reactions in the embedding CRN are

\begin{align*}
\mathcal{R} &= \{r_0: \ce{Y_0 -> X},\\
& \phantom{=0} r_1: \ce{Y_2 + Y_0 + 2X -> 3X + Y_3},\\
& \phantom{=0} \ldots \\
& \phantom{=0} r_g: \ce{Y_{2g} + Y_0 + (2g)X -> (2g + 1)X + Y_{2g+1}}\}.
\end{align*}

For this system, we define the set of processes $\mcl{P}$ to consist of

\[p_0 := \ce{Y_0 -> X}, \hspace{1em} p_j := \ce{Y_{2j} -> Y_{2j+1}},\]

for $j \in [1,g].$
Observe that the set respects the conservation laws and spans the stoichiometric subspace of the embedding CRN. It is easy to see that

\begin{equation}
r_0 = p_0,\hspace{1em} r_j = p_j + p_0 \label{eq:process_decomp_eg}
\end{equation}

for $j \in [1,g].$ This can be used to obtain the process-reaction matrix $\mbf{C}$ as defined earlier in the subsection.

To write each reaction's change in Gibb's energy and EPR, and decompose it using processes, observe that we can write the MAK as a sum over reaction pairs as

\begin{align}
\dot{q} &= \sum_{j=0}^g  k_{2j+1} \left(\ell_j - q \right)q^{2j}, \label{eq:MAK_ellj}
\end{align}

\[  \ell_j := \frac{k_{2j}}{k_{2j+1}} \text{ for } j \in [0,g].\]

Using Eq. <ref>, the change in Gibb's energy for the $j^\text{th}$ reaction is

\begin{align*}
\Delta_j G &= \ln \left( \frac{q}{\ell_j}\right).
\end{align*}

The above equation, along with the decomposition of reactions into processes from Eq. <ref>, yields

\begin{align*}
\Delta_0' G &= \Delta_0 G = \ln \left( \frac{q}{\ell_0}\right),\\
\Delta_j' G &= \Delta_j G - \Delta_0' G = \ln \left( \frac{\ell_0}{\ell_j}\right),
\end{align*}

for $j \in [1,g]$.
For a choice of rate constants where all fixed points are real, using [86], the sequence of the coefficients $k_i$ is log-concave, or

\[ k_i^2 \geq k_{i-1} k_{i+1}.\]

This yields,

\begin{align*}
k_i^2 k_{i-1}^2 &\geq k_{i-1} k_{i+1} k_{i-2} k_{i}\\
k_i k_{i-1} &\geq k_{i+1} k_{i-2} \\
\ell_{i} & \geq \ell_{i-1},
\end{align*}

and, in particular,

\begin{equation}
\ell_0 \leq \ldots \leq \ell_g.
\end{equation}

\[ \Delta_j' G  < 0, \text{ for } j \in [1,g].\]

The set of NESSs to which the open system can relax are given by the fixed-points of the internal CRN.
Let the set of fixed points $\{\eq\}$ of the MAK of the internal CRN in species $X$ satisfy

\[ 0 \leq \eq_1 \leq \ldots \leq \eq_{2g+1},\]

such that

\[ \dot{q} = - \prod_{i=0}^{2g+1} (q - \eq_i).\]

At any of the fixed points, the current through process $0$ that involves $X$ must be zero ($j'_0=0$), and only the remaining processes produce entropy. Using Eq. <ref>, the EPR through process $j \in [1,g]$ at a fixed point $\eq$ is

\[
\sigma_j'(\eq) =  -\left(\Delta'_j G\right)  k_{2j+1} \left(\ell_j - \eq \right)\eq^{2j}.
\]

As shown above, $\Delta'_j G$ is always negative for all $j$. For a process, different regimes of entropy production and free-energy transduction will be achieved when $\eq$ crosses $\ell_j$. In particular, process $j$ produces entropy when $\eq < \ell_j$ and transduces free energy when $\eq > \ell_j$. This means that more processes transduce free energy rather than produce entropy as we move from a fixed point closer from the origin to one that is further. Thus, we conjecture that the model's transduction efficiency at NESSs is monotonic in their distance from the origin.

The transduction efficiency at the fixed-points of the 3 models from e.g. <ref> are shown in the left panel of Fig. <ref>. Recall that the attractors for this model are at $\{1,3,5\}.$ Observe that the transduction efficiency is higher at a fixed point further away from the origin, and we conclude that for these models, the attractor further from the origin transduces free energy more efficiently. Also, from Eq. <ref>, the total EPR at any concentration $q$ is (see right panel in Fig. <ref>)

\[\sigma(q) = \sum_{j=0}^g -\ln\left(\frac{q}{\ell_j}\right)  k_{2j+1} \left(\ell_j - q \right)q^{2j}.\]

Observe that the entropy produced also increases monotonically at the fixed points in their distance from the origin, and we conjecture that the EPR at fixed points is monotonic in the distance of the fixed points from the origin for any choice of rate constants for this model. The implications of this conjecture are discussed in Sec. <ref>.

[1]
Andrei D Linde.
Phase transitions in gauge theories and cosmology.
Reports on Progress in Physics, 42(3):389, 1979.

[2]
Alexander Altland and Ben D Simons.
Condensed matter field theory.
Cambridge university press, 2010.

[3]
Eric Smith and Supriya Krishnamurthy.
Symmetry and collective fluctuations in evolutionary games.
IOP Publishing, 2015.

[4]
Cole Mathis, Tanmoy Bhattacharya, and Sara Imari Walker.
The emergence of life as a first-order phase transition.
Astrobiology, 17(3):266–276, 2017.

[5]
Martin A Nowak and Hisashi Ohtsuki.
Prevolutionary dynamics and the origin of evolution.
Proceedings of the National Academy of Sciences, 105(39):14924–14927, 2008.

[6]
Stuart Kauffman and Andrea Roli.
Is the emergence of life an expected phase transition in the evolving universe?
arXiv preprint arXiv:2401.09514, 2024.

[7]
Freeman J Dyson.
A model for the origin of life.
Journal of Molecular Evolution, 18:344–350, 1982.

[8]
James D Watson and Francis HC Crick.
The structure of dna.
In Cold Spring Harbor symposia on quantitative biology, volume 18, pages 123–131. Cold Spring Harbor Laboratory Press, 1953.

[9]
Eric Smith and Harold J Morowitz.
The origin and nature of life on earth: the emergence of the fourth geosphere.
Cambridge University Press, 2016.

[10]
Manoj Gopalkrishnan.
Catalysis in reaction networks.
Bulletin of mathematical biology, 73(12):2962–2982, 2011.

[11]
Hong Qian.
Open-system nonequilibrium steady state: statistical thermodynamics, fluctuations, and chemical oscillations, 2006.

[12]
David Andrieux and Pierre Gaspard.
Nonequilibrium generation of information in copolymerization processes.
Proceedings of the National Academy of Sciences, 105(28):9516–9521, 2008.

[13]
Pierre Gaspard.
Kinetics and thermodynamics of living copolymerization processes.
Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374(2080):20160147, 2016.

[14]
Praful Gagrani, Victor Blanco, Eric Smith, and David Baum.
Polyhedral geometry and combinatorics of an autocatalytic ecosystem.
Journal of Mathematical Chemistry, pages 1–67, 2024.

[15]
Vinay Ambekar Ranganath and Indrajit Maity.
Artificial homeostasis systems based on feedback reaction networks: Design principles and future promises.
Angewandte Chemie International Edition, page e202318134, 2024.

[16]
Stefan Schuster and Ronny Schuster.
A generalization of wegscheider's condition. implications for properties of steady states and for quasi-steady-state approximation.
Journal of Mathematical Chemistry, 3(1):25–42, 1989.

[17]
Artur Wachtel, Riccardo Rao, and Massimiliano Esposito.
Free-energy transduction in chemical reaction networks: From enzymes to metabolism.
The Journal of Chemical Physics, 157(2), 2022.

[18]
Praful Gagrani and Eric Smith.
Action functional gradient descent algorithm for estimating escape paths in stochastic chemical reaction networks.
Physical Review E, 107(3):034305, 2023.

[19]
Hugo Touchette.
The large deviation approach to statistical mechanics.
Physics Reports, 478(1-3):1–69, 2009.

[20]
David F Anderson, Gheorghe Craciun, and Thomas G Kurtz.
Product-form stationary distributions for deficiency zero chemical reaction networks.
Bulletin of mathematical biology, 72:1947–1970, 2010.

[21]
Linard Hoessly.
Stationary distributions via decomposition of stochastic reaction networks.
Journal of Mathematical Biology, 82(7):67, 2021.

[22]
Frank P Kelly.
Reversibility and stochastic networks.
Cambridge University Press, 2011.

[23]
Daniel T Gillespie.
Stochastic simulation of chemical kinetics.
Annu. Rev. Phys. Chem., 58:35–55, 2007.

[24]
Hector Zenil, N Kiani, and J Tegnér.
Methods and applications of algorithmic complexity, 2020.

[25]
Charles H Bennett.
Logical depth and physical complexity.
Citeseer, 1988.

[26]
Martin Mann, Heinz Ekker, and Christoph Flamm.
The graph grammar library-a generic framework for chemical graph rewrite systems.
In Theory and Practice of Model Transformations: 6th International Conference, ICMT 2013, Budapest, Hungary, June 18-19, 2013. Proceedings 6, pages 52–53. Springer, 2013.

[27]
Jakob L Andersen, Christoph Flamm, Daniel Merkle, and Peter F Stadler.
Inferring chemical reaction patterns using rule composition in graph grammars.
Journal of Systems Chemistry, 4:1–14, 2013.

[28]
Jakob L Andersen, Christoph Flamm, Daniel Merkle, and Peter F Stadler.
Defining autocatalysis in chemical reaction networks.

[29]
Abhishek Sharma, Dániel Czégel, Michael Lachmann, Christopher P Kempes, Sara I Walker, and Leroy Cronin.
Assembly theory explains and quantifies selection and evolution.
Nature, 622(7982):321–328, 2023.

[30]
Christopher Kempes, Sara I Walker, Michael Lachmann, and Leroy Cronin.
Assembly theory and its relationship with computational complexity.
arXiv preprint arXiv:2406.12176, 2024.

[31]
Robert M Alaniz, Bin Fu, Timothy Gomez, Elise Grizzell, Andrew Rodriguez, Robert Schweller, and Tim Wylie.
Reachability in restricted chemical reaction networks.
arXiv preprint arXiv:2211.12603, 2022.

[32]
Jakob L Andersen, Sissel Banke, Rolf Fagerberg, Christoph Flamm, Daniel Merkle, and Peter F Stadler.
On the realisability of chemical pathways.
In International Symposium on Bioinformatics Research and Applications, pages 409–419. Springer, 2023.

[33]
Stephen Jay Gould.
Wonderful life: the Burgess Shale and the nature of history.
WW Norton & Company, 1989.

[34]
Zachary D Blount, Richard E Lenski, and Jonathan B Losos.
Contingency and determinism in evolution: Replaying life’s tape.
Science, 362(6415):eaam5979, 2018.

[35]
David J Gross.
The role of symmetry in fundamental physics.
Proceedings of the National Academy of Sciences, 93(25):14256–14259, 1996.

[36]
Peter Walters.
An introduction to ergodic theory, volume 79.
Springer Science & Business Media, 2000.

[37]
Polly Y Yu and Gheorghe Craciun.
Mathematical analysis of chemical reaction systems.
Israel Journal of Chemistry, 58(6-7):733–741, 2018.

[38]
James R Norris.
Markov chains.
Number 2. Cambridge university press, 1998.

[39]
Massimiliano Esposito, Katja Lindenberg, and Christian Van den Broeck.
Extracting chemical energy by growing disorder: efficiency at maximum power.
Journal of Statistical Mechanics: Theory and Experiment, 2010(01):P01008, 2010.

[40]
Jeremy L England.
Dissipative adaptation in driven self-assembly.
Nature nanotechnology, 10(11):919–923, 2015.

[41]
Christopher Jarzynski.
Nonequilibrium equality for free energy differences.
Physical Review Letters, 78(14):2690, 1997.

[42]
Gavin E Crooks.
Entropy production fluctuation theorem and the nonequilibrium work relation for free energy differences.
Physical Review E, 60(3):2721, 1999.

[43]
Tim Schmiedl and Udo Seifert.
Stochastic thermodynamics of chemical reaction networks.
The Journal of chemical physics, 126(4), 2007.

[44]
Dilip K Kondepudi, Benjamin De Bari, and James A Dixon.
Dissipative structures, organisms and evolution.
Entropy, 22(11):1305, 2020.

[45]
Abhishek Deshpande and Manoj Gopalkrishnan.
Autocatalysis in reaction networks.
Bulletin of mathematical biology, 76(10):2570–2595, 2014.

[46]
Alex Blokhuis, David Lacoste, and Philippe Nghe.
Universal motifs and the diversity of autocatalytic systems.
Proceedings of the National Academy of Sciences, 117(41):25230–25236, 2020.

[47]
Zhen Peng, Alex M Plum, Praful Gagrani, and David A Baum.
An ecological framework for the analysis of prebiotic chemical reaction networks.
Journal of theoretical biology, 507:110451, 2020.

[48]
Alberto P Muñuzuri and Juan Pérez-Mercader.
Unified representation of life's basic properties by a 3-species stochastic cubic autocatalytic reaction-diffusion system of equations.
Physics of Life Reviews, 2022.

[49]
Alex M Plum and David A Baum.
Aces in spaces: Autocatalytic chemical ecosystems in spatial settings.
arXiv preprint arXiv:2212.14445, 2022.

[50]
Nicola Vassena and Peter F Stadler.
Unstable cores are the source of instability in chemical reaction networks.
Proceedings of the Royal Society A, 480(2285):20230694, 2024.

[51]
Armand Despons, Yannick De Decker, and David Lacoste.
Structural constraints limit the regime of optimal flux in autocatalytic reaction networks.
arXiv preprint arXiv:2306.02366, 2023.

[52]
Charles R. Johnson, Ronald L. Smith, and Michael J. Tsatsomeros.
Matrix Positivity.
Cambridge Tract s in Mathematics. Cambridge University Press, 2020.

[53]
Daniel J. Bates, Jonathan D. Hauenstein, Andrew J. Sommese, and Charles W. Wampler.
Bertini: Software for numerical algebraic geometry.
Available at bertini.nd.edu with permanent doi: dx.doi.org/10.7274/R0H41PB5.

[54]
Hong Qian and Elliot L Elson.
Single-molecule enzymology: stochastic michaelis–menten kinetics.
Biophysical chemistry, 101:565–576, 2002.

[55]
Jinsu Kim and German Enciso.
Absolutely robust controllers for chemical reaction networks.
Journal of the Royal Society Interface, 17(166):20200031, 2020.

[56]
Yoshiya J Matsubara, Sandeep Ameta, Shashi Thutupalli, Philippe Nghe, and Sandeep Krishna.
Robustness of compositional heredity to the growth and division dynamics of prebiotic compartments.
arXiv preprint arXiv:2211.03155, 2022.

[57]
Steven H Strogatz.
Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering.
CRC press, 2018.

[58]
Raj Kumar Pathria.
Statistical mechanics.
Elsevier, 2016.

[59]
David F Anderson, Gheorghe Craciun, Manoj Gopalkrishnan, and Carsten Wiuf.
Lyapunov functions, stationary distributions, and non-equilibrium potential for reaction networks.
Bulletin of mathematical biology, 77(9):1744–1767, 2015.

[60]
David A Baum, Zhen Peng, Emily Dolson, Eric Smith, Alex M Plum, and Praful Gagrani.
The ecology–evolution continuum and the origin of life.
Journal of the Royal Society Interface, 20(208):20230346, 2023.

[61]
Martin Feinberg.
Foundations of chemical reaction network theory.

[62]
Eric Smith and Supriya Krishnamurthy.
Flows, scaling, and the control of moment hierarchies for stochastic chemical reaction networks.
Physical Review E, 96(6):062102, 2017.

[63]
Alexandre Lazarescu, Tommaso Cossetto, Gianmaria Falasco, and Massimiliano Esposito.
Large deviations and dynamical phase transitions in stochastic chemical networks.
The Journal of Chemical Physics, 151(6):064117, 2019.

[64]
David F Anderson and Thomas G Kurtz.
Stochastic analysis of biochemical systems, volume 674.
Springer, 2015.

[65]
H Ge, H Qian, W Dubitzky, O Wolkenhauer, KH Cho, and H Yokota.
Chemical master equation.
Encyclopedia of systems biology, pages 396–399, 2013.

[66]
Thomas G Kurtz.
The relationship between stochastic and deterministic models for chemical reactions.
The Journal of Chemical Physics, 57(7):2976–2978, 1972.

[67]
Alex Kamenev.
Field theory of non-equilibrium systems.
Cambridge University Press, 2023.

[68]
Hugo Touchette.
Legendre-fenchel transforms in a nutshell.
URL http://www. maths. qmul. ac. uk/ ht/archive/lfth2. pdf, 2005.

[69]
Adam Shwartz and Alan Weiss.
Large deviations for performance analysis: queues, communication and computing, volume 5.
CRC Press, 1995.

[70]
Andrea Agazzi, Amir Dembo, and Jean-Pierre Eckmann.
Large deviations theory for markov jump models of chemical reaction networks.
The Annals of Applied Probability, 28(3):1821–1855, 2018.

[71]
Herbert Goldstein, Charles Poole, and John Safko.
Classical mechanics, 2002.

[72]
Daniel T Gillespie.
Exact stochastic simulation of coupled chemical reactions.
The journal of physical chemistry, 81(25):2340–2361, 1977.

[73]
Eric Smith.
Intrinsic and extrinsic thermodynamics for stochastic population processes with multi-level large-deviation structure.
Entropy, 22(10):1137, 2020.

[74]
Eric Smith, Harrison B Smith, and Jakob Lykke Andersen.
Rules, hypergraphs, and probabilities: the three-level analysis of chemical reaction systems and other stochastic stoichiometric population processes.
bioRxiv, pages 2023–12, 2023.

[75]
Lawrence C Evans.
Partial differential equations, volume 19.
American Mathematical Society, 2022.

[76]
James Bisgard.
Mountain passes and saddle points.
Siam Review, 57(2):275–292, 2015.

[77]
Mark I Dykman, Eugenia Mori, John Ross, and PM Hunt.
Large fluctuations and optimal paths in chemical kinetics.
The Journal of chemical physics, 100(8):5735–5750, 1994.

[78]
Ruben Zakine and Eric Vanden-Eijnden.
Minimum-action method for nonequilibrium phase transitions.
Physical Review X, 13(4):041044, 2023.

[79]
Fritz Horn and Roy Jackson.
General mass action kinetics.
Archive for rational mechanics and analysis, 47:81–116, 1972.

[80]
Manoj Gopalkrishnan.
On the lyapunov function for complex-balanced mass-action systems.
arXiv preprint arXiv:1312.3043, 2013.

[81]
Ėrnest Borisovich Vinberg.
A course in algebra.
Number 56. American Mathematical Soc., 2003.

[82]
Riccardo Rao and Massimiliano Esposito.
Nonequilibrium thermodynamics of chemical reaction networks: Wisdom from stochastic thermodynamics.
Physical Review X, 6(4):041064, 2016.

[83]
Kohei Yoshimura, Artemy Kolchinsky, Andreas Dechant, and Sosuke Ito.
Housekeeping and excess entropy production for general nonlinear dynamics.
Physical Review Research, 5(1):013017, 2023.

[84]
Matteo Polettini and Massimiliano Esposito.
Irreversible thermodynamics of open chemical networks. i. emergent cycles and broken conservation laws.
The Journal of chemical physics, 141(2), 2014.

[85]
Yuji Hirono, Takashi Okada, Hiroyasu Miyazaki, and Yoshimasa Hidaka.
Structural reduction of chemical reaction networks based on topology.
Physical Review Research, 3(4):043123, 2021.

[86]
Petter Brändén.
Unimodality, log-concavity, real-rootedness and beyond.
Handbook of enumerative combinatorics, 87:437, 2015.

We note that the

In particular, we find that a sufficient condition for multistability is the functional species must be a boundary species, an autocatalytic species

and explain that, for multistability, these species must be autocatalytic and boundary species and part of an autocatalytic set in Sec. <ref> and <ref>, respectively. Furthermore, in Sec. <ref>, we also remark on the role of entropy production in our formalism.

§.§ Boundary species

* Without chemostatting or some form of control by the environment, the system can only have one detaield-balancede attractor.

* Multistability requires non-detailed balanced, and hence there must be a species in the system that interacts with the environment. This is what we term as boundary species.

* Thus, the first functional role that a species can have is that of a boundary species that extracts chemical work from the environment to maintain the system at a NESS.

In particular, we prove that this indeed is the case for our class of models in Thm. <ref>, and the maximum number of attractors possible under any choice of kinetics is equal to the number of functional species in our system plus one.

* The efficiency depends on :

* Rate constants

* Boundary CRN - embedding of the internal CRN into the complete CRN - define processes and calculate efficiency

* In the absence of chemostatting, Different systems can be hooked up to provide control on each other. Thus efficiency, while not intrinsic to the dynamics, is intrinsic to the physics.

* Say complexity is the same as dissipation in our case. Different ways of chemostatting/setting boundary changes the attractors as well. COnjecture that efficiency will be monotonic in total complexity.

§.§ Autocatalytic species

We call such species as a functional species

* Model specific observations

* Functional species must transduce to make possible a NESS

* Functioal species must be autocatalytic to give rise to another fixed point in the internal CRN by Descarte's rule of signs.

* These can be generalized and are make sense that if there are infinite functional species then the system is open-ended. This is the ratcheting mechanism, when new species are created, new flows get opened up and new possibilities of interaction with environment and autocatalysis open up.

Mention Transience

other words, when the model is started at the origin, the CRN can progressively evo

First, we use
The conceptual contribution of this work is to clarify the role played by the non-equilibrium potential (NEP) in determining the stability of an attractor,

* Mathematical contribution: Proof that the number of attractors are increased by introducing functionality to internal species. That is making them boundary-internal species that interact with boundary-external species needed for energy and mass conservation as well as autocatalytic driving. In the NESS.

* Number of zeros of polynomials counted by Descarte's rule of signs.

* The lyapunov function must have phases for rate constants - different phases select how the reaction network changes and different parts of it become prominent.

* In the external world, there are many different instances of the same chemistry, but different parts of the network are activated.

* Phenomenology: phases are interpreted as these particular roots of polynomials with these deterministic properties. The stochastic properties become important for temporal predictions in biology or stochastic computation.

* Hamilton-Jacobi and physics have a deep connection. Using modern methods, together we can also make it the formulation of biochemistry.

Recipe. Introduce autocatalytic subnetwork in internal species. Due to mass conservation, some external species needs to be introduced. Then change rate parametrs so that the saddle point is closer to the less complex attractor. Then the system can evolve.

In Sec. <ref>, we construct simple models of abstract polymer chemistry demonstrating three types of evolution namely: progressive, open-ended, and historically-contingent evolution. In Sec. <ref>, through a series of definitions, we generalize these models into a mathematical framework and explain how one may ascertain whether or not any CRN is capable of a particular type of evolution. In Sec. <ref>, we explain the role of autocatalysis in our models and remark on how it can be used to generate more realistic models with desired features. In Sec. <ref> we calculate the dissipation and transduction efficiency of a model and demonstrate that dissipation or transduction efficiency is not always correlated with evolution by our criterion. We conclude with a discussion of the key features in our framework and its contributions to understanding the origins of life in  Sec. <ref>.

* We explained how to summarize the stochastic behavior of a mulstistable rule-generated Markov process such that it is amenable to an evolutionary analysis.

* attractor-transition graph obtained by finding instantons.

* Complexity grpah is obtained by getting an algorithmic measure (or relative ranking) of complexity on CRN.