  * November 2022

# Dynamical mean-field theory: from ecosystems to reaction networks

Eric De Giuli1 and Camille Scalliet2 1 Department of Physics, Toronto
Metropolitan University111formerly Ryerson University, M5B 2K3, Toronto,
Canada  2 Department of Applied Mathematics and Theoretical Physics,
University of Cambridge, Wilberforce Road, Cambridge, CB3 0WA, United Kingdom
edegiuli<EMAIL_ADDRESS>
###### Abstract

Both natural ecosystems and biochemical reaction networks involve populations
of heterogeneous agents whose cooperative and competitive interactions lead to
a rich dynamics of species’ abundances, albeit at vastly different scales. The
maintenance of diversity in large ecosystems is a longstanding puzzle, towards
which recent progress has been made by the derivation of dynamical mean-field
theories of random models. In particular, it has recently been shown that
these random models have a chaotic phase in which abundances display wild
fluctuations. When modest spatial structure is included, these fluctuations
are stabilized and diversity is maintained. If and how these phenomena have
parallels in biochemical reaction networks is currently unknown. Making this
connection is of interest since life requires cooperation among a large number
of molecular species. In this work, we find a reaction network whose large-
scale behavior recovers the random Lotka-Volterra model recently considered in
theoretical ecology. We clarify the assumptions necessary to derive its large-
scale description, and reveal the underlying assumptions made on the noise to
recover previous dynamical mean-field theories. Then, we show how local
detailed balance and the positivity of reaction rates, which are key physical
requirements of chemical reaction networks, provide obstructions towards the
construction of an associated dynamical mean-field theory of biochemical
reaction networks. Finally, we outline prospects and challenges for the
future.

††:

  * Keywords: Disordered systems, Statistical field theory, Theoretical ecology, Reaction networks.

  *

##  1 Introduction

The development of a theoretical framework to describe life as a physical
system, with relevant order parameters that distinguish animate from inanimate
matter, is a major goal of 21st century physics [1, 2]. First, it has long
been accepted that such a framework should deal with conditions far from
thermodynamic equilibrium [3]. The absence of a Boltzmann distribution forces
one to generally consider the full dynamics of out-of-equilibrium systems.
Ideally, the theoretical framework should be dynamic in nature. A second key
aspect of life is that it involves, at any scale, the coordinated action of a
large number $N$ of chemical species, which interact via a complex set of
chemical reactions. This opens the possibility of exploiting a large-species
$N\to\infty$ limit to obtain a tractable theory, with the prospect of
obtaining universal results. This presence of variability at the level of
species chemical interactions also adds motivation to consider simplified
models with random interactions.

Dynamical mean-field theory (DMFT) is essentially the unique framework that
combines both of the above essential aspects of the living state. It is a
theoretical approach to disordered statistical models that starts from a
mesoscopic description. The disorder is typically introduced via the couplings
between degrees of freedom. The method relies on the assumption that the
number of degrees of freedom in the system is large; formally the
thermodynamic limit of an infinite number of species is taken. This dynamic
mean-field theory describes the time evolution of a typical degree of freedom
after the average over the quenched disorder has been carried out, and yields
a closed description of the dynamics in terms of an effective species and a
small number of order parameters. First developed to describe the dynamics of
spin-glasses [4], DMFT has since been successful in identifying and describing
phase transitions in other areas of physics, such as neural networks [5, 6],
structural glasses [7, 8, 9], and ecosystems [10, 11, 12]. Dynamical mean-
field theory thus appears as a promising framework to describe life as a
physical system characterized by distinct dynamical phases.

Recently, substantial progress was made to obtain a dynamical mean-field
theory of ecosystems, using the generalized Lotka-Volterra model as a starting
point [10, 11, 12]. This theory has shed light on the ‘diversity paradox’,
namely the counterintuitive fact that many ecosystems are much more diverse
than what is expected from niche-filling. Importantly, the theory has revealed
that dynamical models for ecosystems can exist in a variety of phases,
characterized either by a single equilibrium, multiple marginally-stable
equilibria, or chaos [13, 14]. The framework has thus allowed extension of the
notion of phase away from thermodynamic equilibrium.

Compared to ecosystems, the organization of life occurs primarily at much
smaller scales, where molecules are the main actors. As stressed by Anderson
[15], life requires a diversity of chemical species, which are often involved
in complex chemical reactions networks (CRN). Even in the constant temperature
case common in biology, CRNs can exhibit rich dynamical behavior: unstable
steady states, multiple steady states, sustained composition oscillations, and
chaotic dynamics [16, 17, 18]. These behaviors have a parallel in model
ecosystems, where their existence has been firmly established by DMFT.

Our article builds naturally on past works employing a reaction network
representation of ecological models [19, 20, 21, 22, 23]. Our interest is
however not only to use CRN as a means to study model ecosystems, but to
investigate to what extent the recent advances made possible by the
development of DMFT in theoretical ecology can advance our understanding of
emergent behavior in CRNs. We specifically investigate the behavior of CRNs in
the limit of many chemical species which, to the best of our knowledge, has
not been considered in the past.

The overarching goal of this article is to ask: to what extent can one develop
a DMFT for general chemical reaction networks? Our strategy is to first find a
chemical reaction network whose large-scale behavior yields the generalized
Lotka-Volterra (GLV) model studied in ecology. We then examine which features
of such a GLV reaction network are universal, and which features of general
reaction networks are absent from the mapping. We will conclude that the
present DMFTs for ecology do not apply to generic chemical reaction networks,
suggesting the existence of universality classes of DMFTs yet to be explored.

The article is organized as follows. In section 2, we review the two existing
analytical techniques on which our work is based: dynamical mean field theory,
and the Doi-Peliti path integral formalism for chemical reaction networks. In
section 3 we find a microscopic reaction network which maps to the generalized
Lotka-Volterra model studied in the ecology literature. We show that the
derivation of the ecology model is not systematic and requires stringent, and
in general unjustified, assumptions. We then derive a generalized dynamical
mean-field theory starting from the path integral, and recover the large-
species DMFT obtained previously. In section 4, we explore the consequences of
the mapping between chemical reaction networks and theoretical ecology. In
section 5, we conclude by discussing future directions for the development of
new classes of DMFT.

##  2 Background

In this section, we review two existing analytical techniques that we will
combine in this article: (i) dynamical mean-field theory, and (ii) the Doi-
Peliti path integral formulation for chemical reaction networks.

###  2.1 A brief review of dynamical mean-field theory

Dynamical mean-field theory (DMFT) is a general technique which allows to
describe the dynamical behaviour of mean-field disordered systems with a
minimal number of relevant order parameters. It is distinct from the dynamical
mean-field theory of correlated electron systems [24], which describes the
behaviour of a quantum-mechanical system in terms of an effective impurity.
DMFT was first introduced to study the dynamics of mean-field spin glasses
[4], and later applied to the spherical $p$-spin model [25, 26]. The long-time
asymptotic solution to the equations for the latter exhibits the aging
phenomenon and explains it in terms of an effective temperature, a great
achievement for DMFT [26]. Since then, DMFT has successfully been applied to
neural networks [5, 6], structural glasses [7, 8, 9], and ecosystems [10, 11,
12]. To the best of our knowledge, DMFT has not previously been applied to
chemical reaction networks.

The main ingredients to construct a dynamical-mean field theory are the
following. Consider a generic theory with dynamical variables $\sigma_{i}(t)$
with $i=1,\ldots,N$ in the presence of disorder. Depending on the context,
these could be spin variables, or species abundances, for which disorder is
introduced via random coupling coefficients, or random inter-species
interactions. Dynamical mean-field theories generally begin with a Langevin
dynamics for the variables $\sigma_{i}(t)$. Either by path-integral methods
[27, 28, 4], or the dynamical cavity method [29], one performs a disorder
average to eliminate it. The many-species dynamical problem is then reduced to
an effective process for a representative degree of freedom. A consequence of
the disorder average is however to introduce two-time quantities. The most
basic of these are the correlation function
$C(t,t^{\prime})=\mbox{$\frac{1}{N}$}\sum_{i}\langle\sigma_{i}(t)\sigma_{i}(t^{\prime})\rangle$
and the response function
$R(t,t^{\prime})=\mbox{$\frac{1}{N}$}\sum_{i}\langle\left.\delta\sigma_{i}(t)/\delta
h_{i}(t^{\prime})\right|_{h=0}\rangle$, where $h_{i}(t)$ is a field conjugate
to the dynamical variable $\sigma_{i}(t)$. Here $\langle\cdot\rangle$ denotes
a thermal average subject to appropriate initial conditions. In the simplest
class of DMFTs, one obtains a pair of integro-differential equations involving
only $C(t,t^{\prime})$ and $R(t,t^{\prime})$ [4, 25]. One can also find
theories for which the equations for the evolution of $C$ and $R$ involve a
memory kernel that must be determined from a self-consistent 1D stochastic
process [10, 7, 11, 12]. More generally, we include in the DMFT category any
theory that reduces the number of degrees of freedom from $N\to\infty$ down to
a finite number. We will distinguish multiple levels of complexity in our
work.

###  2.2 A brief review of the Hamiltonian formulation of chemical reaction
networks

We consider a closed volume $\Omega$ containing a mixture of $N$ chemical
species denoted $X_{i}$ with $i=1,...,N$. We note $n_{i}(t)$ the number of
molecules $X_{i}$ at time $t$. We focus on the well-mixed limit and consider
homogeneous concentrations, expressed in units of $c^{\circ}$, the standard
concentration 1 mol/L. The state of the mixture is fully specified by the
vector $\vec{n}(t)$ of the number of molecules of each species, which changes
when a chemical reaction occurs. We consider chemical reactions, labelled by
$\alpha$, of the form

$\displaystyle\sum_{i}p_{\alpha
i}X_{i}\xrightharpoonup{k_{\alpha}}\sum_{i}q_{\alpha i}X_{i},$ (1)

where $k_{\alpha}$ is the rate of reaction $\alpha$. By convention, $p_{\alpha
i}$ and $q_{\alpha i}$ are the stoichiometric coefficients for the molecular
species $X_{i}$ as reactant and product, respectively. The complete set of
reactions is understood to contain each reaction in both forward and backward
directions. In full generality, reactions need not conserve the number of
molecules, thus allowing interactions with external reservoirs.

####  2.2.1 Chemical master equation

While the rate equations describe the time evolution of the average number of
molecules, they do not capture the fluctuations intrinsic to the stochastic
nature of chemical reactions. To obtain a complete and closed description of
the dynamics, one must track the evolution of the probability distribution
$P(\vec{n},t)$, as described in Ref. [30]. This yields the chemical master
equation

$\displaystyle\frac{\partial P(\vec{n},t)}{\partial
t}=\sum_{\alpha}\tilde{k}_{\alpha}\left(\prod_{i}\mathscr{E}_{i}^{p_{\alpha
i}-q_{\alpha i}}-1\right)\prod_{j}\left(\frac{n_{j}!}{(n_{j}-p_{\alpha
j})!}\right)P(\vec{n},t),$ (2)

where the products run over the $N$ species, and $p_{\alpha j}=0$ if $X_{j}$
does not intervene in reaction $\alpha$. The probability for a collision
involving $p_{\alpha j}$ molecules $X_{j}$ is proportional to
$n_{j}(n_{j}-1)...(n_{j}-p_{\alpha j}+1)/(c^{\circ}\Omega)^{p_{\alpha j}}$.
This gives rise to the factorial terms in equation (2), and the rescaled
reaction rate
$\tilde{k}_{\alpha}=k_{\alpha}(c^{\circ}\Omega)^{1-\sum_{i}p_{\alpha i}}$. The
effect of each reaction $\alpha$ is to transform $n_{j}$ into $n_{j}+q_{\alpha
j}-p_{\alpha j}$. This is conveniently captured by the step operator
$\mathscr{E}_{i}$ which acts on any arbitrary function $f(\vec{n})$ as

$\displaystyle\mathscr{E}_{i}f(\vec{n})=f(\vec{n}+\vec{1}^{i}),\quad\mathscr{E}^{-1}_{i}f(\vec{n})=f(\vec{n}-\vec{1}^{i}),$

with $1^{i}_{j}=\delta_{ij}$. We assume a well-stirred setting in which there
is no spatial dependence of concentration, although the derivation
straightforwardly extends to the more general case [31, 32, 33].

However, the chemical master equation (2) cannot be solved analytically for
generic chemical reaction networks.

####  2.2.2 Doi formalism

To make progress, we adopt the Hamiltonian formalism developed by Doi [31]
which consists in an exact rewriting of the chemical master equation (2). Doi
observed that the step operators in (2) and associated factorial terms
$n_{j}!/(n_{j}-p_{\alpha j})!$ have a natural representation in terms of
bosonic creation and annihilation operators. Reaction systems thus have a
natural interpretation in terms of quantum mechanics through second
quantization.

In this formalism, one introduces a non-Hermitian quasi-Hamiltonian, or
Liouvillian, $H$. Unlike the rate equations, the Hamiltonian $H$ completely
specifies the reaction network, as the chemical master equation does. Doi’s
construction and the corresponding path integral formulation are reviewed in
A. The resulting Hamiltonian $H$ is a sum of contributions from each reaction,
$H=\sum_{\alpha}H_{\alpha}$, with

$\displaystyle H_{\alpha}(\vec{\phi},\vec{\phi}^{*})$
$\displaystyle=\tilde{k}_{\alpha}\left[\prod_{j}\phi^{*}_{j}{}^{q_{\alpha
j}}-\prod_{j}\phi^{*}_{j}{}^{p_{\alpha j}}\right]\prod_{i}\phi_{i}^{p_{\alpha
i}}.$ (3)

The arguments $\vec{\phi},\vec{\phi}^{*}$ are abstract quantities that get
introduced first via creation and annihilation operators, and then as complex
variables in the path-integral construction. Only their product,
$\phi^{*}_{j}\phi_{j}=n_{j}$, has a direct physical interpretation as the
number of molecules of $X_{j}$. We thus perform a Hopf-Cole transformation
[34, 35, 36, 37, 38]

$\displaystyle\phi^{*}_{j}=e^{\nu_{j}},\quad\phi_{j}=n_{j}e^{-\nu_{j}},$ (4)

in terms of which we have

$\displaystyle H_{\alpha}(\vec{n},\vec{\nu})$
$\displaystyle=\left[e^{\vec{\nu}\cdot\vec{s}_{\alpha}}-1\right]\underbrace{\tilde{k}_{\alpha}\prod_{i}n_{i}^{p_{\alpha
i}}}_{\equiv F_{\alpha}(\vec{n})},$ (5)

where we introduced the stoichiometric matrix ${s}_{\alpha j}=q_{\alpha
j}-p_{\alpha j}$. Here we are abusing notation slightly to write $H$ as a
function of its new arguments $\vec{n}$ and $\vec{\nu}$, in place of the
previous ones. In what follows we always intend this form for $H$. The
doubling of degrees-of-freedom in Eqs.(3),(5) is characteristic of dynamic
problems [39, 40, 34].

The $\vec{\nu}$ variables do not have a straightforward interpretation. They
can be considered as a per-species bias or ‘tilt’ [38]. As explained below,
$\vec{\nu}=0$ corresponds to mean-field trajectories.

Several features of (5) are important. First, the $\vec{\nu}$-dependent factor
is independent of the kinetics: conservation of probability requires that
$H_{\alpha}(\vec{n},0)=0$ for any $\vec{n}$, and the specific form in (5)
follows from the Doi-Peliti construction for reaction networks with a
stoichiometric matrix ${s}_{\alpha j}$. Second, the kinetics are encoded in
the function $F_{\alpha}(\vec{n})$, which can be written as

$\displaystyle
F_{\alpha}(\vec{n})=k_{\alpha}(c^{\circ}\Omega)\prod_{i}\left(\frac{n_{i}}{c^{\circ}\Omega}\right)^{p_{\alpha
i}}.$ (6)

The specific form in (6) corresponds to the usual mass-action kinetics, but
other choices are possible, so long as $F_{\alpha}(\vec{n})\geq 0$. Finally,
the Hamiltonian (5) specifies only the dynamics of the system, and does not
contain any information about thermodynamics. We will return to the
consequences of thermodynamics, in particular the detailed balance condition,
below.

####  2.2.3 Peliti integral

Using Doi’s Hamiltonian formulation, Peliti proposed a path-integral
description of reaction systems, constructed as in quantum mechanics [32]. In
principle, one can obtain the complete statistics of species abundances by
evaluating the path integral

$\displaystyle Z=\int{\mathscr{D}}n\int{\mathscr{D}}\nu\;e^{-S},$

with the action

$\displaystyle
S=\int_{0}^{t_{f}}dt\left[\vec{\nu}\cdot\partial_{t}\vec{n}-H(\vec{n},\vec{\nu})\right]+S_{BC},$
(7)

where

$\displaystyle
S_{BC}=[\vec{n}(0)-\vec{n}^{0}]\cdot\vec{\nu}(0)+\sum_{j}n_{j}(t_{f})[-e^{-\nu_{j}(t_{f})}+1-\nu_{j}(t_{f})],$

encodes the boundary conditions at initial $t=0$, and final $t_{f}$ times.
Here we consider $n_{i}(0)=n^{0}_{i}$ as initial condition.

In practice, observables are computed by adding sources, or fields, to the
partition function $Z$. For that purpose, we consider the generating function

$\displaystyle
Z(\vec{z},t_{f})=\sum_{\\{\vec{n}\\}}z_{1}^{n_{1}}z_{2}^{n_{2}}\cdots
z_{N}^{n_{N}}P(\vec{n},t_{f}).$

The partition function is recovered for $\vec{z}=1$, and moments are computed
by derivatives at $\vec{z}=1$. Elgart and Kamenev [41] showed that values
$\vec{z}$ different from unity are necessary to understand rare fluctuations.
They also showed that the path-integral representation of the generating
function is that of the original partition function ($\vec{z}=1$), equation
(7), with differences only in the boundary terms

$\displaystyle
S_{BC}=\sum_{j}\left[[\vec{n}(0)-\vec{n}^{0}]\cdot\vec{\nu}(0)+n_{j}(t_{f})[-z_{j}e^{-\nu_{j}(t_{f})}+1-\nu_{j}(t_{f})]\right].$

This gives boundary conditions $\vec{\nu}(t_{f})=\log\vec{z}$,
$\vec{n}(0)=\vec{n}^{0}$. The particle number statistics are extracted from
the generating function by appropriate contour integrals. For example, the
marginal distribution of the $j^{th}$ species is

$\displaystyle\rho_{j}(m_{j})$ $\displaystyle=\oint\frac{dz_{j}}{2\pi
i}\frac{1}{z_{j}^{m_{j}+1}}Z(\vec{1}+(z_{j}-1)\hat{e}_{j},t)$

where $\vec{1}=(1,1,\ldots,1)$ and $\hat{e}_{j}=(0,\ldots,0,1,0,\ldots,0)$
with a 1 in the $j^{th}$ position, and the integral is around the unit circle.
Physical quantities are thus determined by $Z(\vec{z},t_{f})$ in the vicinity
of $\vec{z}=\vec{1}$.

To evaluate the path integral $Z$ in practice, one resorts to approximations.
Since $\vec{\nu}(t_{f})=\log\vec{z}$, and we need to consider $Z$ in the
vicinity of $\vec{z}=\vec{1}$, this motivates a small-$\nu$ expansion. From
conservation of probability, the action begins with $\mathscr{O}(\nu)$. To
leading order, $\mathscr{O}(\nu)$ in $S$, one obtains the classical rate
equations. To the next order, $\mathscr{O}(\nu^{2})$, one obtains the chemical
Langevin equation, derived below.

####  2.2.4 Instanton equations

Instead of making a small-$\nu$ expansion of the action, one can look for
saddle-points, i.e. paths that give the largest contribution to the integral.
To motivate this, we consider a naive scaling analysis of the action. Suppose
that the species numbers are large, $n\gg 1$. Since the number of reactants
$p_{\alpha i}$ is generally varying from reaction to reaction, we see from (5)
and (6) that in order to have nontrivial dynamics for general reaction
networks, we should have $c^{\circ}\Omega\gg 1$, with the dimensionless
concentrations $n_{i}/(c^{\circ}\Omega)=\mathscr{O}(1)$. In this case $H\sim
c^{\circ}\Omega\gg 1$ and a saddle-point analysis of the path integral is
formally justified. Following the nomenclature of field theory we call these
saddle-point paths ‘instantons’ [42]. The instanton equations, derived in B,
are

$\displaystyle\frac{\partial\vec{n}}{\partial t}$
$\displaystyle=+\frac{\partial
H}{\partial\vec{\nu}}=\sum_{\alpha}F_{\alpha}(\vec{n})\vec{s}_{\alpha}e^{\vec{\nu}\cdot\vec{s}_{\alpha}},$
(8a) $\displaystyle\frac{\partial\vec{\nu}}{\partial t}$
$\displaystyle=-\frac{\partial H}{\partial\vec{n}}=\sum_{\alpha}\frac{\partial
F_{\alpha}(\vec{n})}{\partial\vec{n}}\left[e^{\vec{\nu}\cdot\vec{s}_{\alpha}}-1\right].$
(8b)

We see immediately that $\vec{\nu}=0$, corresponding to mean-field
trajectories, is always a solution. In that case (8a) reduces to the reaction
rate equations.

Eqs.(8a),(8b) take the form of Hamilton’s equations, so several properties
familiar from classical mechanics also hold here. In particular, along an
instanton the time-evolution of the Hamiltonian is

$\displaystyle\frac{dH}{dt}=\frac{\partial
H}{\partial\vec{\nu}}\cdot\frac{\partial\vec{\nu}}{\partial t}+\frac{\partial
H}{\partial\vec{n}}\cdot\frac{\partial\vec{n}}{\partial t}+\frac{\partial
H}{\partial t}=\frac{\partial H}{\partial t}.$

We see that if the Hamiltonian is time-independent, it is conserved by the
dynamics along the instantons.

####  2.2.5 Expansion to obtain the Langevin equation

As explained in C, keeping terms to order $\mathscr{O}(\nu^{2})$ in the
action, one obtains the chemical Langevin equation

$\displaystyle\partial_{t}\vec{n}=\vec{h}(\vec{n},t)+B(\vec{n},t)^{1/2}\cdot\vec{\xi\;},$
(9)

where the drift $\vec{h}$ and correlation matrix $B$ are related to the
Hamiltonian

$\displaystyle\vec{h}(\vec{n},t)$
$\displaystyle=\partial_{\vec{\nu}}H(\vec{n}(t),0),$ (10) $\displaystyle
B(\vec{n},t)$
$\displaystyle=\partial_{\vec{\nu}}\partial_{\vec{\nu}}H(\vec{n}(t),0),$ (11)

and $\vec{\xi\;}$ is a white noise

$\displaystyle\langle\xi_{j}(t)\rangle$
$\displaystyle=0,\quad\langle\xi_{j}(t)\xi_{k}(t^{\prime})\rangle=\delta(t-t^{\prime})\delta_{jk}.$

The boundary conditions are $\vec{n}(0)=\vec{n}^{0}$ and trajectories are
weighted with $e^{-S}$, with

$\displaystyle
S=\mbox{$\frac{1}{2}$}\sum_{j}\frac{n_{j}(t_{f})(1-z_{j}^{2})}{z_{j}}+\mbox{$\frac{1}{2}$}\sum_{j}\log(z_{j}n_{j}(t_{f})).$

The late-time boundary contribution is specific to the generating function;
alternative boundary conditions are discussed in the C.

From the derivation of the Langevin equation we can infer that $z-1\sim
1/\sqrt{n}$, thus $\nu\sim\log z\sim 1/\sqrt{n}$. This justifies a posteriori
the small-$\nu$ expansion for large systems, and is also consistent with the
saddle-point analysis, since the action scales as $n/\sqrt{n}\gg 1$. More
precisely, for a species $j$ the small-$\nu_{j}$ expansion is justified when
$n_{j}$ is large. This small-$\nu$ expansion is thus the system-size expansion
in path-integral language.

####  2.2.6 Thermodynamics and detailed balance

To this point, we have only considered the kinetics of reaction networks. We
now turn to the constraints imposed by thermodynamics. Consistency with
thermodynamics requires that if a reaction $\alpha$ occurs with rate
$k^{+}_{\alpha}$, then its reverse must occur with rate $k^{-}_{\alpha}$. The
ratio of the rates is fixed by $(\Delta
G)_{\alpha}=(\vec{q}_{\alpha}-\vec{p}_{\alpha})\cdot\vec{G}$, the difference
in molar Gibbs free energy between products and reactants, expressed in units
of $RT$

$\displaystyle\frac{k^{+}_{\alpha}}{k^{-}_{\alpha}}=e^{-(\Delta G)_{\alpha}},$
(12)

which is known as local detailed balance [43]. Using modern transition state
theory [44], one can predict from first principles

$\displaystyle k^{+}_{\alpha}=k_{0}\;e^{-(\delta G)_{\alpha}},$ (13)

where we introduced $(\delta
G)_{\alpha}=G_{A_{\alpha}}-\vec{p}_{\alpha}\cdot\vec{G}$, the difference in
molar Gibbs free energy between the activated complex $A_{\alpha}$ of reaction
$\alpha$ and the reactants, and $k_{0}=1/(2\pi\beta\hbar)\approx 6\times
10^{12}$ Hz at 300K. Eq. (13) holds for both forward and backward reactions,
mutatis mutandis so that (12) holds

$\displaystyle\frac{k^{+}_{\alpha}}{k^{-}_{\alpha}}=e^{-(G_{A_{\alpha}}-\sum_{i}p_{\alpha
i}G_{i})}e^{+(G_{A_{\alpha}}-\sum_{i}q_{\alpha i}G_{i})_{\alpha}}=e^{-(\Delta
G)_{\alpha}}.$

This allows us to rewrite the kinetic part of the Hamiltonian as

$\displaystyle
F_{\alpha}(\vec{n})=ke^{-G_{A_{\alpha}}}\prod_{i}\left(\frac{n_{i}}{c^{eq}_{i}\Omega}\right)^{p_{\alpha
i}},$ (14)

with $k=k_{0}c^{\circ}\Omega$ and $c^{eq}_{i}=c^{\circ}e^{-G_{i}}$.

When a reaction occurs in a forward-reverse pair $\alpha_{+}$ and
$\alpha_{-}$, its contribution to the right-hand side of the $\partial
n/\partial t$ instanton equation is

$\displaystyle
F_{\alpha_{+}}\vec{s}_{\alpha_{+}}e^{\vec{\nu}\cdot\vec{s}_{\alpha_{+}}}+F_{\alpha_{-}}\vec{s}_{\alpha_{-}}e^{\vec{\nu}\cdot\vec{s}_{\alpha_{-}}}=\vec{s}_{\alpha_{+}}ke^{-G_{A_{\alpha}}}\left[e^{\vec{\nu}\cdot\vec{s}_{\alpha_{+}}}\prod_{i}\left(\frac{n_{i}}{c^{eq}_{i}\Omega}\right)^{p_{\alpha
i}}-e^{\vec{\nu}\cdot\vec{s}_{\alpha_{-}}}\prod_{i}\left(\frac{n_{i}}{c^{eq}_{i}\Omega}\right)^{q_{\alpha
i}}\right]$

The bracket vanishes if $n_{i}=c^{eq}_{i}\Omega$ and $\nu_{i}=0$ for all $i$.
Thus if all reactions occur in such pairs, this is a steady-state solution.
This is the well-known expression for thermal equilibrium in a CRN [30].

The fact that chemical reactions satisfy detailed balance has profound
physical and conceptual consequences. In particular, this makes the system
time-reversal symmetric, at least statistically. For any field theory,
symmetries of the system appear as symmetries in the action. There is thus a
symmetry operation for the action that reverses time, as detailed in D. When
all reactions occur in forward-reverse pairs satisfying local detailed
balance, the transformation

$\displaystyle\tilde{\nu}_{i}(t_{f}-t)$
$\displaystyle=-\nu_{i}(t)+\log(n_{i}(t)/(c^{eq}_{i}\Omega)),$
$\displaystyle\tilde{n}_{i}(t_{f}-t)$ $\displaystyle=n_{i}(t),$

is a symmetry of the action. Note that the transformation is nonlocal in time.
Steady solutions, that are invariant under time-reversal, satisfy
$2\nu_{i}=\log(n_{i}/(c^{eq}_{i}\Omega))$. The thermal equilibrium solution
$n_{i}/\Omega=c^{eq}_{i}$ is recovered when $\nu_{i}=0$, as expected [30].

####  2.2.7 Existence of Lyapunov function

We have shown that $\vec{\nu}=0$ along mean-field solutions, hence $H=0$ from
equation (5). Since $H$ is conserved by the instanton dynamics, we also have
$H=0$ on relaxations toward mean-field solutions. The locus of points defined
by $H=0$ thus plays a special role [41, 34, 38]. Let us implicitly define a
functional $\Psi(\vec{n})$ by considering the equation
$H[\vec{\mu},\vec{n}]=0$ where $\vec{\mu}=\partial\Psi/\partial\vec{n}$. This
is a time-independent Hamilton-Jacobi equation [38, 45]. Let us compute the
derivative of $\Psi$ along the deterministic dynamics $\vec{\nu}=0$ [45]

$\displaystyle\frac{d\Psi}{dt}=\left.\frac{\partial\Psi}{\partial\vec{n}}\cdot\frac{d\vec{n}}{dt}\right|_{\vec{\nu}=0}=\left.\vec{\mu}\cdot\frac{\partial
H}{\partial\vec{\nu}}\right|_{\vec{\nu}=0}=\vec{\mu}\cdot\sum_{\alpha}F_{\alpha}(\vec{n})\vec{s}_{\alpha},$

where we have used Hamilton’s equation (8a). Subtracting
$H[\vec{\mu},\vec{n}]=0$ (5), we find that $\Psi$ is non-increasing under the
deterministic flow

$\displaystyle\frac{d\Psi}{dt}$
$\displaystyle=\sum_{\alpha}F_{\alpha}(\vec{n})\left[\vec{\mu}\cdot\vec{s}_{\alpha}-e^{\vec{\mu}\cdot\vec{s}_{\alpha}}+1\right]\leq
0,$

since $e^{a}-1-a\geq 0$ for all real $a$, and $F_{\alpha}\geq 0$. We have
shown that $\Psi$ is a Lyapunov function for the deterministic dynamics. The
quantity $-\Psi$ has the interpretation of a non-equilibrium entropy
functional, and $-d\Psi/dt$ corresponds to the total entropy change, including
both entropy production and entropy exchanges [45].

The existence of a Lyapunov function for chemical reaction networks in the
deterministic limit has important consequences. In fact, complex chemical
reaction networks may have multiple attractors [46, 47, 48, 37, 38]. In the
limit of small systems, hopping between these attractors may occur, giving
rise to interesting dynamical behavior. The existence of a Lyapunov function
in the deterministic limit, recovered for large system size, however prevents
such hopping dynamics.

A Lyapunov function is also a practical tool, since an explicit expression for
it can be used to study the statistics of attractors. When noise is present,
then at the level of the Langevin dynamics, there may be no Lyapunov function,
but still there is an effective Hamiltonian governing the long-time behaviour
[49]. For the Lotka-Volterra model studied below, such an effective
Hamiltonian has been used [14].

##  3 Generalized Lotka-Volterra model as a reaction network

###  3.1 Generalized Lotka-Volterra model

Consider an isolated ecosystem with $N$ species $i=1,\ldots,N$ and abundances
$n_{i}(t)$ of each species. We assume that the ecosystem is well-mixed and do
not consider spatial dependences of the abundances. A minimal model of
population dynamics in the ecosystem must include a birth rate and death rate
for each species, along with inter-species interactions. The generalized
Lotka-Volterra (GLV) model includes these main features. It also includes a
small immigration rate $\lambda$ which prevents extinctions. The stochastic
nature of births and deaths, causing populations to fluctuate even without
external noise, is captured within the GLV by demographic noise. We are
eventually interested in large ecosystems, whose complexity and inherent
randomness is modeled by taking random inter-species interaction strengths
[50]. We consider the GLV model in a Langevin form

$\displaystyle\frac{dn_{i}}{dt}=\frac{r_{i}}{\kappa_{i}}(n_{i}\kappa_{i}-n_{i}^{2})-\sum_{j\neq
i}S_{ij}n_{i}n_{j}+\lambda+\sqrt{n_{i}}\xi_{i},$ (15)

where $r_{i}$ and $\kappa_{i}$ are respectively the birth rate and carrying
capacity of species $i$, and $\xi_{i}$ is a white noise with
$\langle\xi_{i}\rangle=0,\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=2\omega^{2}\delta_{ij}\delta(t-t^{\prime})$.
The interaction matrix elements $S_{ij}$ are Gaussian random variables, with
mean, variance and asymmetry given by

$\displaystyle\overline{S_{ij}}=s/N,\quad\overline{(S_{ij}-\overline{S_{ij}})^{2}}=\sigma^{2}/N,\quad\overline{(S_{ij}-\overline{S_{ij}})(S_{ji}-\overline{S_{ij}})}=\gamma\sigma^{2}/N,$
(16)

with a scaling of the cumulants ensuring a proper large ecosystem limit $N\gg
1$. The parameter $-1\leq\gamma\leq 1$ characterizes the degree of asymmetry
in the inter-species interactions.

The generalized Lotka-Volterra model Eq. (15) has been studied in the limit of
many interacting species $N\gg 1$, both in the deterministic limit where
$\omega=0$ [13, 51] for general $r_{i}$ and $\kappa_{i}$, and with finite
demographic noise [52] with parameters $r_{i}=1,\kappa_{i}=1$. The phase
diagram was also established with and without demographic noise in the limit
of symmetric strongly competitive interactions [53]. The species-local part of
the right-hand side ($n_{i}\kappa_{i}-n_{i}^{2}$) can be generalized to other
polynomials $f_{i}(n_{i})$, see [14, 52]. Recently, a DMFT has also been built
for the GLV considering structured, non-random, interactions [54].

Although one might expect demographic noise to be irrelevant in large
populations, it is crucial when the underlying system has multiple attractors,
as discussed below. Moreover in related models, its importance has been
emphasized in its effect on spatial patterning [55, 19], temporal oscillations
[20], ecology-evolution interactions [21], and on the full distribution of
abundances [56]. Quite generally, in nonlinear complex systems, one expects
that the limits $t\to\infty$ and volume $\Omega\to\infty$ (or carrying
capacity $\kappa\to\infty$) do not commute [47]: if $t\to\infty$ is taken at
finite $\Omega$, then the system can explore multiple attractors and will
sample states with an effective “thermodynamic” distribution, while if
$\Omega\to\infty$ is taken at finite $t$, then instead the system can be stuck
in one of multiple attractors.

In both deterministic ($\omega=0$) and stochastic ($\omega>0$) forms, the GLV
model is necessarily approximate because population abundances are discrete
numbers, while the $n_{i}(t)$ in Eq. (15) are continuous variables. While this
difference is not expected to be critical when $n_{i}(t)\gg 1$, theory and
simulations of Eq. (15) show that abundances can vary over a wide range, with
many becoming small [13, 51, 14, 12]. This is particularly so in a phase with
multiple attractors, for which the dynamics is chaotic [14, 12]. In these
scenarios it is important to generalize Eq. (15) to a model that accounts for
the fundamental discreteness of species abundances.

Our goal here is to find a reaction network whose macroscopic behaviour yields
Eq. (15) in some limit. By doing so, the tools of chemical reaction network
theory reviewed in section 2 can be brought to bear on ecosystems. In the few-
species case, the GLV model has a standard interpretation as a reaction
system, for example between predator and prey, and reaction network
representations have been used elsewhere in ecology [19, 20, 21, 22, 23].
Likewise the Doi-Peliti path-integral formulation for CRNs has been adopted
for ecological models (reviewed in [57]), but not in the many-species limit of
interest here. To our knowledge there is no exposition of the conditions under
which the GLV model Eq. (15), for which DMFT has been built, can be obtained
from a microscopic model. The reaction network interpretation is important
because it holds for an arbitrary number of individuals, and can therefore be
used to disambiguate possible continuum limits. Moreover this mapping will
allow us to discuss the generality of the ecology DMFT, and whether it has an
analogue for true chemical reaction networks.

###  3.2 Deterministic limit

We search for a reaction network whose dynamics coincides with that of the
ecology GLV model (15). As a first step, we consider the deterministic limit
and leave the treatment of the noise term to the next section. In order to
find the relevant set of reactions, we compare the deterministic terms in the
right-hand side of (15) with the general reaction equation (8a), evaluated in
the deterministic limit $\vec{\nu}=0$. We compare the chemical reaction (1)
with the corresponding reaction rate equations (8a) and assume, without loss
of generality, the form from transition state theory (14), with $k_{0}$
obviously arbitrary. Indeed, we are free to consider a reaction and not its
detailed-balanced reverse one, if necessary. We see that the power of $n_{j}$
in the reaction rate equation for species $j$ and reaction $\alpha$ indicates
the number of times it appears as a reactant in the reaction.

Given that the deterministic part of the GLV model (15) contains a positive
linear and a negative quadratic term in $n_{j}$, we consider the reaction pair

$\displaystyle A_{j}\xrightleftharpoons[k^{-}_{j}]{k^{+}_{j}}2A_{j}.$ (17)

These reactions have a clear ecological interpretation as asexual reproduction
in the $k^{+}$ direction, and competitive death in the $k^{-}$ direction. By
identifying the kinetic factor (14) with the GLV model (15), we see that the
pair of reactions (17) appears if we take
$ke^{-G_{j}}/(c^{eq}_{j}\Omega)=r_{j}$ and
$ke^{-G_{j}}/(c^{eq}_{j}\Omega)^{2}=r_{j}/\kappa_{j}$. This yields

$\displaystyle\kappa_{j}=c^{eq}_{j}\Omega,\qquad k^{+}_{j}=r_{j},\qquad
k^{-}_{j}=(r_{j}c^{\circ}\Omega)/\kappa_{j}.$

Therefore $\kappa_{j}$ is the carrying capacity of species $j$ in the absence
of interactions, $r_{j}$ is the birth rate, and
$(r_{j}c^{\circ}\Omega)/\kappa_{j}$ is the rate of intra-species competition.
These reactions are in detailed balance.

We now consider the deterministic part which models the inter-species
interactions, $\sum_{j}S_{ij}n_{i}n_{j}$. Contrary to the sum of linear and
quadratic terms, which can be coupled into a pair of reactions satisfying
detailed balance, the terms $S_{ij}n_{i}n_{j}$ have no reason to be paired.
The reactions which give rise to the inter-species interaction term will not
satisfy detailed balance, in general. In fact, consider the reactions

$\displaystyle A_{j}+A_{l}$
$\displaystyle\overset{\;v_{jl}}{\rightarrow}2A_{j}+A_{l}$ (18) $\displaystyle
A_{j}+A_{l}$ $\displaystyle\overset{\;w_{jl}\;}{\rightarrow}0+A_{l},$ (19)

which give a contribution to the reaction rate equation for $j$ of the form
$(v_{jl}-w_{jl})n_{j}n_{l}$. This corresponds to a term
$ke^{-G_{jl}}n_{j}n_{l}/(c^{eq}_{j}c^{eq}_{l}\Omega^{2})$ in (14) where
therefore $ke^{-G_{jl}}/(c^{eq}_{j}c^{eq}_{l}\Omega^{2})=v_{jl}-w_{jl}$. The
inter-species interaction terms in (15) are recovered if
$S_{jl}=w_{jl}-v_{jl}$. The rates $v_{jl}$ and $w_{jl}$ must each be
nonnegative. If $S_{jl}>0$, then (19) is needed, but (18) is not, while if
$S_{jl}<0$, then (18) is needed, but (19) is not. Since the GLV model studied
in [52, 14, 12] takes $S_{jl}$ to be random numbers with both signs, then in
general both reactions are needed, and we can take

$\displaystyle v_{jl}$ $\displaystyle=|S_{jl}|\Theta(-S_{jl})$ $\displaystyle
w_{jl}$ $\displaystyle=|S_{jl}|\Theta(+S_{jl}),$

since then $v_{jl}-w_{jl}=|S_{jl}|[\Theta(-S_{jl})-\Theta(+S_{jl})]=-S_{jl}$
as required. Here $\Theta(x)$ is the Heaviside function. In the ecology
context, the reactions (18) and (19) are interpreted as cooperative growth and
competitive death, respectively.

Finally, a spontaneous birth reaction is needed to recover the immigration
term in (15)

$\displaystyle 0$ $\displaystyle\overset{\lambda_{j}}{\rightarrow}A_{j}.$ (20)

We have thus found a reaction network, defined by reactions
(17),(18),(19),(20), whose rate equations recover exactly the deterministic
part of the GLV model (15). This completes the first part of the mapping.

The Hamiltonian for this GLV reaction network is given by (5)

$\displaystyle
H=\sum_{j}r_{j}[e^{\nu_{j}}-1]\left[n_{j}-\frac{n_{j}^{2}e^{-\nu_{j}}}{\kappa_{j}}\right]+\sum_{i,(j\neq
i)}S_{ij}[e^{-\nu_{i}}-1]n_{i}n_{j}+\sum_{i,(j\neq
i)}v_{ij}[e^{\frac{1}{2}\nu_{i}}-e^{-\frac{1}{2}\nu_{i}}]^{2}n_{i}n_{j}+\sum_{j}\lambda_{j}[e^{\nu_{j}}-1].$
(21)

###  3.3 Langevin regime

We now consider the finite-$\omega$ regime where demographic noise plays a
role. Using the reaction network (17),(18),(19), we investigate under which
conditions the noise term in (15) is recovered. Building on the results of
section 2.2.5, the Langevin equation for the reaction system is obtained by a
small-$\nu$ expansion of the Hamiltonian (21)

$\displaystyle H$
$\displaystyle=\sum_{j}r_{j}[\nu_{j}+\mbox{$\frac{1}{2}$}\nu_{j}^{2}]\left[n_{j}-\frac{n_{j}^{2}(1-\nu_{j})}{\kappa_{j}}\right]+\sum_{i,(j\neq
i)}S_{ij}[-\nu_{i}+\mbox{$\frac{1}{2}$}\nu_{i}^{2}]n_{i}n_{j}+\sum_{i,(j\neq
i)}v_{ij}\nu_{i}^{2}n_{i}n_{j}$
$\displaystyle\qquad+\sum_{j}\lambda_{j}[\nu_{j}+\mbox{$\frac{1}{2}$}\nu_{j}^{2}]+\mathscr{O}(\nu^{3}).$

We obtain the Langevin equation for the abundance of species $j$ using
equations (9), (10), (11)

$\displaystyle\partial_{t}n_{j}=r_{j}\left[n_{j}-\frac{n_{j}^{2}}{\kappa_{j}}\right]-\sum_{i\neq
j}S_{ji}n_{i}n_{j}+\lambda_{j}+\left[r_{j}(n_{j}+n_{j}^{2}/\kappa_{j})+n_{j}\sum_{k\neq
j}[w_{jk}+v_{jk}]n_{k}+\lambda_{j}\right]^{1/2}\xi_{j},$ (22)

which resembles that of the GLV model (15), except that the noise has richer
correlations.

In order to connect with recent works on the GLV model, we consider that all
species have the same birth rate and carrying capacity, $r_{j}=r$ and
$\kappa_{j}=\kappa$. We perform the scaling transformation

$\displaystyle n_{j}=\kappa n_{j}^{\prime},\quad t=t^{\prime}/r,\quad
w_{ij}=rw_{ij}^{\prime}/\kappa,\quad
v_{ij}=rv_{ij}^{\prime}/\kappa,\quad\lambda_{j}=r\kappa\lambda_{j}^{\prime},\quad\xi_{j}=r^{1/2}\xi^{\prime}_{j},$

under which the Langevin equation (22) becomes

$\displaystyle\partial_{t^{\prime}}n_{j}^{\prime}=n_{j}^{\prime}(1-n_{j}^{\prime})-\sum_{i\neq
j}S^{\prime}_{ji}n_{i}^{\prime}n_{j}^{\prime}+\lambda_{j}^{\prime}+\frac{1}{\sqrt{\kappa}}\left[n_{j}^{\prime}+n_{j}^{\prime
2}+n_{j}^{\prime}\sum_{k\neq
j}[w^{\prime}_{jk}+v^{\prime}_{jk}]n_{k}^{\prime}+\lambda_{j}^{\prime}\right]^{1/2}\xi^{\prime}_{j},$
(23)

with
$\langle\xi^{\prime}_{j}(t^{\prime})\xi_{k}^{\prime}(0)\rangle=\delta_{jk}\delta(t^{\prime})$.

The GLV model Eq. (15) considered in [52, 53] is recovered assuming that the
bracket in (23) can be approximated by $\alpha n_{j}^{\prime}$, the noise
strength then being given by

$\displaystyle 2\omega^{2}=\alpha/\kappa.$ (24)

To the best of our knowledge, this is the first derivation of the GLV model
starting from a microscopic description. In the next section, we critically
expose the assumptions underlying such an approximation, and detail their
consequences.

Before doing so, let us briefly discuss the difference between (23) and the
reduced model (15) with simplified noise. Consider (23) for a generic pair of
species $j$ and $k$. Suppose the abundance of $k$ is abnormally large:
$n_{k}^{\prime}\gg 1$. This will affect the dynamics of species $j$ in two
ways: first, and most directly, it will increase the rate of reactions in
which these species are reactants, i.e. increase the magnitude of the
$S^{\prime}_{jk}n_{j}^{\prime}n_{k}^{\prime}$ term. But secondly, and more
subtlely, it will increase the prefactor of the noise $\xi^{\prime}_{j}$. Thus
species $j$ becomes more susceptible to noise as a result of the fluctuation
in species $k$. This susceptibility is specific, because it is multiplied by
the corresponding rates $w^{\prime}_{jk}$ and $v^{\prime}_{jk}$. In principle,
it allows for fine-tuning of states to these noise-mediated interactions. When
the noise amplitude is approximated by $\alpha n_{j}^{\prime}$, this
susceptibility is lost.

###  3.4 Regime of validity of the generalized Lotka-Volterra model

The derivation of the GLV model Eq. (15) is not systematic because we need to
justify (i) the approximation of the noise, and (ii) the small-$\nu$
expansion. We expose the conditions under which (i) and (ii) are justified,
and therefore clarify the regime of validity of the GLV model studied recently
[52, 53].

To justify (i), consider first the case of infinitesimal immigration,
$\lambda=0^{+}$. The regime of interest is large populations, $\kappa\gg 1$.
Note that the equation of motion is scaled such that $n_{j}^{\prime}\sim 1$,
typically. For such species the noise is important when
$\xi^{\prime}_{j}/\sqrt{\kappa}\sim 1$. Since the noise is scaled to be of
unit magnitude, if $\kappa\gg 1$, then $\xi^{\prime}_{j}/\sqrt{\kappa}$ is
small, and noise is typically irrelevant. However, if the abundance of one
species gets very small $n_{j}^{\prime}\ll 1$, one can have
$n_{j}^{\prime}\sim\sqrt{n_{j}^{\prime}/\kappa}$ even for large $\kappa$, i.e.
when $n_{j}^{\prime}\sim 1/\kappa$. Thus for $\kappa\gg 1$ the noise term can
be approximated by its value for $n_{j}^{\prime}\ll 1$. In this case the
bracket in (23) is approximated by

$\displaystyle n_{j}^{\prime}\left[1+\sum_{k\neq
j}[w^{\prime}_{jk}+v^{\prime}_{jk}]n_{k}^{\prime}\right].$

The rates $w^{\prime}$ and $v^{\prime}$ are nonnegative, so in general the
second term does not cancel and must be retained. However, if we consider the
large species limit $N\gg 1$ where all $w^{\prime}_{jk}$ and $v^{\prime}_{jk}$
are identically and independently distributed, say with mean $s/N$, then the
leading term as $N\to\infty$ will be

$\displaystyle n_{j}^{\prime}\left[1+2s\overline{n}(t)+\ldots\right],$

where $\overline{n}(t)=\mbox{$\frac{1}{N}$}\sum_{k}n_{k}^{\prime}(t)$ is the
mean rescaled abundance. For the fluctuations to be irrelevant we must assume
that the $n_{k}^{\prime}$ are not strongly correlated with
$w^{\prime}_{jk}-s/N$ and $v^{\prime}_{jk}-s/N$. Assuming all this, we obtain
a noise strength of the form (24) except that $\alpha=1+2s\overline{n}(t)$ is
time-dependent. Now, if the system fluctuates around a non-equilibrium steady
state, or if the total population is explicitly fixed as in [11], then this
time-dependence disappears and (i) is justified. The argument holds for finite
immigration if it is subdominant in the noise, i.e.
$1/\kappa\gg\lambda^{\prime}$, which is equivalent to $\lambda\ll r$.

Let us try to justify (ii). As noted above, the derivation of the Langevin
equation implies that $\nu_{j}\sim 1/\sqrt{n_{j}}$. A necessary condition for
$\nu_{j}$ to take only infinitesimal values is $\kappa\to\infty$. This cannot
be sufficient, because when $n_{j}^{\prime}\sim 1/\kappa$ then clearly
$\nu_{j}$ is not small. Thus the small-$\nu$ expansion, while justified for
species of typical (large) abundance, has no clear justification if any
species tends to extinction.

To summarize, for the chemical reaction network composed of reactions (17),
(18), (19), (20), if we make the following assumptions (A)

(A1): $\displaystyle k^{+}_{j}=r,\quad
k^{-}_{j}=(rc^{\circ}\Omega)/\kappa,\quad\lambda_{j}=\lambda$ (A2):
$\displaystyle\kappa\to\infty,$

then we obtain the Langevin equation (23) in rescaled variables
$n_{j}^{\prime}=n_{j}/\kappa,t^{\prime}=t\;r$. The derivation is fully self-
consistent only if all rescaled abundances are $\mathscr{O}(1)$, that is,
there are no rare species. In this case, the noise is typically irrelevant.

The noise becomes important if any abundance becomes small,
$n_{j}^{\prime}\sim 1/\kappa$. In general the noise cannot be simplified
beyond (23), except omitting the $n_{j}^{\prime 2}$ term. It is then necessary
to make these additional assumptions

(A3): steady state (A4): $\displaystyle N\to\infty$ (A5):
$\displaystyle\lambda\ll r,$

in order to justify the noise amplitude approximation. If the GLV model with
demographic noise studied in the literature [52, 53] is intended to represent
the macroscopic limit of the microscopic process given by the above
‘reactions’, then the implicit assumptions which we expose in (A1-5) are
necessary for its validity. Of course, it may be possible that the GLV model
is justified by a different argument. For example, one could apply the
renormalization group to see which terms are nontrivially irrelevant in the
macroscopic limit; to our knowledge such an analysis has not been performed.

From the derivation we learn that the noise strength $\omega$ satisfies

$\displaystyle\omega^{2}=\frac{1+2s\overline{n}}{2\kappa},$

where $\overline{n}$ is the mean dimensionless abundance, $s/N$ is the mean of
rescaled reaction rates, and $\kappa$ is the carrying capacity of the species,
assumed all to be equal.

The above derivation highlights several features of the GLV model. First, its
main limitation is that the small-$\nu$ expansion is not justified when any
species has a small population. Since demographic noise is only important when
a population is small, any phenomenology present in the stochastic model which
is absent in the deterministic model must be carefully scrutinized to check
that it is insensitive to the correlations in the noise, which are customarily
neglected. This is more than an academic point since if noise is neglected
entirely, then there is a Lyapunov function, as discussed in 2.2.7. With a
Lyapunov function there can be no hopping between basins, so the dynamical
relevance of the multiple-attractors phase discussed in [14, 12, 52, 53]
relies upon the existence of noise.

A second feature is also highlighted, namely that when the coupling $S_{ij}$
takes both signs, it is a sum of non-negative and non-positive quantities,
representing different physical processes. We will get back to this point in
the discussion.

###  3.5 Universality of the deterministic limit

While the above derivation highlights subtleties in the treatment of
demographic noise, the deterministic limit, obtained by taking
$\kappa\to\infty$ and assuming that all rescaled abundances are
$\mathscr{O}(1)$, is well-motivated. We now show that it is also universal:
additional reactions can be present microscopically, without changing the GLV
model.

In particular, the reactions (18) and (19) are the simplest that lead to a
general asymmetric interaction $S_{jl}$, but they are not unique. For example,
we could add the symmetric reactions

$\displaystyle A_{j}+A_{l}$
$\displaystyle\overset{t_{jl}}{\rightarrow}2A_{j}+2A_{l}$ (25) $\displaystyle
A_{j}+A_{l}$ $\displaystyle\overset{\;u_{jl}\;}{\rightarrow}0$ (26)

which would also end up subsumed into $S_{jl}$. They would be distinguished at
the level of the full fluctuation matrix $B_{ij}$, but under the assumptions
(A1-5), and truncating the small$-\nu$ expansion at $\mathscr{O}(\nu^{2})$, we
will obtain again (15). This is a form of universality.

Eqs. (25) and (26) have an ecological interpretation in terms of cooperative
reproduction, and cooperative death, respectively. It is possible to add even
more reactions that maintain the universal deterministic limit, but the
ecological interpretation would be dubious. For example,
$A_{j}+A_{l}\rightarrow m(A_{j}+A_{l})$ also gets subsumed into the same
interaction term, but does not have an obvious ecological interpretation for
$m>2$.

Insensitivity to microscopic reactions is one element of universality; others
have been discussed in [13, 51].

###  3.6 Dynamical mean-field theory from the Doi-Peliti integral

The DMFT for (15) was originally derived with the dynamical cavity technique
[12]. Complementing this, we show here that the same DMFT can be derived from
the Doi-Peliti integral. Unlike the cavity method, which works directly in the
large number of species $N\to\infty$ limit, our technique is more general, and
highlights the essential features necessary to obtain a DMFT.

The key advantage of treating disordered systems with dynamics is that the
path integral, which is equal to one, can be directly averaged over the
disorder, without replicas [39, 27]. The path integral has an action $S$ given
by (7), with the reaction network Hamiltonian (21). We will initially set
$v_{ij}=0$ in (21), so that disorder only appears in the second term of the
Hamiltonian via the random variables $S_{ij}$. We thus need to compute the
disorder average (noted $\overline{\;\cdot\;}$) of

$\displaystyle e^{\int_{t}\sum_{i<j}[S_{ij}A_{ij}+S_{ji}A_{ji}]}$

where $A_{ij}(t)=[e^{-\nu_{i}(t)}-1]n_{i}(t)n_{j}(t)$. Following previous
works, we assume that the $\\{S_{ij}\\}$ are Gaussian, with mean and variance
given in (16). Using standard methods, outlined in E, the disorder-averaged
generating function is exactly transformed into a functional of several order
parameters, each a sum over contributions from each species. The order
parameters are:

$\displaystyle P(t,t^{\prime})$
$\displaystyle=\frac{1}{N}\sum_{i}e^{-\nu_{i}(t)}e^{-\nu_{i}(t^{\prime})}n_{i}(t)n_{i}(t^{\prime})$
$\displaystyle Q(t,t^{\prime})$
$\displaystyle=\frac{1}{N}\sum_{i}e^{-\nu_{i}(t)}n_{i}(t)n_{i}(t^{\prime})$
$\displaystyle C(t,t^{\prime})$
$\displaystyle=\frac{1}{N}\sum_{i}n_{i}(t)n_{i}(t^{\prime})$ $\displaystyle
p(t)$ $\displaystyle=\frac{1}{N}\sum_{i}e^{-\nu_{i}(t)}n_{i}(t)$
$\displaystyle m(t)$ $\displaystyle=\frac{1}{N}\sum_{i}n_{i}(t),$

and each is accompanied by a Lagrange multiplier, labelled
$\Lambda_{P}(t,t^{\prime}),\Lambda_{Q}(t,t^{\prime}),\Lambda_{C}(t,t^{\prime}),\lambda_{p}(t),$
and $\lambda_{m}(t)$, respectively. Although $C$ and $m$ have a clear physical
meaning in terms of time-correlation and average of the species abundancies,
the other order parameters do not have straightforward interpretations. Any
order parameter that involves $\nu$ measures the sensitivity of abundances to
fluctuations, and is therefore a type of response function. The standard,
simple response function will be obtained only in the limit of small $\nu$, as
explained below.

The final disorder-averaged generating function takes the form

$\displaystyle\overline{Z}$
$\displaystyle=\int{\mathscr{D}}[\cdots]e^{-NF}\prod_{j}\left(\int{\mathscr{D}}n_{j}\int{\mathscr{D}}\nu_{j}e^{-S_{j}}\right),$
(27)

where

$\displaystyle\int{\mathscr{D}}[\cdots]=\int{\mathscr{D}}\Lambda_{P}\int{\mathscr{D}}\Lambda_{Q}\int{\mathscr{D}}\Lambda_{C}\int{\mathscr{D}}\lambda_{p}\int{\mathscr{D}}\lambda_{m}\int{\mathscr{D}}P\int{\mathscr{D}}Q\int{\mathscr{D}}C\int{\mathscr{D}}p\int{\mathscr{D}}m,$

and

$\displaystyle F$
$\displaystyle=-\frac{\sigma^{2}}{2}\int_{t,t^{\prime}}\left[P(t,t^{\prime})C(t,t^{\prime})+\gamma
Q(t,t^{\prime})Q(t^{\prime},t)+(1+\gamma)C(t,t^{\prime})[-Q(t,t^{\prime})-Q(t^{\prime},t)+C(t,t^{\prime})]\right]$
$\displaystyle\quad+\int_{t,t^{\prime}}[\Lambda_{P}(t,t^{\prime})P(t,t^{\prime})+\Lambda_{Q}(t,t^{\prime})Q(t,t^{\prime})+\Lambda_{C}(t,t^{\prime})C(t,t^{\prime})]$
$\displaystyle\quad+\int_{t}[\lambda_{p}(t)p(t)+\lambda_{m}(t)m(t)-s(p(t)-m(t))m(t)],$

and

$\displaystyle S_{j}$
$\displaystyle=-\int_{t}V_{j}+S_{BC,j}+\int_{t}\nu_{j}\partial_{t}n_{j}-\int_{t,t^{\prime}}[\Lambda_{P}(t,t^{\prime})e^{-\nu_{j}(t)}e^{-\nu_{j}(t^{\prime})}+\Lambda_{Q}(t,t^{\prime})e^{-\nu_{j}(t)}+\Lambda_{C}(t,t^{\prime})]n_{j}(t)n_{j}(t^{\prime})$
$\displaystyle-\int_{t}[\lambda_{p}(t)e^{-\nu_{j}(t)}+\lambda_{m}(t)]n_{j}(t)+\frac{s}{N}\int_{t}A_{jj}(t)+\frac{\sigma^{2}}{2N}(1+\gamma)\int_{t,t^{\prime}}A_{jj}(t)A_{jj}(t^{\prime}),$

and where $V_{j}$ is the non-interacting part of the Hamiltonian (21)

$\displaystyle
V_{j}(n,\nu)=r_{j}[e^{\nu}-1]n\left[1-n\frac{e^{-\nu}}{\kappa_{j}}\right]+\lambda_{j}[e^{\nu}-1].$

For simplicity choose $S_{BC}$ so that it does not depend on $j$. Then the
single-species dependence in $S_{j}$ comes from $r,\kappa,$ and $\lambda$. Let

$\displaystyle
W(r,\kappa,\lambda)=\int{\mathscr{D}}n\int{\mathscr{D}}\nu\;e^{-S(r,\kappa,\lambda)}$
(28)

be the single-species generating function, suppressing the dependence on all
the Lagrange multipliers. Then we have

$\displaystyle\prod_{j}W(r_{j},\kappa_{j},\lambda_{j})$
$\displaystyle=\exp\left(\int dr\int d\kappa\int
d\lambda\;\sum_{j}\delta(r-r_{j})\delta(\kappa-\kappa_{j})\delta(\lambda-\lambda_{j})\log
W(r,\kappa,\lambda)\right)$ $\displaystyle\equiv\exp\left(N\int dr\int
d\kappa\int d\lambda\;\tilde{\rho}(r,\kappa,\lambda)\log
W(r,\kappa,\lambda)\right),$

defining the empirical density $\tilde{\rho}$ of the single-species
parameters. If these parameters are iid, then as $N\to\infty$, $\tilde{\rho}$
converges to the disorder-average density $\rho$ (more precisely the
cumulative distributions converge). (At finite $N$ we can consider the simpler
model where all $r_{j},\kappa_{j}$, and $\lambda_{j}$ are independent of $j$.)
Then we can write

$\displaystyle\prod_{j}W(r_{j},\kappa_{j},\lambda_{j})$
$\displaystyle\to\exp\left(N\;\overline{\log W}\right)$
$\displaystyle=\exp\left(N\;\overline{\log\int{\mathscr{D}}n\int{\mathscr{D}}\nu\;e^{-S}}\right)$
(29)

so that

$\displaystyle\overline{Z}$
$\displaystyle=\int{\mathscr{D}}[\cdots]\exp\left(-NF+N\;\overline{\log\int{\mathscr{D}}n\int{\mathscr{D}}\nu\;e^{-S}}\right)$

This is interesting because we have already reduced the path integral from
$2N$ dynamical variables to 12 dynamical variables. We did not need to assume
large number of species limit $N\to\infty$, and we did not need to make any
small-$\nu$ expansion. Thus this path integral still includes the effect of
rare fluctuations and demographic noise to all orders in $\nu$, and could be
useful to understand rare events [41]. In principle, it may be subject to
other analytical approximation methods that do not require large $N$.

In practice, the DMFT derivation proceeds by considering the limit
$N\to\infty$, in which case the action is proportional to $N$ and the saddle-
point method can be applied. The correlation and response functions are then
determined by saddle-point equations

$\displaystyle P(t,t^{\prime})$ $\displaystyle=\langle
e^{-\nu(t)}e^{-\nu(t^{\prime})}n(t)n(t^{\prime})\rangle$ (30a) $\displaystyle
Q(t,t^{\prime})$ $\displaystyle=\langle e^{-\nu(t)}n(t)n(t^{\prime})\rangle,$
(30b) $\displaystyle C(t,t^{\prime})$ $\displaystyle=\langle
n(t)n(t^{\prime})\rangle$ (30c) $\displaystyle p(t)$ $\displaystyle=\langle
e^{-\nu(t)}n(t)\rangle,$ (30d) $\displaystyle m(t)$ $\displaystyle=\langle
n(t)\rangle,$ (30e) $\displaystyle\Lambda_{P}(t,t^{\prime})$
$\displaystyle=\mbox{$\frac{\sigma^{2}}{2}$}C(t,t^{\prime}),$ (30f)
$\displaystyle\Lambda_{Q}(t,t^{\prime})$
$\displaystyle=\mbox{$\frac{\sigma^{2}}{2}$}\left[2\gamma
Q(t^{\prime},t)-2(1+\gamma)C(t,t^{\prime})]\right],$ (30g)
$\displaystyle\Lambda_{C}(t,t^{\prime})$
$\displaystyle=\mbox{$\frac{\sigma^{2}}{2}$}\left[P(t,t^{\prime})+(1+\gamma)[-Q(t,t^{\prime})-Q(t^{\prime},t)+C(t,t^{\prime})]\right],$
(30h) $\displaystyle\lambda_{p}(t)$ $\displaystyle=sm(t),$ (30i)
$\displaystyle\lambda_{m}(t)$ $\displaystyle=sp(t)-2sm(t),$ (30j)

where the expectation values are taken with respect to the single-species
action $S_{1}$

$\displaystyle\langle
O(\nu(t),n(t^{\prime}))\rangle\equiv\overline{\left\\{\frac{\int{\mathscr{D}}n\int{\mathscr{D}}\nu\;e^{-S_{1}}O(\nu(t),n(t^{\prime}))}{\int{\mathscr{D}}n\int{\mathscr{D}}\nu\;e^{-S_{1}}}\right\\}},$

for a general observable $O(\nu(t),n(t^{\prime}))$. This reduction assumes
only $N\to\infty$ with no further assumption regarding large populations.
Note, however, that evaluation of the single-species expectation values
requires evaluating a full 1D path integral over the effective fields $n(t)$
and $\nu(t)$.

Since it does not require large populations, this theory generalizes the DMFT
obtained previously [12], and still includes rare events. However, it has a
problem: we have derived this theory ignoring the $v_{ij}$ terms in (21), so
that $S_{ij}=w_{ij}$. The distribution of $w_{ij}$ assumed in (16) follows
that of [12]. This is equivalent to assuming
$w_{ij}=\mbox{$\frac{1}{N}$}[s+\sqrt{N}\sigma W_{ij}]$ where
$W_{ij}\sim\mathscr{N}(0,1)$, i.e. $\langle W_{ij}\rangle=0,\langle
W_{ij}^{2}\rangle=1$. But now we see a problem: for large $N$, $w_{ij}$ will
be dominated by the random part,
$w_{ij}\approx\mbox{$\frac{\sigma}{\sqrt{N}}$}W_{ij}$, and therefore it will
be negative for nearly half of the pairs, contradicting the fact that
$w_{ij}\geq 0$ (unless $\sigma=0$).

We face a similar problem if we begin with (21). Although we can let
$S_{ij}=w_{ij}-v_{ij}$ take both signs, the term involving $v_{ij}$ will lead
to the same problem. This result is not specific to the choice of particular
reactions to compose $S_{ij}$. It is a simple consequence of the fact that all
the reaction rates contribute positively to the noise, so there can be no
compensation from positive and negative quantities needed for the above
scaling.

It follows that to obtain a dynamical mean-field theory for a reaction
network, with the above scaling for the variance compared to the mean of the
couplings, the aforementioned simplification of the noise is essential. Since
this requires the small-$\nu$ expansion, to complete the derivation we now
show that by expanding $S_{1}$ we can recover the DMFT of [12], and the
closely related one of [11].

Ignoring terms of $\mathscr{O}(1/N)$, we have

$\displaystyle S_{1}$
$\displaystyle=-\int_{t,t^{\prime}}[\Lambda_{P}(t,t^{\prime})+\Lambda_{Q}(t,t^{\prime})+\Lambda_{C}(t,t^{\prime})]n(t)n(t^{\prime})-\int_{t}[\lambda_{p}(t)+\lambda_{m}(t)]n(t)+\mathscr{O}(\nu^{3})$
$\displaystyle+\int_{t}\nu(t)[-rn(t)(1-n(t)/\kappa)-\lambda+\partial_{t}n(t)+n(t)\lambda_{p}(t)]-\int_{t,t^{\prime}}\nu(t)\chi(t,t^{\prime})n(t)n(t^{\prime})$
$\displaystyle-\mbox{$\frac{1}{2}$}\int_{t,t^{\prime}}\nu(t)\tilde{K}(t,t^{\prime})\nu(t^{\prime})+\mathscr{O}(\nu^{3})$

where

$\displaystyle\chi(t,t^{\prime})=-2\Lambda_{P}(t,t^{\prime})-\Lambda_{Q}(t,t^{\prime})$

and where the noise kernel is

$\displaystyle\tilde{K}(t,t^{\prime})$
$\displaystyle=\delta(t-t^{\prime})\left[rn(1+2n/\kappa)+\lambda+n\lambda_{p}-n(t)\int_{t^{\prime\prime}}\chi(t,t^{\prime\prime})n(t^{\prime\prime})\right]+2\Lambda_{P}(t,t^{\prime})n(t)n(t^{\prime}),$

which is coloured and multiplicative, in general.

Now we seek a regime in which the small-$\nu$ expansion is justified. Above we
saw that the correct scaling ansatz is $n\sim\kappa,t\sim 1/r$. We obtain the
scaling of $\nu$ by equating the linear and quadratic terms, giving $\nu
trn\sim t^{2}\nu^{2}\tilde{K}$, i.e. $\nu\sim r^{2}\kappa/\tilde{K}$. In
$\tilde{K}$, we have the time-local term, scaling as $(1/t)rn\sim
r^{2}\kappa$, and the nonlocal term, scaling as
$\sigma^{2}n^{4}\sim\sigma^{2}\kappa^{4}$. The nonlocal term will dominate if
$\sigma^{2}\kappa^{4}\gg r^{2}\kappa$, i.e. $\kappa^{3}\gg r^{2}/\sigma^{2}$.
Assuming this, we obtain $\nu\sim r^{2}/(\sigma^{2}\kappa^{3})$. Thus if
$\sigma\sim 1,r\sim 1$, and $\kappa\gg 1$, then $\nu$ will indeed be small,
and moreover the kernel will be dominated by its nonlocal part

$\displaystyle K(t,t^{\prime})$
$\displaystyle=2\Lambda_{P}(t,t^{\prime})n(t)n(t^{\prime})$

We now replace $\tilde{K}(t,t^{\prime})$ by $K(t,t^{\prime})$. We will check a
posteriori that the time-local response term in $\tilde{K}$ is indeed
subdominant.

Now perform a Hubbard-Stratonovich transformation

$\displaystyle
e^{\frac{1}{2}\int_{t,t^{\prime}}K(t,t^{\prime})\nu(t)\nu(t^{\prime})}\propto|\det
K|^{-1/2}\int{\mathscr{D}}\xi\;e^{-\frac{1}{2}\int_{t,t^{\prime}}K^{-1}(t,t^{\prime})\xi(t)\xi(t^{\prime})}e^{\int_{t}\xi(t)\nu(t)}$

so that $S_{1}$ is now linear in $\nu$. Integrating $\nu$ out, we obtain the
nonlocal DMFT equation

$\displaystyle\partial_{t}n=[rn(1-n/\kappa)+\lambda-n\lambda_{p}+\xi]|_{t}+\int_{t^{\prime}}\chi(t,t^{\prime})n(t)n(t^{\prime})$

where the correlators that appear here must be determined self-consistently
using (30). Expanding the correlators in small $\nu$ we find

$\displaystyle\Lambda_{P}(t,t^{\prime})$
$\displaystyle=\mbox{$\frac{\sigma^{2}}{2}$}C(t,t^{\prime})$
$\displaystyle\Lambda_{Q}(t,t^{\prime})$
$\displaystyle=\mbox{$\frac{\sigma^{2}}{2}$}\left[-2C(t,t^{\prime})-2\gamma\langle
n_{t}n_{t^{\prime}}\big{(}\nu_{t^{\prime}}+\ldots\big{)}\rangle_{\xi}\right]$

so that

$\displaystyle\chi(t,t^{\prime})$
$\displaystyle=-\mbox{$\frac{\sigma^{2}}{2}$}\left[2C(t,t^{\prime})-2C(t,t^{\prime})-2\gamma\langle
n_{t}n_{t^{\prime}}\nu_{t^{\prime}}\rangle_{\xi}+\ldots\right]$
$\displaystyle=\sigma^{2}\gamma\langle
n_{t}n_{t^{\prime}}\nu_{t^{\prime}}\rangle_{\xi}+\ldots$

For a causal solution we should have $\chi(t,t^{\prime})=0$ for
$t^{\prime}>t$. If the average were computed using the bare propagator, this
would follow from the causality rules for the Doi-Peliti path integral [33].
We assume that it extends to DMFT. To rewrite $\chi$ in a more convenient
form, note that

$\displaystyle\int{\mathscr{D}}\xi\;e^{-\frac{1}{2}\int_{t,t^{\prime}}K^{-1}(t,t^{\prime})\xi(t)\xi(t^{\prime})}e^{\int_{t}\xi(t)\nu(t)}\nu(t^{\prime\prime})=\int{\mathscr{D}}\xi\;e^{-\frac{1}{2}\int_{t,t^{\prime}}K^{-1}(t,t^{\prime})\xi(t)\xi(t^{\prime})}\frac{\delta}{\delta\xi(t^{\prime\prime})}e^{\int_{t}\xi(t)\nu(t)}$
$\displaystyle\quad=-\int{\mathscr{D}}\xi\;\frac{\delta}{\delta\xi(t^{\prime\prime})}\left[e^{-\frac{1}{2}\int_{t,t^{\prime}}K^{-1}(t,t^{\prime})\xi(t)\xi(t^{\prime})}\right]e^{\int_{t}\xi(t)\nu(t)}$
$\displaystyle\quad=\int{\mathscr{D}}\xi\;e^{-\frac{1}{2}\int_{t,t^{\prime}}K^{-1}(t,t^{\prime})\xi(t)\xi(t^{\prime})}e^{\int_{t}\xi(t)\nu(t)}\int_{t}K^{-1}(t,t^{\prime\prime})\xi(t).$

This implies that $\nu(t^{\prime\prime})$ has the same distribution as
$\int_{t}K^{-1}(t,t^{\prime\prime})\xi(t)$ and hence

$\displaystyle\chi(t,t^{\prime})=\sigma^{2}\gamma\left\langle
n_{t}n_{t^{\prime}}\int_{t^{\prime\prime}}K^{-1}(t^{\prime},t^{\prime\prime})\xi(t^{\prime\prime})\right\rangle_{\xi}.$

If we rescale $\xi(t)=\sigma n(t)\tilde{\xi}(t)$, we get

$\displaystyle\chi(t,t^{\prime})=\sigma\gamma\left\langle
n_{t}\int_{t^{\prime\prime}}C^{-1}(t^{\prime},t^{\prime\prime})\tilde{\xi}(t^{\prime\prime})\right\rangle_{\xi}.$

The DMFT equation becomes

$\displaystyle\partial_{t}n_{t}=n_{t}\left[1-n_{t}-s\langle
n_{t}\rangle+\sigma\tilde{\xi}_{t}+\int_{0}^{t}dt^{\prime}\chi(t,t^{\prime})n_{t^{\prime}}\right]+\lambda,$
(31)

where $\langle\tilde{\xi}_{t}\tilde{\xi}_{s}\rangle=C(t,s)$. The obtained Eq.
31 is the same as the DMFT equation of [58], where we identify $\chi$ as the
response function.

We should check that the time-local response term in $\tilde{K}$ is
subdominant. It scales as $(1/t)t\chi
n^{2}\sim\sigma^{2}n^{4}\nu\sim\sigma^{2}\kappa^{4}r^{2}(\sigma^{2}\kappa^{3})^{-1}\sim
r^{2}\kappa$. This needs to be smaller than $\sigma^{2}\kappa^{4}$, which
requires $\kappa^{3}\gg r^{2}/(\sigma^{2})$, which certainly holds if
$\sigma\sim 1$ and $r\sim 1$, as assumed.

We note in passing that the DMFT can also be derived by first expanding the
Doi-Peliti action to quadratic order in $\nu$, and then performing the random
average and introducing new variables.

##  4 Using the mapping

We have found a chemical reaction network (CRN) whose large-scale behavior
maps to the generalized Lotka-Volterra (GLV) model recently considered in
theoretical ecology. By doing so, we have discussed the underlying assumptions
for the GLV model to be valid. Assuming conditions under which the reduction
from a chemical reaction network (CRN) to the generalized Lotka-Volterra (GLV)
model is valid, we ask how (i) this mapping may be useful in the ecology
context; and (ii) how the extensive knowledge of the GLV may find application
in the study of CRNs.

###  4.1 From reaction networks to ecosystems

We discuss here how the CRN theory might prove useful in the ecology context.
First, as emphasized above, the existence of a Lyapunov function is a generic
feature of CRNs in the deterministic limit. For the GLV model, a Lyapunov
function was known [14]; from the mapping we see that this is not an accident.

Second, there is debate about the correct numerical treatment of Langevin
equations with multiplicative noise [59, 60], such as the GLV model. In
particular, Ref. [52] discusses how delicate the treatment of the immigration
reaction in numerical simulations is. The presence of an immigration term is
important to prevent extinctions. The noise, and the immigration reaction in
particular, are important when a species abundance is small. In this article,
we demonstrated that the reduction from a full path integral to a Langevin
equation is not generally justified in this regime. Thus, the GLV CRN can be
used to benchmark any treatment of the Langevin model. To this end the
extensive literature on approximation and simulation of CRNs may be useful
[61, 62, 63, 64].

Third, we showed above that one can obtain a DMFT for the GLV CRN that does
not require any assumption about large populations. This DMFT, while
sophisticated, could be used to understand rare events in ecosystems, such as
extinctions, where the Hamiltonian framework has already proved useful [65,
66].

###  4.2 From ecosystems to reaction networks

We now consider, conversely, what the study of the GLV model can tell us about
chemical reaction networks? A key outcome of recent works is the phase diagram
of the GLV model as a function of typical interaction $s$ [13, 14, 58],
interaction fluctuation strength $\sigma$ [13, 14, 58], interaction symmetry
$\gamma$ [58], and noise strength $\omega$ [52, 53].

In the limit of high disorder (large $\sigma$), some species grow without
bound and the system does not reach a steady state. This ‘unbounded growth’
phase is however non-physical and a pathology of the model. Indeed, a large
variance $\sigma$ implies the presence of strongly cooperative interactions
between some of the species, that are stronger than the species’ logistic
saturation and lead to unbounded growth. This non-physical behavior can be
cured by a saturation stronger than quadratic [58].

For smaller $\sigma$ there are three dynamically stable phases, characterized
by the number and type of the attractor states. The simplest is the single
equilibrium phase, which holds for large enough noise strength $\omega$ [13,
14]. In that phase, the system relaxes toward a unique stationary state with
fluctuating equilibrium dynamics, irrespective of the initial condition. When
decreasing further the demographic noise $\omega$, or considering larger
populations, two new phases appear, both containing multiple attractors of the
dynamics [13, 14, 52]. In these multiple-attractor phases, the equilibrium
state of the system depends on the initial condition. The most remarkable
phase is the one discovered at low $\omega$ in models with symmetric
interactions $\gamma=1$ [52]. In this region of the phase diagram,
interactions play an important role, and their disordered nature gives rise to
a so-called Gardner, marginally stable, phase. This phase was first discovered
in the context of mean-field spin glasses [29], and recently shown to play an
important role in the physics of jamming and amorphous materials [67, 68, 69].
In the Gardner phase, there is an exponential (in system size) number of
equilibrium states, with a hierarchical organization in configuration space.
Importantly, the states are marginally stable, i.e. poised at the edge of
stability [14], and separated by arbitrarily small barriers. The presence of a
fractal landscape of marginally stable minima has profound implications on the
dynamics of the system. In spin and structural glasses, this gives rise to
large responses (avalanches) to small perturbations, aging dynamics [70], as
well as rejuvenation and memory effects [71].

Let us suppose, for the sake of argument, that this phase structure carries
over to large biochemical reaction networks; later, we will take a critical
stance. The prospect of a Parisi-Gardner phase in biochemical reaction
networks is tantalizing. Indeed, the notion of a landscape with many
attractors has a long history in biology, going back to Waddington’s notion of
epigenetic landscapes [72, 48]. While many authors have argued for a multiple-
attractor landscape in biology [46, 47, 48, 37, 38], explicit models of
multiple attractors usually focus on the simple case of two, or few,
attractors. This is the case of the Schlögl model [73, 74], which is perhaps
the simplest reaction network to support multiple equilibria, and many of its
properties can be computed exactly [74]. However, in a Parisi-Gardner phase,
the number of equilibria scales exponentially with $N$, and one must account
for the statistics of the barriers. We are not aware of any studies on
chemical reaction networks where the notion of exponentially many attractor
states has been entertained. Therefore it is crucial to see if the GLV DMFT,
where the existence of this phase can be established, has any analogue in real
CRNs.

##  5 Discussion: Towards a new universality class of DMFT

The success of DMFT in understanding large ecosystems opens the possibility to
describe the behaviour of large chemical reaction networks with similar
methods. However, while the general goal is clear, crucial features of
chemical reactions give rise to technical obstacles to overcome. Several of
these have been alluded to above, and we add another below. Addressing them in
turn, we argue that the current ecology DMFT is inadequate to capture generic
behaviour in real CRNs. This presents an important opportunity for future
work.

###  5.1 Incompatibility with local detailed balance

In extant life, many reactions are catalyzed so that the reverse reactions
rarely occur, and are often neglected in modeling. Indeed, since life involves
energy fluxes it must exist far from thermal equilibrium and it is expected
that its current state is strongly time-asymmetric, as is now being confirmed
experimentally [75, 76]. However, a complete theoretical description for
chemical reactions must include each reaction as part of a pair in local
detailed balance.

We noted above that the birth-death reactions in the GLV model occur in
detailed balance, but that the interactions in general do not. Let us now ask
how one could modify the DMFT to include the reverse reactions.

Consider the reaction $A_{j}+A_{l}\overset{\;v_{jl}}{\rightarrow}2A_{j}+A_{l}$
and let us now add the reverse reaction
$2A_{j}+A_{l}\overset{\;\tilde{v}_{jl}}{\rightarrow}A_{j}+A_{l}$. Local
detailed balance requires that the reaction rates satisfy a constraint of the
form (12), i.e.
$v_{jl}/\tilde{v}_{jl}=e^{-(\vec{q}_{\alpha}-\vec{p}_{\alpha})\cdot\vec{G}}$.
Without loss of generality we can adopt the form from transition state theory
used for CRNs (13). For the reactions with rates $v_{jl}$ and
$\tilde{v}_{jl}$, let us denote the corresponding ‘activation energies’ as
$G^{v}_{jl}$. For the GLV model we need also the reactions
$A_{j}+A_{l}\overset{\;w_{jl}}{\rightarrow}A_{l}$ and
$A_{j}+A_{l}\overset{\;\tilde{w}_{jl}}{\leftarrow}A_{l}$. Altogether these 4
reactions will contribute to the Hamiltonian

$\displaystyle
H_{ij}=k_{0}e^{-G_{ij}^{v}}(e^{\nu_{i}}-1)x_{i}x_{j}+k_{0}e^{-G_{ij}^{w}}(e^{-\nu_{i}}-1)x_{i}x_{j}+k_{0}e^{-G_{ij}^{v}}(e^{-\nu_{i}}-1)x_{i}^{2}x_{j}+k_{0}e^{-G_{ij}^{w}}(e^{\nu_{i}}-1)x_{j}$
(32)

where $x_{i}=n_{i}/e^{-G_{i}}$ in appropriate units. It will suffice to
consider this to leading order in $\nu$:

$\displaystyle
H_{ij}=k_{0}\nu_{i}x_{j}\left[e^{-G_{ij}^{v}}x_{i}-e^{-G_{ij}^{w}}x_{i}-e^{-G_{ij}^{v}}x_{i}^{2}+e^{-G_{ij}^{w}}\right]+\mathscr{O}(\nu^{2})$
(33)

From here we can appreciate the problem. Earlier we combined the rates so that
$v_{ij}-w_{ij}$ takes both signs and can therefore have a Gaussian
distribution. This is not possible here, because the reverse reactions have
different powers of $x_{i}$. Although for the $x_{i}^{2}x_{j}$ we could add
another reaction with the same factor, this is not possible with the $x_{j}$
reaction; physically this corresponds to the fact that extinction reactions
can never be in local detailed balance. It is thus not possible to combine
multiple reaction rates so that their combination has a Gaussian distribution.

In order to derive a true CRN DMFT, one must therefore deal with non-Gaussian
distributions and average terms of the form

$\displaystyle\overline{e^{\int_{t}e^{-G_{ij}^{v}}A_{ij}}}$ (34)

over $G_{ij}^{v}$, where for example
$A_{ij}=k_{0}(x_{i}x_{j}(e^{\nu_{i}}-1)+x_{j}(e^{-\nu_{i}}-1))$. This is the
moment-generating function of the random variable $y=e^{-G_{ij}^{v}}$,
evaluated at $\int_{t}A_{ij}$. Further development of the theory uses the
cumulant-generating function

$\displaystyle L(X)=\log\overline{\exp\left(yX\right)}$ (35)

evaluated at $X=\int_{t}A_{ij}$. If $y$ were Gaussian, then $L$ would be a
quadratic polynomial in $\int_{t}A_{ij}$; this leads to the two-time
correlation and response quantities defined above. One might hope that an
appropriate distribution would lead to a finite low-order polynomial for the
non-Gaussian case. However, it is known from probability theory that no such
distribution exists: there are no distributions whose cumulant generating
function is a finite polynomial of order greater than two [77]. Thus $L$ will
have terms of all orders in $\int_{t}A_{ij}$, leading to a theory with
quantities depending on arbitrarily many different times. This appears
completely intractable.

Summarizing this subsection, in a theory with local detailed balance, every
reaction appears along with its reverse. Once these reversed reactions are
included in the GLV model, it is impossible to group together multiple
reactions in order for a generalized rate to have a Gaussian distribution.
Non-Gaussian distributions (which are necessary to respect positivity of
reaction rates) lead to an intractable DMFT.

###  5.2 Non-positivity of reaction rates

We have seen that imposing local detailed balance leads to problems because of
the positivity of reaction rates. The latter can happen even without imposing
local detailed balance. Indeed, the grouping together of multiple reactions to
form an effective rate $v_{ij}-w_{ij}$ that takes both signs, and therefore
can be given a Gaussian distribution, is a very special feature of the GLV
model. As shown above, it does not even hold beyond the leading order in $\nu$
(i.e. system size). Thus the requirement that reaction rates are positive
creates an even more fundamental roadblock to construction of DMFT by path
integrals than local detailed balance. A similar difficulty appears in random
language models, where positivity of weights leads to the same issue, in a
somewhat simpler context [78]. A solution to this problem might therefore have
many applications.

In real ecosystems too, going beyond the Gaussian approximation is crucial to
describe purely competitive systems and overcome issues related to both
positive and negative interactions that arise from Gaussian distributions with
finite variance.

It is worth noting that the dynamical cavity method, which can also be used to
construct DFMTs, is not subject to the same limitations on non-Gaussian
distributions. The cavity method and the path-integral method are
complementary: while the former works directly in the limit $N\to\infty$, the
latter can be used at general $N$, and can be analyzed with the sophisticated
machinery of the renormalization group. In previous work [79], the cavity
method was used with non-Gaussian distributions for interactions in the
MacArthur resource model [80, 81], another key model of theoretical ecology.

###  5.3 Absence of conserved quantities in stoichiometry

Finally, the GLV model, when considered as a CRN, lacks the most conspicuous
feature of real CRNs: nontrivial stoichiometry, as characterized by the
stoichiometric matrix $S_{\alpha j}=q_{\alpha j}-p_{\alpha j}$. Since (non-
nuclear) chemical reactions preserve the number of each element, stoichiometry
gives a large number of constraints on possible reactions for a given set of
species. However, it is easy to see that reactions of the form (17),(18),(19)
cannot have conserved quantities. For example, if $A+B\rightarrow A$, then $B$
cannot carry any conserved ‘elements’. This should be compared with a typical
real reaction, such as Ca(OH)2 \+ CO${}_{2}\rightarrow$ CaCO3 \+ H2O, which
conserves the number of calcium, oxygen, hydrogen, and carbon atoms.

In some contexts, it can be permissible to consider chemical reactions that do
not conserve elements. In particular, in biochemistry, reactions are often
expressed up to numbers of water atoms, since water is abundant in the
cellular environment. If the GLV model is to be considered as a CRN, then it
too must be interpreted in a context where elements do not give stringent
conservation laws.

Work on random CRNs that feature conservation laws would thus be welcome. It
was already shown in [82] that adding conserved quantities (called moieties)
to random CRNs strongly affects the possible growth of the network.

###  5.4 Looking ahead

Summarizing this section, we have seen that the GLV model lacks several
crucial features of real CRNs. This of course does not mitigate the importance
of the GLV model in ecology, but instead points to its location in the phase
space of DMFT as firmly in the emergent, mesoscopic regime. Its utility there
has been demonstrated: in particular, the DMFT has shown that this
nonequilibrium system can exist in different phases, and universality with
respect to microscopic parameters has been shown. Tackling the challenges
exposed above will prove useful to describe emergent behavior in many-species
chemical reaction networks.

##  6 Conclusion

Dynamical mean-field theory has recently shown its power to identify dynamical
phases in a model ecosystem, the generalized Lotka-Volterra (GLV) model. In
this article, we explored to what extent these recent successes can help the
development of a theoretical framework for chemical reaction networks (CRN),
focusing on the limit of many chemical species.

To guide the mapping between ecology and CRN at the DMFT level, we have shown
how one can construct a CRN whose infinite-size limit precisely recovers the
GLV model using the Hamiltonian formalism of CRNs. We critically discussed the
conditions under which this limit is obtained, and also addressed the large-
but-finite-system size limit in which a Langevin description is obtained. We
found that it is not possible to obtain the GLV Langevin model without quite
stringent, and in general unjustified, assumptions. Then, we showed how the
GLV DMFT can be obtained from the Doi-Peliti path integral, and again
critically examined necessary assumptions. Finally, we discussed how the GLV-
CRN mapping could be used, in both directions. Several specific features of
the GLV DMFT suggest that it is not directly applicable to real microscopic
CRNs: it does not have all reactions in local detailed balance, and there are
no conserved moieties. These results suggest that new classes of DMFT should
be developed to investigate emergent behavior in chemical reactions networks,
which are relevant to describe living matter.

We are grateful to Masanari Shimada for a critical reading of the manuscript.
This work was funded by NSERC Discovery Grant RGPIN-2020-04762 (to E. De
Giuli), a Herchel Smith Fellowship, University of Cambridge (to C. Scalliet),
and a Ramon Jenkins Research Fellowship from Sidney Sussex College, Cambridge
(to C. Scalliet).

##  Appendix A Doi-Peliti path integral

###  A.1 Doi formalism

Doi [31] observed that the step operators in (2), together with associated
numerical factors $n_{j}!/(n_{j}-p_{\alpha j})!$, have a natural
representation in terms of bosonic creation and annihilation operators, thus
connecting (2) to quantum mechanics through second quantization. The vacuum,
with no molecules, is denoted by $|0\rangle$, and by definition satisfies
$\hat{a}_{i}|0\rangle=0$ for all $i$. A molecule of species $i$ is created by
$\hat{a}^{\dagger}_{i}$. The operators are normalized such that

$\displaystyle\hat{a}_{i}|\vec{n}\rangle$
$\displaystyle=n_{i}|\vec{n}-\vec{1}^{i}\rangle,$ (36)
$\displaystyle\hat{a}^{\dagger}_{i}|\vec{n}\rangle$
$\displaystyle=|\vec{n}+\vec{1}^{i}\rangle,$ (37)

and

$\displaystyle[\hat{a}_{i},\hat{a}^{\dagger}_{j}]=\delta_{ij},\qquad[\hat{a}_{i},\hat{a}_{j}]=[\hat{a}^{\dagger}_{i},\hat{a}^{\dagger}_{j}]=0.$
(38)

From these relations one can prove that [34]

$\displaystyle
e^{\hat{a}_{i}}f(\hat{a}_{i},\hat{a}^{\dagger}_{i})=f(\hat{a}_{i},\hat{a}^{\dagger}_{i}+1)e^{\hat{a}_{i}},$
(39)

for arbitrary functions $f$.

The state of the entire system is then encoded in the state vector

$\displaystyle|\phi(t)\rangle=\sum_{\\{\vec{n}\\}}P(\vec{n},t)\prod_{i}(\hat{a}^{\dagger}_{i})^{n_{i}}|0\rangle$
(40)

The master equation is then equivalent to an imaginary-time Schrödinger
equation

$\displaystyle\frac{\partial}{\partial
t}|\phi(t)\rangle=\hat{H}|\phi(t)\rangle,$ (41)

with a non-Hermitian quasi-Hamiltonian (or Liouvillian)

$\displaystyle\hat{H}(\vec{a}^{\dagger},\vec{a})=\sum_{\alpha}\tilde{k}_{\alpha}\left[\prod_{j}(\hat{a}^{\dagger}_{j})^{q_{\alpha
j}}-\prod_{j}(\hat{a}^{\dagger}_{j})^{p_{\alpha
j}}\right]\prod_{i}\hat{a}_{i}^{p_{\alpha i}}$ (42)

Notice that the entropic factors $n_{j}!/(n_{j}-p_{\alpha j})!$ are absent
from Eqs.(41),(3).

In Eq. (42), the terms involving $(\hat{a}^{\dagger}_{j})^{q_{\alpha j}}$
straightforwardly represent the effect of each reaction on molecule $j$, while
those involving $(\hat{a}^{\dagger}_{j})^{p_{\alpha j}}$ are instead present
to ensure that probability normalization is preserved.

Eq. (41) is solved by $|\phi(t)\rangle=e^{\hat{H}t}|\phi(0)\rangle$. The Born
rule does not apply to this construction. Instead, observables are extracted
with the coherent state

$\displaystyle\langle\phi_{0}|\equiv\langle 0|\Pi_{i}\;e^{\hat{a}_{i}},$ (43)

which satisfies $\langle\phi_{0}|\prod_{i}n_{i}\rangle=1$ for any $\vec{n}$,
using Eq. (39). This implies
$\langle\phi_{0}|\hat{a}^{\dagger}_{i}=\langle\phi_{0}|$, since
$\langle\phi_{0}|\hat{a}^{\dagger}_{i}|n_{i}\rangle=\langle\phi_{0}|n_{i}+1\rangle=1=\langle\phi_{0}|n_{i}\rangle$.
The expectation value of an observable $\mathscr{O}$ is given by projection
onto the coherent state:

$\displaystyle\langle\hat{\mathscr{O}}\rangle=\langle\phi_{0}|\hat{\mathscr{O}}|\phi(t)\rangle$
(44)

For example, for the number operator $\hat{a}^{\dagger}_{i}\hat{a}_{i}$ we
have

$\displaystyle\langle n_{i}(t)\rangle$
$\displaystyle=\langle\phi_{0}|\hat{a}^{\dagger}_{i}\hat{a}_{i}|\phi(t)\rangle=\langle\phi_{0}|\hat{a}^{\dagger}_{i}\hat{a}_{i}\Sigma_{\\{\vec{n}\\}}P(\vec{n},t)|\vec{n}\rangle=\langle\phi_{0}|\Sigma_{\\{\vec{n}\\}}n_{i}P(\vec{n},t)|\vec{n}\rangle=\Sigma_{\\{\vec{n}\\}}n_{i}P(\vec{n},t),$
(45)

which is indeed the expectation value of $n_{i}$. Let us check normalization

$\displaystyle\Sigma_{\\{\vec{n}\\}}P(\vec{n},t)$
$\displaystyle=\langle\phi_{0}|\phi(t)\rangle=\langle
0|\Pi_{i}\;e^{\hat{a}_{i}}|e^{\hat{H}(\vec{a}^{\dagger},\vec{a})t}|\phi(0)\rangle=\langle
0|e^{\hat{H}(\vec{a}^{\dagger}+1,\vec{a})t}\Pi_{i}\;e^{\hat{a}_{i}}|\phi(0)\rangle$
$\displaystyle=\langle
0|e^{\hat{H}(1,\vec{a})t}\Pi_{i}\;e^{\hat{a}_{i}}|\phi(0)\rangle.$ (46)

This will equal unity if probability is initially normalized and

$\displaystyle\hat{H}(1,\vec{a})=0,$ (47)

which is indeed satisfied by Eq. (42).

###  A.2 Path integral formulation

The path-integral representation is built just as in quantum mechanics. In the
absence of sources, the partition function

$\displaystyle Z=\langle
0|e^{\hat{H}(\vec{a}^{\dagger}+1,\vec{a})t_{f}}\Pi_{i}\;e^{\hat{a}_{i}}|\phi(0)\rangle,$
(48)

is equal to unity $Z=1$ due to conservation of probability. This is
reminiscent of the De Dominicis-Janssen field theory [39, 40]. The path
integral representation is built by dividing the time interval $[0,t_{f}]$
into $N\rightarrow\infty$ slices, and introducing the coherent state
resolution of unity at each time slice. Following this standard procedure, one
obtains the path-integral representation of the partition function [34, 32,
33]

$\displaystyle Z=\int{\mathscr{D}}\phi\int{\mathscr{D}}\phi^{*}\;e^{-S},$ (49)

with

$\displaystyle
S=\sum_{i}\left[-\phi_{i}(t_{f})-n^{0}_{i}\log\phi^{*}_{i}(0)+\phi_{i}(0)\phi^{*}_{i}(0)\right]+\int_{0}^{t_{f}}dt\left[\sum_{i}\phi^{*}_{i}\partial_{t}\phi_{i}-H(\vec{\phi}^{*},\vec{\phi})\right].$
(50)

Here the Hamiltonian becomes a function, evaluated on the fields
$\vec{\phi}^{*}$ and $\vec{\phi}$ rather than $\vec{a}^{\dagger}$ and
$\vec{a}$. This is for initial conditions $n_{i}(0)=n^{0}_{i}$; if instead we
have Poisson-distributed initial conditions, then
$n^{0}_{i}\log\phi^{*}_{i}(0)$ becomes $n^{0}_{i}(\phi^{*}_{i}(0)-1)$.

The doubling of degrees-of-freedom in Eq. (50) is characteristic of dynamic
problems [39, 40, 34]. Loosely, $\phi$ can be considered as the classical part
and $\phi^{*}$ as the quantum part of the concentration field [34]. However,
beyond the mean-field limit, the relationship between $\phi,\phi^{*}$ and the
true density field is very subtle [34, 35, 36, 37]. For example, the Langevin
equation obtained from the Doi-Peliti action in a naive semi-classical limit
has imaginary noise. This can be avoided by a Hopf-Cole transformation
$\phi^{*}=e^{\nu},\phi=ne^{-\nu}$, which makes the number operator
$a^{\dagger}a\to\phi^{*}\phi=n$.

The Doi-Peliti action (50) is an exact rewriting of the original master
equation and can be used to obtain the behaviour of arbitrary observables. It
is the basis for renormalization group treatments of reaction-diffusion and
related systems, when diffusion terms are added [83, 33].

###  A.3 Generating function

In order to extract observables, sources are added to $Z$ by introducing the
generating function

$\displaystyle
Z(\vec{z},t)=\sum_{\\{\vec{n}\\}}z_{1}^{n_{1}}z_{2}^{n_{2}}\cdots
z_{N}^{n_{N}}P(\vec{n},t).$ (51)

The partition function is recovered for $\vec{z}=1$, and moments are computed
by derivatives at $\vec{z}=1$. Elgart and Kamenev [41] showed that the path
integral representation of the generating function Eq. (51) is similar to that
of the partition function (no sources, $\vec{z}=1$), with differences only in
the boundary terms

$\displaystyle
S[\phi,\phi^{*};\overline{n},\vec{z}]=\sum_{i}\left[\phi_{i}(0)\phi^{*}_{i}(0)-z_{i}\phi_{i}(t_{f})-n^{0}_{i}\log\phi^{*}_{i}(0)\right]+\int_{0}^{t_{f}}dt\left[\sum_{i}\phi^{*}_{i}\partial_{t}\phi_{i}-H(\vec{\phi}^{*},\vec{\phi})\right],$
(52)

where $\overline{n}$ is the mean trajectory. As shown by Elgart and Kamenev
[41], values of $\vec{z}$ different from unity are necessary to understand
rare fluctuations.

It is convenient to consider the Hopf-Cole transformation
$\phi_{j}=n_{j}e^{-\nu_{j}},\phi^{*}_{j}=e^{\nu_{j}}$. This is a canonical
transformation with unit Jacobian. The kinetic term in the action transforms
as

$\displaystyle\int_{t}\phi^{*}_{j}\partial_{t}\phi_{j}=\int_{t}e^{\nu_{j}}e^{-\nu_{j}}\left[\partial_{t}n_{j}-n_{j}\partial_{t}\nu_{j}\right]=n_{j}(1-\nu_{j})|_{0}^{t_{f}}+\int_{t}\nu_{j}\partial_{t}n_{j}$
(53)

The first term can be absorbed into the boundary conditions, i.e. we now have

$\displaystyle
S[n,\nu;\overline{n},\vec{z}]=\sum_{j}\left[\phi_{j}(0)\phi^{*}_{j}(0)-z_{j}\phi_{j}(t)-n^{0}_{j}\log\phi^{*}_{j}(0)\right]+\sum_{j}n_{j}(1-\nu_{j})|_{0}^{t_{f}}+\int_{0}^{t_{f}}dt\left[\vec{\nu}\cdot\partial_{t}\vec{n}-H(\vec{n},\vec{\nu})\right]$
$\displaystyle=\sum_{j}\left[n_{j}(0)-z_{j}n_{j}(t_{f})e^{-\nu_{j}(t_{f})}-n^{0}_{j}\nu_{j}(0)+n_{j}(t_{f})(1-\nu_{j}(t_{f}))-n_{j}(0)(1-\nu_{j}(0))\right]+\int_{0}^{t_{f}}dt\left[\vec{\nu}\cdot\partial_{t}\vec{n}-H(\vec{n},\vec{\nu})\right]$
$\displaystyle=\underbrace{\sum_{j}\left[n_{j}(0)\nu_{j}(0)+n_{j}(t_{f})[-z_{j}e^{-\nu_{j}(t_{f})}+1-\nu_{j}(t_{f})]-n^{0}_{j}\nu_{j}(0)\right]}_{S_{BC}}+\int_{0}^{t_{f}}dt\left[\vec{\nu}\cdot\partial_{t}\vec{n}-H(\vec{n},\vec{\nu})\right]$
(54)

The particle number statistics are extracted from the generating function by
contour integrals

$\displaystyle P(\vec{m},t_{f})=\prod_{j}\oint\frac{dz_{j}}{2\pi
i}\frac{1}{z_{j}^{m_{j}+1}}Z(\vec{z},t_{f})$ (55)

For example, from $P(\vec{m},t_{f})$ the marginal distribution of the $j^{th}$
species is

$\displaystyle\rho_{j}(m_{j})$ $\displaystyle=\prod_{l\neq j}\sum_{m_{l}\geq
0}\prod_{k}\oint\frac{dz_{k}}{2\pi
i}\frac{1}{z_{k}^{m_{k}+1}}Z(\vec{z},t_{f})$ (56)
$\displaystyle=\oint\frac{dz_{j}}{2\pi i}\frac{1}{z_{j}^{m_{j}+1}}\prod_{k\neq
j}\oint\frac{dz_{k}}{2\pi
i}\frac{1}{z_{k}}\frac{1}{1-z_{k}^{-1}}Z(\vec{z},t_{f}),$ (57)

where we take an integral contour $|z_{k}|>1$, for all $k\neq j$. Since the
$z_{k}$-dependence in $Z(\vec{z},t_{f})$ is restricted to a term
$e^{z_{k}\phi_{k}(t_{f})}$, we will have integrals of the form $(2\pi
i)^{-1}\oint dze^{z\phi}/(z-1)=\text{Res}[e^{z\phi}/(z-1),z=1]=e^{\phi}$ so
that

$\displaystyle\rho_{j}(m_{j})$ $\displaystyle=\oint\frac{dz_{j}}{2\pi
i}\frac{1}{z_{j}^{m_{j}+1}}Z((1,\ldots,1,z_{j},1,\ldots,1),t_{f})$
$\displaystyle=\oint\frac{dz_{j}}{2\pi
i}\frac{1}{z_{j}^{m_{j}+1}}Z(\vec{1}+\delta z_{j}\hat{e}_{j},t_{f})$ (58)

where $\vec{1}=(1,1,\ldots,1)$ and $\hat{e}_{j}=(0,\ldots,1,0,\ldots,0)$ with
a 1 in the $j^{th}$ position, and $\delta z_{j}=z_{j}-1$. Therefore we need to
know $Z$ in the vicinity of $\vec{z}=\vec{1}$.

##  Appendix B Instantons

As explained in the main text, we will need the stationary points of the
action, i.e. instantons.

For later use, we have $\partial S/\partial
z_{j}=-n_{j}(t_{f})e^{-\nu_{j}(t_{f})}$, and the first variation is

$\displaystyle\delta S[n,\nu;\overline{n},\vec{z}]=\sum_{j}\left[\delta
n_{j}(0)\nu_{j}(0)+\delta\nu_{j}(0)[n_{j}(0)-n^{0}_{j}]+\delta
n_{j}(t_{f})[-z_{j}e^{-\nu_{j}(t_{f})}+1-\nu_{j}(t_{f})]\right]$
$\displaystyle\qquad+\sum_{j}\delta\nu_{j}(t_{f})n_{j}(t_{f})[z_{j}e^{-\nu_{j}(t_{f})}-1]+\int_{0}^{t_{f}}dt\left[\sum_{i}[\delta\nu_{j}\partial_{t}n_{j}+\nu_{j}\partial_{t}\delta
n_{j}]-\delta H(\vec{n},\vec{\nu})\right]$ $\displaystyle=\sum_{j}\left[\delta
n_{j}(0)\nu_{j}(0)+\delta\nu_{j}(0)[n_{j}(0)-n^{0}_{j}]+\delta
n_{j}(t_{f})[-z_{j}e^{-\nu_{j}(t_{f})}+1-\nu_{j}(t_{f})]\right]$
$\displaystyle\qquad+\sum_{j}\delta\nu_{j}(t_{f})n_{j}(t_{f})[z_{j}e^{-\nu_{j}(t_{f})}-1]$
$\displaystyle\qquad+\left[\sum_{j}[\int_{0}^{t_{f}}dt\delta\nu_{j}\partial_{t}n_{j}+\nu_{j}(t_{f})\delta
n_{j}(t_{f})-\nu_{j}(0)\delta n_{j}(0)-\int_{0}^{t_{f}}dt\delta
n_{j}\partial_{t}\nu_{j}-\int_{0}^{t_{f}}dt\delta H(\vec{n},\vec{\nu})\right]$
$\displaystyle=\sum_{j}\left[\delta\nu_{j}(0)[n_{j}(0)-n^{0}_{j}]+\delta
n_{j}(t_{f})[-z_{j}e^{-\nu_{j}(t_{f})}+1]+\delta\nu_{j}(t_{f})n_{j}(t_{f})[z_{j}e^{-\nu_{j}(t_{f})}-1]\right]$
$\displaystyle\qquad+\int_{0}^{t_{f}}dt\sum_{j}\left[\delta\nu_{j}\partial_{t}n_{j}-\delta
n_{j}\partial_{t}\nu_{j}-\frac{\partial H}{\partial n_{j}}\delta
n_{j}-\frac{\partial H}{\partial\nu_{j}}\delta\nu_{j}\right]$ (59)

The stationary solutions must solve

$\displaystyle\frac{\partial n_{j}}{\partial t}$
$\displaystyle=+\frac{\partial H}{\partial\nu_{j}}$ (60)
$\displaystyle\frac{\partial\nu_{j}}{\partial t}$
$\displaystyle=-\frac{\partial H}{\partial n_{j}}$ (61)

with boundary conditions

$\displaystyle 0$ $\displaystyle=n_{j}(0)-n^{0}_{j}$ (62) $\displaystyle 1$
$\displaystyle=z_{j}e^{-\nu_{j}(t_{f})}$ (63) $\displaystyle 0$
$\displaystyle=n_{j}(t_{f})[z_{j}e^{-\nu_{j}(t_{f})}-1]$ (64)

These are obviously degenerate and reduce to $\nu_{j}(t_{f})=\log z_{j}$,
$n_{j}(0)=n^{0}_{j}$.

Since (60),(61) take the form of Hamilton’s equations, they conserve $H$:

$\displaystyle\frac{dH}{dt}$ $\displaystyle=\frac{\partial H}{\partial
t}+\frac{\partial H}{\partial\vec{n}}\cdot\frac{\partial\vec{n}}{\partial
t}+\frac{\partial H}{\partial\vec{\nu}}\cdot\frac{\partial\vec{\nu}}{\partial
t}$ $\displaystyle=\frac{\partial H}{\partial t}+\frac{\partial
H}{\partial\vec{n}}\cdot\frac{\partial H}{\partial\vec{\nu}}-\frac{\partial
H}{\partial\vec{\nu}}\cdot\frac{\partial H}{\partial\vec{n}}$
$\displaystyle=\frac{\partial H}{\partial t}$ (65)

###  B.1 Mean-field instantons minimize the action

Along saddle-points, the action is

$\displaystyle S_{c}$ $\displaystyle=-\sum_{j}n_{j}(t_{f})\log
z_{j}+\int_{0}^{t_{f}}dt\left[\vec{\nu}\cdot\partial_{t}\vec{n}-H(\vec{n},\vec{\nu})\right]$
(66)

Define the function $f(x)=1+e^{x}(x-1)$ and note that the convexity inequality
$x-1\geq-e^{-x}$ implies $f(x)\geq 1+e^{x}(-e^{-x})=0$. Using the form

$\displaystyle H_{\alpha}(\vec{n},\vec{\nu})$
$\displaystyle=\left[e^{\vec{\nu}\cdot(\vec{q}_{\alpha}-\vec{p}_{\alpha})}-1\right]\underbrace{\tilde{k}_{\alpha}\prod_{i}n_{i}^{p_{\alpha
i}}}_{\equiv F_{\alpha}(\vec{n})}$ (67)

we have

$\displaystyle\vec{\nu}\cdot\partial_{t}\vec{n}-H(\vec{n},\vec{\nu})$
$\displaystyle=\vec{\nu}\cdot\frac{\partial H}{\partial\vec{\nu}}-H$
$\displaystyle=\sum_{\alpha}F_{\alpha}(\vec{n})\left[\vec{\nu}\cdot(\vec{q}_{\alpha}-\vec{p}_{\alpha})e^{\vec{\nu}\cdot(\vec{q}_{\alpha}-\vec{p}_{\alpha})}-e^{\vec{\nu}\cdot(\vec{q}_{\alpha}-\vec{p}_{\alpha})}+1\right]$
$\displaystyle=\sum_{\alpha}F_{\alpha}(\vec{n})f\left[\vec{\nu}\cdot(\vec{q}_{\alpha}-\vec{p}_{\alpha})\right]$
$\displaystyle\geq 0,$

where we are assuming that $\vec{\nu}\in\mathbb{R}^{N}$. Therefore

$\displaystyle S_{c}+\vec{n}(t_{f})\cdot\log\vec{z}$ $\displaystyle\geq 0$
(68)

The term in the action $-\vec{n}(t_{f})\cdot\log\vec{z}$ is used to choose
observables. Aside from this term (which has no definite sign since in general
$\vec{z}$ is complex), the action is non-negative. Equality is obtained only
when for each $\alpha$ and each time $t$, we have either
$F_{\alpha}(\vec{n}(t))=0$, or
$\vec{\nu}(t)\cdot(\vec{q}_{\alpha}-\vec{p}_{\alpha})=0$.

Let us note that mean-field solutions correspond to $\nu_{j}(t)=0,\forall j,t$
and $z_{j}=1\forall j$, in which case $S_{c}=0$. Thus we have shown that mean-
field solutions minimize the action.

###  B.2 Fluctuations around instantons

The leading fluctuations around the instanton are obtained by computing the
second variation of $S$:

$\displaystyle\delta^{2}S$ $\displaystyle=\sum_{j}\left[\delta\nu_{j}(0)\delta
n_{j}(0)+\delta
n_{j}(t_{f})z_{j}\delta\nu_{j}(t_{f})e^{-\nu_{j}(t_{f})}+\delta\nu_{j}(t_{f})\delta
n_{j}(t_{f})[z_{j}e^{-\nu_{j}(t_{f})}-1]-\delta\nu_{j}(t_{f})^{2}n_{j}(t_{f})z_{j}e^{-\nu_{j}(t_{f})}\right]$
$\displaystyle\quad+\int_{0}^{t_{f}}dt\left[\delta\vec{\nu}\cdot\partial_{t}\delta\vec{n}-\delta\vec{n}\cdot\partial_{t}\delta\vec{\nu}\right]-\int_{0}^{t_{f}}dt\left[\frac{\partial^{2}H}{\partial\vec{n}\partial\vec{n}}:\delta\vec{n}\delta\vec{n}+2\frac{\partial^{2}H}{\partial\vec{\nu}\partial\vec{n}}:\delta\vec{\nu}\delta\vec{n}+\frac{\partial
H}{\partial\vec{\nu}\partial\vec{\nu}}:\delta\vec{\nu}\delta\vec{\nu}\right]$
$\displaystyle=\sum_{j}\left[\delta\nu_{j}(0)\delta n_{j}(0)+\delta
n_{j}(t_{f})\delta\nu_{j}(t_{f})-\delta\nu_{j}(t_{f})^{2}n_{j}(t_{f})\right]$
$\displaystyle\quad+\int_{0}^{t_{f}}dt\left[2\delta\vec{\nu}\cdot\partial_{t}\delta\vec{n}-\partial_{t}[\delta\vec{n}\cdot\delta\vec{\nu}]\right]-\int_{0}^{t_{f}}dt\left[\frac{\partial^{2}H}{\partial\vec{n}\partial\vec{n}}:\delta\vec{n}\delta\vec{n}+2\frac{\partial^{2}H}{\partial\vec{\nu}\partial\vec{n}}:\delta\vec{\nu}\delta\vec{n}+\frac{\partial
H}{\partial\vec{\nu}\partial\vec{\nu}}:\delta\vec{\nu}\delta\vec{\nu}\right]$
$\displaystyle=\sum_{j}\left[2\delta\nu_{j}(0)\delta
n_{j}(0)-\delta\nu_{j}(t_{f})^{2}n_{j}(t_{f})\right]$
$\displaystyle\quad+\int_{0}^{t_{f}}dt\left[2\delta\vec{\nu}\cdot\partial_{t}\delta\vec{n}\right]-\int_{0}^{t_{f}}dt\left[\frac{\partial^{2}H}{\partial\vec{n}\partial\vec{n}}:\delta\vec{n}\delta\vec{n}+2\frac{\partial^{2}H}{\partial\vec{\nu}\partial\vec{n}}:\delta\vec{\nu}\delta\vec{n}+\frac{\partial
H}{\partial\vec{\nu}\partial\vec{\nu}}:\delta\vec{\nu}\delta\vec{\nu}\right]$
(69)

where in the second step we evaluated this on the instanton. We can write the
bulk part of this as a Hessian operator

$\displaystyle\mathscr{H}=\begin{bmatrix}-\partial_{\vec{n}}\partial_{\vec{n}}H&-\partial_{t}\hat{\delta}-\partial_{\vec{n}}\partial_{\vec{\nu}}H\\\
\partial_{t}\hat{\delta}-\partial_{\vec{\nu}}\partial_{\vec{n}}H&-\partial_{\vec{\nu}}\partial_{\vec{\nu}}H\end{bmatrix},$
(70)

acting in the space $(\vec{n},\vec{\nu})$, where
$(\hat{\delta})_{ij}=\delta_{ij}$.

##  Appendix C Langevin approximation

To develop the Langevin level of approximation, we need to expand $S$ to
$\mathscr{O}(\nu^{2})$. We have

$\displaystyle S[n,\nu;\overline{n},\vec{z}]$
$\displaystyle=\sum_{j}\left[n_{j}(0)\nu_{j}(0)+n_{j}(t_{f})[1-z_{j}+\nu_{j}(t_{f})(z_{j}-1)-\mbox{$\frac{1}{2}$}z_{j}\nu_{j}(t_{f})^{2}]-n^{0}_{j}\nu_{j}(0)\right]$
$\displaystyle\quad+\int_{0}^{t_{f}}dt\left[\vec{\nu}\cdot\partial_{t}\vec{n}-\vec{\nu}\cdot\partial_{\vec{\nu}}H(\vec{n},0)-\mbox{$\frac{1}{2}$}\sum_{k,l}\nu_{k}\nu_{l}\partial_{\nu_{l}}\partial_{\nu_{k}}H(\vec{n},0)\right]+\mathscr{O}(\nu^{3})$
(71)

Now define $B_{ij}(t)=\partial_{\nu_{i}}\partial_{\nu_{j}}H(\vec{n}(t),0)$ and
note that it is positive semi-definite: for any $\vec{x}$,

$\displaystyle\vec{x}\cdot
B\cdot\vec{x}=\sum_{\alpha}F_{\alpha}(\vec{n})(\vec{x}\cdot(\vec{q}_{\alpha}-\vec{p}_{\alpha}))^{2}\geq
0$ (72)

so that it has a unique positive semi-definite square root. Now perform a
Hubbard-Stratonovich transformation

$\displaystyle
e^{\frac{1}{2}\int_{t}\sum_{i,j}B_{ij}\nu_{i}\nu_{j}}\propto\int{\mathscr{D}}\xi\;e^{-\frac{1}{2}\int_{t}\vec{\xi\;}^{2}}e^{\int_{t}\vec{\xi\;}\cdot
B^{1/2}\cdot\vec{\nu}}$ (73)

and similarly

$\displaystyle e^{\frac{1}{2}z_{i}n_{i}(t_{f})\nu_{i}(t_{f})^{2}}\propto
z_{i}^{1/2}\int d\xi^{1}_{i}\
\;e^{-\frac{1}{2}z_{i}(\xi^{1}_{i})^{2}}e^{(n_{i}(t_{f}))^{1/2}z_{i}\xi^{1}_{i}\nu_{i}(t_{f})}$
(74)

so that the path integral is

$\displaystyle
Z=\int{\mathscr{D}}n\int{\mathscr{D}}\nu\int{\mathscr{D}}\xi\;e^{-S_{BC}^{0}}e^{-\int_{0}^{t_{f}}dt\left[\vec{\nu}\cdot\partial_{t}\vec{n}-\vec{\nu}\cdot\partial_{\vec{\nu}}H(\vec{n},0)\right]}e^{-\frac{1}{2}\int_{t}\vec{\xi\;}^{2}}e^{\int_{t}\vec{\xi\;}\cdot
B^{1/2}\cdot\vec{\nu}}.$ (75)

with

$\displaystyle S_{BC}^{0}$
$\displaystyle=\sum_{j}\left[\nu_{j}(0)\big{[}n_{j}(0)-n^{0}_{j}\big{]}+\nu_{j}(t_{f})\big{[}n_{j}(t_{f})(z_{j}-1)-(n_{j}(t_{f}))^{1/2}z_{j}\xi^{1}_{j}\big{]}\right.$
$\displaystyle\qquad\left.-\mbox{$\frac{1}{2}$}\log
z_{j}+n_{j}(t_{f})[1-z_{j}]+\frac{1}{2}z_{j}(\xi^{1}_{j})^{2}\right]$ (76)

The path integral has an action linear in $\nu$, so it can be integrated out
to obtain the Langevin equation

$\displaystyle
0=\partial_{t}\vec{n}-\partial_{\vec{\nu}}H(\vec{n},0)-B^{1/2}\cdot\vec{\xi\;}$
(77)

with boundary conditions

$\displaystyle 0=n_{j}(0)-n^{0}_{j},\qquad
0=n_{j}(t_{f})(z_{j}-1)-\sqrt{n_{j}(t_{f})}z_{j}\xi^{1}_{j}.$ (78)

Since $\xi^{1}\sim\mathscr{O}(1)$, from the second boundary condition we infer
that $n_{j}(t_{f})\sim z_{j}^{2}/(z_{j}-1)^{2}$. Moreover, since
$\vec{\xi\;}^{1}$ has no dynamics, this boundary condition can immediately be
applied to eliminate the variable, to obtain

$\displaystyle S_{BC}^{0}$
$\displaystyle\to\sum_{j}\left[n_{j}(t_{f})[1-z_{j}]+\mbox{$\frac{1}{2}$}\frac{n_{j}(t_{f})(z_{j}-1)^{2}}{z_{j}}-\mbox{$\frac{1}{2}$}\log
z_{j}+\log z_{j}+\mbox{$\frac{1}{2}$}\log(n_{j}(t_{f}))\right]$
$\displaystyle=\mbox{$\frac{1}{2}$}\sum_{j}n_{j}(t_{f})\frac{(1-z_{j}^{2})}{z_{j}}+\mbox{$\frac{1}{2}$}\sum_{j}\log(z_{j}n_{j}(t_{f})).$
(79)

The logarithmic terms come from Jacobians and will be ignored in what follows,
since they are of relative importance $\mathscr{O}(\log n/n)$. The noise has
correlations

$\displaystyle\langle\xi_{j}(t)\rangle$
$\displaystyle=0,\quad\langle\xi_{j}(t)\xi_{k}(t^{\prime})\rangle=\delta(t-t^{\prime})\delta_{jk}$
(80)

###  C.1 Alternative boundary conditions.

Suppose that instead of considering a generating function with $z_{j}$, we
instead ask for the probability that $n_{j}(t)=m_{j}$. Consider this before
introducing $\xi^{1}_{j}$. The $z_{j}$ dependence in $Z(\vec{z},t)$ is then
fully contained in a contribution $-z_{j}n_{j}(t_{f})e^{-\nu_{j}(t_{f})}$ in
the action. We need

$\displaystyle\oint\frac{dz_{j}}{2\pi
i}\frac{1}{z_{j}^{m_{j}+1}}e^{z_{j}n_{j}(t_{f})e^{-\nu_{j}(t_{f})}}$
$\displaystyle=\frac{1}{m_{j}!}(n_{j}(t_{f})e^{-\nu_{j}(t_{f})})^{m_{j}}$
$\displaystyle=\frac{1}{m_{j}!}n_{j}(t_{f})^{m_{j}}e^{-m_{j}\nu_{j}(t_{f})}$
(81)

There is no longer any need to introduce $\xi^{1}_{j}$ and the final-time
boundary condition gets modified to $0=-n_{j}(t_{f})+m_{j}$, as expected. In
this case we will have

$\displaystyle S_{BC}^{0}$ $\displaystyle\to\sum_{j}[n_{j}(t_{f})-m_{j}\log
n_{j}(t_{f})+\log m_{j}!]$
$\displaystyle=\mbox{$\frac{1}{2}$}\sum_{j}\log(m_{j})+\mathscr{O}(1)$ (82)

where the second line uses the Stirling approximation and the boundary
condition $\vec{n}(t_{f})=\vec{m}$.

A useful variation on this is to ask that at time $t_{f}$ the number of
species $j$ is given by $n_{j}(t_{f})=m_{j}+\sqrt{m_{j}}\xi^{1}$, where
$\langle\xi^{1}\rangle=0,\langle(\xi^{1})^{2}\rangle=1$. Then we simply
replace $m_{j}$ by $m_{j}+\sqrt{m_{j}}\xi^{1}$ in the previous result, and
integrate against $e^{-\frac{1}{2}(\xi^{1})^{2}}$. Then to obtain the modified
boundary condition we have

$\displaystyle\int d\nu_{j}e^{\nu_{j}(t_{f})n_{j}(t_{f})}\int
d\xi^{1}e^{-\frac{1}{2}(\xi^{1})^{2}}\frac{1}{\sqrt{m_{j}}}e^{-(m_{j}+\sqrt{m_{j}}\xi^{1})\nu_{j}(t_{f})}$
$\displaystyle=\int
d\xi^{1}e^{-\frac{1}{2}(\xi^{1})^{2}}\frac{1}{\sqrt{m_{j}}}\delta(n_{j}(t_{f})-m_{j}-\sqrt{m_{j}}\xi^{1})$
$\displaystyle=\frac{1}{m_{j}}e^{-\frac{1}{2}\left(\frac{n_{j}(t_{f})-m_{j}}{\sqrt{m_{j}}}\right)^{2}}$
(83)

The dependence on $n_{j}(t_{f})$ is Gaussian, which can be useful.

###  C.2 Scaling with system size

In some applications it is convenient to work with concentrations rather than
particle numbers. Moreover we noted above that $H$ carries a factor of volume.
Let $H=\Omega H^{\prime}$ and $\vec{n}=\Omega\vec{c}$,
$B^{\prime}=\partial_{\vec{\nu}}\partial_{\vec{\nu}}H^{\prime}$. The Langevin
equation becomes
$0=\partial_{t^{\prime}}\vec{c}-\partial_{\vec{\nu}}H^{\prime}(\Omega\vec{c},0)-\Omega^{-1/2}B^{\prime}{}^{1/2}\cdot\vec{\xi\;}$.
This explicitly shows that the noise appears to be sub-dominant for large
systems.

##  Appendix D Time-reversal symmetry

For a forward-reverse reaction pair in local detailed balance, its
contribution to $H$ takes the form

$\displaystyle H^{DB}_{\alpha}(\vec{n},\vec{\nu})$
$\displaystyle=\tilde{k}^{+}_{\alpha}\left[e^{\sum_{j}\nu_{j}(\vec{q}_{\alpha}-\vec{p}_{\alpha})}-1\right]\prod_{i}n_{i}^{p_{\alpha
i}}+\tilde{k}^{-}_{\alpha}\left[e^{-\sum_{j}\nu_{j}(\vec{q}_{\alpha}-\vec{p}_{\alpha})}-1\right]\prod_{i}n_{i}^{q_{\alpha
i}}$ $\displaystyle=\tilde{k}^{+}_{\alpha}\left[e^{\sum_{j}\nu_{j}q_{\alpha
j}}-e^{\sum_{j}\nu_{j}p_{\alpha
j}}\right]\prod_{i}(n_{i}e^{-\nu_{i}})^{p_{\alpha
i}}+\tilde{k}^{-}_{\alpha}\left[e^{\sum_{j}\nu_{j}p_{\alpha
j}}-e^{\sum_{j}\nu_{j}q_{\alpha
j}}\right]\prod_{i}(n_{i}e^{-\nu_{i}})^{q_{\alpha i}}$
$\displaystyle=k_{0}c^{\circ}\Omega
e^{-G_{A_{\alpha}}}\left[e^{\sum_{j}\nu_{j}q_{\alpha
j}}-e^{\sum_{j}\nu_{j}p_{\alpha
j}}\right]\left[\prod_{i}\left(\frac{n_{i}e^{-\nu_{i}}}{c_{i}^{eq}\Omega}\right)^{p_{\alpha
i}}-\prod_{i}\left(\frac{n_{i}e^{-\nu_{i}}}{c^{eq}_{i}\Omega}\right)^{q_{\alpha
i}}\right]$

Now consider the (nonlocal in time) transformation

$\displaystyle\tilde{\nu}_{i}(t_{f}-t)$
$\displaystyle=-\nu_{i}(t)+\log(n_{i}(t)/(c^{eq}_{i}\Omega))$ (84)
$\displaystyle\tilde{n}_{i}(t_{f}-t)$ $\displaystyle=n_{i}(t)$ (85)

Then
$\tilde{n}_{i}(t_{f}-t)e^{-\tilde{\nu}_{i}(t_{f}-t)}/(c^{eq}_{i}\Omega)=e^{\nu_{i}(t)}$
and $e^{\tilde{\nu}_{i}(t_{f}-t)}=n_{i}(t)e^{-\nu_{i}(t)}/(c^{eq}_{i}\Omega)$
so that

$\displaystyle H^{DB}_{\alpha}(\tilde{\vec{n}},\tilde{\vec{\nu}})|_{t_{f}-t}$
$\displaystyle=k_{0}c^{\circ}\Omega
e^{-G_{A_{\alpha}}}\left[\prod_{i}\left(\frac{n_{i}(t)e^{-\nu_{i}(t)}}{c^{eq}_{i}\Omega}\right)^{q_{\alpha
i}}-\prod_{i}\left(\frac{n_{i}(t)e^{-\nu_{i}(t)}}{c^{eq}_{i}\Omega}\right)^{p_{\alpha
i}}\right]\left[\prod_{i}(e^{\nu_{i}(t)})^{p_{\alpha
i}}-\prod_{i}(e^{\nu_{i}(t)})^{q_{\alpha i}}\right]$
$\displaystyle=H^{DB}_{\alpha}(\vec{n},\vec{\nu})|_{t}$ (86)

Moreover

$\displaystyle\int_{0}^{t_{f}}d\tau\tilde{\nu}_{j}(\tau)\partial_{\tau}\tilde{n}_{j}(\tau)$
$\displaystyle=-\int_{0}^{t_{f}}dt\tilde{\nu}_{j}(t_{f}-t)\partial_{t}\tilde{n}_{j}(t_{f}-t)$
$\displaystyle=-\int_{0}^{t_{f}}dt[-\nu_{j}(t)+\log
n_{j}(t)/(c^{eq}_{j}\Omega)]\partial_{t}n_{j}(t)$
$\displaystyle=+\int_{0}^{t_{f}}dt\nu_{j}(t)\partial_{t}n_{j}(t)-\left.\left[\log(n_{j}(t)/(c^{eq}_{j}\Omega))n_{j}(t)\right]\right|_{0}^{t_{f}}+\int_{0}^{t_{f}}dt\partial_{t}n_{j}(t)$
$\displaystyle=+\int_{0}^{t_{f}}dt\nu_{j}(t)\partial_{t}n_{j}(t)-n_{j}(t)\log\frac{n_{j}(t)}{ec^{eq}_{j}\Omega}+n_{j}(0)\log\frac{n_{j}(0)}{ec^{eq}_{j}\Omega}$

Since both $H$ and the kinetic term are invariant up to boundary terms, the
transformation (84) is indeed a symmetry of the action, when all reactions are
in detailed balance.

##  Appendix E Derivation of DMFT for GLV:

We need to compute the disorder average (noted $\overline{\;\cdot\;}$) of

$\displaystyle e^{\int_{t}\sum_{i<j}[S_{ij}A_{ij}+S_{ji}A_{ji}]}$ (87)

where $A_{ij}(t)=[e^{-\nu_{i}(t)}-1]n_{i}(t)n_{j}(t)$. Assuming a Gaussian
distribution for the $\\{S_{ij}\\}$ the disorder-averaged part of the action
is equal to

$\displaystyle\overline{e^{\int_{t}\sum_{i<j}[S_{ij}A_{ij}+S_{ji}A_{ji}]}}=e^{\frac{s}{N}\int_{t}\sum_{i<j}[A_{ij}+A_{ji}]}e^{\frac{\sigma^{2}}{2N}\sum_{i<j}\int_{t,t^{\prime}}[A_{ij}(t)A_{ij}(t^{\prime})+A_{ji}(t)A_{ji}(t^{\prime})+2\gamma
A_{ij}(t)A_{ji}(t^{\prime})]}$ (88)

We have

$\displaystyle\sum_{i<j}\int_{t,t^{\prime}}[A_{ij}(t)A_{ij}(t^{\prime})+A_{ji}(t)A_{ji}(t^{\prime})]=\int_{t,t^{\prime}}\left[\sum_{i,j}A_{ij}(t)A_{ij}(t^{\prime})-\sum_{i}A_{ii}(t)A_{ii}(t^{\prime})\right]$
$\displaystyle=\int_{t,t^{\prime}}\left[\sum_{i,j}[e^{-\nu_{i}(t)}-1][e^{-\nu_{i}(t^{\prime})}-1]n_{i}(t)n_{j}(t)n_{i}(t^{\prime})n_{j}(t^{\prime})-\sum_{i}A_{ii}(t)A_{ii}(t^{\prime})\right]$
$\displaystyle=\int_{t,t^{\prime}}N^{2}\left[P(t,t^{\prime})C(t,t^{\prime})-Q(t,t^{\prime})C(t,t^{\prime})-Q(t^{\prime},t)C(t,t^{\prime})+C(t,t^{\prime})^{2}-\frac{1}{N^{2}}\sum_{i}A_{ii}(t)A_{ii}(t^{\prime})\right],$

where

$\displaystyle P(t,t^{\prime})$
$\displaystyle=\frac{1}{N}\sum_{i}e^{-\nu_{i}(t)}e^{-\nu_{i}(t^{\prime})}n_{i}(t)n_{i}(t^{\prime})$
(89) $\displaystyle Q(t,t^{\prime})$
$\displaystyle=\frac{1}{N}\sum_{i}e^{-\nu_{i}(t)}n_{i}(t)n_{i}(t^{\prime})$
(90) $\displaystyle C(t,t^{\prime})$
$\displaystyle=\frac{1}{N}\sum_{i}n_{i}(t)n_{i}(t^{\prime}),$ (91)

and

$\displaystyle\sum_{i<j}\int_{t,t^{\prime}}A_{ij}(t)A_{ji}(t^{\prime})$
$\displaystyle=\mbox{$\frac{1}{2}$}\int_{t,t^{\prime}}\left[\sum_{i,j}A_{ij}(t)A_{ji}(t^{\prime})-\sum_{i}A_{ii}(t)A_{ii}(t^{\prime})\right]$
$\displaystyle=\mbox{$\frac{1}{2}$}\int_{t,t^{\prime}}\left[\sum_{i,j}[e^{-\nu_{i}(t)}-1][e^{-\nu_{j}(t^{\prime})}-1]n_{i}(t)n_{j}(t)n_{i}(t^{\prime})n_{j}(t^{\prime})-\sum_{i}A_{ii}(t)A_{ii}(t^{\prime})\right]$
$\displaystyle=\mbox{$\frac{1}{2}$}\int_{t,t^{\prime}}N^{2}\Bigg{[}Q(t,t^{\prime})Q(t^{\prime},t)-Q(t,t^{\prime})C(t,t^{\prime})-Q(t^{\prime},t)C(t,t^{\prime})+C(t,t^{\prime})^{2}$
$\displaystyle-\frac{1}{N^{2}}\sum_{i}A_{ii}(t)A_{ii}(t^{\prime})\Bigg{]},$

and

$\displaystyle\sum_{i<j}[A_{ij}(t)+A_{ji}(t)]=N^{2}[p(t)-m(t)]m(t)-\sum_{i}A_{ii}(t),$

with

$\displaystyle p(t)$
$\displaystyle=\frac{1}{N}\sum_{i}e^{-\nu_{i}(t)}n_{i}(t)$ (92) $\displaystyle
m(t)$ $\displaystyle=\frac{1}{N}\sum_{i}n_{i}(t).$ (93)

The path integral is nontrivial due to interactions between different species,
represented here by nonlinear dependence of the action on the above quantities
$P,Q,C,p$ and $m$. To disentangle this dependence we introduce these named
quantities as new variables in the path integral, with corresponding Lagrange
multipliers. Thus, we use

$\displaystyle 1$
$\displaystyle\propto\int{\mathscr{D}}[\cdots]e^{-N\int_{t,t^{\prime}}[\Lambda_{P}(t,t^{\prime})(P(t,t^{\prime})-\frac{1}{N}\sum_{i}e^{-\nu_{i}(t)}e^{-\nu_{i}(t^{\prime})}n_{i}(t)n_{i}(t^{\prime}))]}$
$\displaystyle\quad
e^{-N\int_{t,t^{\prime}}[\Lambda_{Q}(t,t^{\prime})(Q(t,t^{\prime})-\frac{1}{N}\sum_{i}e^{-\nu_{i}(t)}n_{i}(t)n_{i}(t^{\prime}))+\Lambda_{C}(t,t^{\prime})(C(t,t^{\prime})-\frac{1}{N}\sum_{i}n_{i}(t)n_{i}(t^{\prime}))]}$
$\displaystyle\quad
e^{-N\int_{t}[\lambda_{p}(t)(p(t)-\frac{1}{N}\sum_{i}e^{-\nu_{i}(t)}n_{i}(t))+\lambda_{m}(t)(m(t)-\frac{1}{N}\sum_{i}n_{i}(t))]},$
(94)

where

$\displaystyle\int{\mathscr{D}}[\cdots]=\int{\mathscr{D}}\Lambda_{P}\int{\mathscr{D}}\Lambda_{Q}\int{\mathscr{D}}\Lambda_{C}\int{\mathscr{D}}\lambda_{p}\int{\mathscr{D}}\lambda_{m}\int{\mathscr{D}}P\int{\mathscr{D}}Q\int{\mathscr{D}}C\int{\mathscr{D}}p\int{\mathscr{D}}m.$
(95)

The disorder-averaged generating function is then

$\displaystyle\overline{Z}$
$\displaystyle=\int{\mathscr{D}}n\int{\mathscr{D}}\nu
e^{-S_{BC}}e^{-\int_{t}\sum_{j}\nu_{j}\partial_{t}n_{j}}\overline{e^{\int_{t}H}}$
$\displaystyle=\int{\mathscr{D}}n\int{\mathscr{D}}\nu
e^{-S_{BC}}e^{-\int_{t}\sum_{j}\nu_{j}\partial_{t}n_{j}}e^{\sum_{i}\int_{t}V_{i}(n_{i},\nu_{i})}\overline{e^{\int_{t}\sum_{i\neq
j}S_{ij}A_{ij}}}$
$\displaystyle=\int{\mathscr{D}}[\cdots]e^{-NF}\prod_{j}\left(\int{\mathscr{D}}n_{j}\int{\mathscr{D}}\nu_{j}e^{-S_{j}}\right)$
(96)

where $S_{j},F,$ and $V_{i}$ are written out in the main text. The
introduction of new quantities has allowed us to “diagonalize” the species
dependence, such that it takes a product form, as seen in the final expression
for $\overline{Z}$. This sequence of manipulations is essential to obtain a
path integral that can be evaluated with the saddle-point method.

## References

  * [1] Laughlin R B and Pines D 2000 Proceedings of the National Academy of Sciences 97 28–31
  * [2] Goldenfeld N and Woese C 2011 Annu. Rev. Condens. Matter Phys. 2 375–399
  * [3] Schrödinger E and Penrose R 1992 What is Life?: With Mind and Matter and Autobiographical Sketches Canto (Cambridge University Press)
  * [4] Sompolinsky H and Zippelius A 1982 Physical Review B 25 6860
  * [5] Sompolinsky H, Crisanti A and Sommers H J 1988 Physical review letters 61 259
  * [6] Kadmon J and Sompolinsky H 2015 Physical Review X 5 041030
  * [7] Maimbourg T, Kurchan J and Zamponi F 2016 Physical review letters 116 015902
  * [8] Szamel G 2017 Phys. Rev. Lett. 119(15) 155502 URL https://link.aps.org/doi/10.1103/PhysRevLett.119.155502
  * [9] Agoritsas E, Maimbourg T and Zamponi F 2019 Journal of Physics A: Mathematical and Theoretical 52 334001 URL https://doi.org/10.1088/1751-8121/ab2b68
  * [10] Opper M and Diederich S 1992 Physical review letters 69 1616
  * [11] Pearce M T, Agarwala A and Fisher D S 2020 Proceedings of the National Academy of Sciences 117 14572–14583 URL https://www.pnas.org/doi/abs/10.1073/pnas.1915313117
  * [12] Roy F, Barbier M, Biroli G and Bunin G 2020 PLoS computational biology 16 e1007827
  * [13] Bunin G 2017 Physical Review E 95 042414
  * [14] Biroli G, Bunin G and Cammarota C 2018 New Journal of Physics 20 083051
  * [15] Anderson P W 1983 Proceedings of the National Academy of Sciences 80 3386–3390
  * [16] Craciun G and Feinberg M 2006 SIAM Journal on Applied Mathematics 66 1321–1338
  * [17] Bray W C 1921 Journal of the american chemical society 43 1262–1267
  * [18] Briggs T S and Rauscher W C 1973 Journal of chemical Education 50 496
  * [19] Butler T and Goldenfeld N 2009 Physical Review E 80 030902
  * [20] Shih H Y and Goldenfeld N 2014 Physical Review E 90 050702
  * [21] Xue C and Goldenfeld N 2017 Physical review letters 119 268101
  * [22] Constable G W, Rogers T, McKane A J and Tarnita C E 2016 Proceedings of the National Academy of Sciences 113 E4745–E4754
  * [23] Constable G W and McKane A J 2017 Physical Review E 96 022416
  * [24] Georges A, Kotliar G, Krauth W and Rozenberg M J 1996 Reviews of Modern Physics 68 13
  * [25] Kirkpatrick T R and Thirumalai D 1987 Phys. Rev. B 36(10) 5388–5397 URL http://link.aps.org/doi/10.1103/PhysRevB.36.5388
  * [26] Cugliandolo L and Kurchan J 1993 Physical review letters 71 173
  * [27] De Dominicis C 1978 Physical Review B 18 4913
  * [28] Martin P C, Siggia E D and Rose H A 1973 Phys. Rev. A 8(1) 423–437 URL https://link.aps.org/doi/10.1103/PhysRevA.8.423
  * [29] Mézard M, Parisi G and Virasoro M A 1987 Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications vol 9 (World Scientific Publishing Company)
  * [30] Van Kampen N G 1992 Stochastic processes in physics and chemistry vol 1
  * [31] Doi M 1976 Journal of Physics A: Mathematical and General 9 1465
  * [32] Peliti L 1985 Journal de Physique 46 1469–1483
  * [33] Täuber U C 2014 Critical dynamics: a field theory approach to equilibrium and non-equilibrium scaling behavior (Cambridge University Press)
  * [34] Kamenev A 2002 Keldysh and Doi-Peliti techniques for out-of-equilibrium systems (Springer) pp 313–340
  * [35] Andreanov A, Biroli G, Bouchaud J P and Lefevre A 2006 Physical Review E 74 030101
  * [36] Lefevre A and Biroli G 2007 Journal of Statistical Mechanics: Theory and Experiment 2007 P07024
  * [37] Smith E 2011 Reports on Progress in Physics 74 046601
  * [38] Smith E 2020 Entropy 22 1137
  * [39] De Dominicis C 1976 J. Phys. Colloques 37 C1–247–C1–253 URL https://doi.org/10.1051/jphyscol:1976138
  * [40] Janssen H K 1976 Zeitschrift für Physik B Condensed Matter 23 377–380
  * [41] Elgart V and Kamenev A 2004 Physical Review E 70 041106
  * [42] Coleman S 1988 Aspects of symmetry: selected Erice lectures (Cambridge University Press)
  * [43] Van den Broeck C 2013 Stochastic thermodynamics: A brief introduction vol 184 (IOS Press Amsterdam) pp 155–193
  * [44] Gilbert R G and Smith S C 1990
  * [45] Hong L and Qian H 2020 Journal of Physics A: Mathematical and Theoretical 53 425202
  * [46] Fang X, Kruse K, Lu T and Wang J 2019 Reviews of Modern Physics 91 045004
  * [47] Qian H, Ao P, Tu Y and Wang J 2016 Chemical Physics Letters 665 153–161
  * [48] Wang J, Zhang K, Xu L and Wang E 2011 Proceedings of the National Academy of Sciences 108 8257–8262
  * [49] Parisi G 1988 Statistical field theory (Addison-Wesley)
  * [50] May R M 1972 Nature 238 413–414
  * [51] Barbier M, Arnoldi J F, Bunin G and Loreau M 2018 Proceedings of the National Academy of Sciences 115 2156–2161 URL https://www.pnas.org/doi/abs/10.1073/pnas.1710352115
  * [52] Altieri A, Roy F, Cammarota C and Biroli G 2021 Phys. Rev. Lett. 126(25) 258301 URL https://link.aps.org/doi/10.1103/PhysRevLett.126.258301
  * [53] Garcia Lorenzana G and Altieri A 2022 Phys. Rev. E 105(2) 024307 URL https://link.aps.org/doi/10.1103/PhysRevE.105.024307
  * [54] Poley L, Baron J W and Galla T 2022 arXiv preprint arXiv:2208.01569
  * [55] Durrett R and Levin S 1994 Theoretical population biology 46 363–394
  * [56] Ottino-Löffler B and Kardar M 2020 Physical Review E 102 052106
  * [57] Dobramysl U, Mobilia M, Pleimling M and Täuber U C 2018 Journal of Physics A: Mathematical and Theoretical 51 063001
  * [58] Roy F, Biroli G, Bunin G and Cammarota C 2019 Journal of Physics A: Mathematical and Theoretical 52 484001
  * [59] Dornic I, Chaté H and Munoz M A 2005 Physical review letters 94 100601
  * [60] Ceccato A and Frezzato D 2018 The Journal of chemical physics 148 064114
  * [61] Grima R, Thomas P and Straube A V 2011 The Journal of chemical physics 135 084103
  * [62] Schnoerr D, Sanguinetti G and Grima R 2014 The Journal of chemical physics 141 07B606
  * [63] Grima R 2015 Physical Review E 92 042124
  * [64] Wallace E, Gillespie D, Sanft K and Petzold L 2012 IET systems biology 6 102–115
  * [65] Kamenev A, Meerson B and Shklovskii B 2008 Physical review letters 101 268103
  * [66] Kamenev A and Meerson B 2008 Physical Review E 77 061107
  * [67] Charbonneau P, Kurchan J, Parisi G, Urbani P and Zamponi F 2014 Nature communications 5 1–6
  * [68] Scalliet C, Berthier L and Zamponi F 2019 Nature communications 10 1–10
  * [69] Berthier L, Biroli G, Charbonneau P, Corwin E I, Franz S and Zamponi F 2019 The Journal of Chemical Physics 151 010901 URL https://doi.org/10.1063/1.5097175
  * [70] Vincent E, Hammann J, Ocio M, Bouchaud J P and Cugliandolo L F 1997 Slow dynamics and aging in spin glasses Complex Behaviour of Glassy Systems (Springer) pp 184–219
  * [71] Scalliet C and Berthier L 2019 Physical review letters 122 255502
  * [72] Waddington C H 2014 The strategy of the genes (Routledge)
  * [73] Schlögl F 1972 Zeitschrift für physik 253 147–161
  * [74] Vellela M and Qian H 2009 Journal of The Royal Society Interface 6 925–940
  * [75] Battle C, Broedersz C P, Fakhri N, Geyer V F, Howard J, Schmidt C F and MacKintosh F C 2016 Science 352 604–607
  * [76] Tan T H, Watson G A, Chao Y C, Li J, Gingrich T R, Horowitz J M and Fakhri N 2021 arXiv preprint arXiv:2107.05701
  * [77] Lukacs E 1960 Characteristic functions 5 (C. Griffin)
  * [78] De Giuli E 2019 Journal of Physics A: Mathematical and Theoretical 52 504001
  * [79] Advani M, Bunin G and Mehta P 2018 Journal of Statistical Mechanics: Theory and Experiment 2018 033406
  * [80] MacArthur R 1970 Theoretical population biology 1 1–11
  * [81] Tikhonov M and Monasson R 2017 Physical review letters 118 048103
  * [82] De Martino A, Martelli C and Massucci F A 2009 EPL (Europhysics Letters) 85 38007
  * [83] Täuber U C 2012 Nuclear Physics B-Proceedings Supplements 228 7–34