# Effective one-dimension reduction of multi-compartment complex systems
dynamics

Giorgio Vittorio Visco<EMAIL_ADDRESS>Department of
Physics and Astronomy “Galileo Galilei”, University of Padua, Via F. Marzolo
8, 315126 Padova, Italy  Oriol Artime  Departament de Física de la Matèria
Condensada, Universitat de Barcelona, 08028 Barcelona, Spain  University of
Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona,
08028 Barcelona, Spain  Universitat de les Illes Balears, 07122 Palma, Spain
Johannes Nauta  Department of Physics and Astronomy “Galileo Galilei”,
University of Padua, Via F. Marzolo 8, 315126 Padova, Italy  Tomas Scagliarini
Department of Physics and Astronomy “Galileo Galilei”, University of Padua,
Via F. Marzolo 8, 315126 Padova, Italy  Manlio De Domenico
<EMAIL_ADDRESS>Department of Physics and Astronomy “Galileo
Galilei”, University of Padua, Via F. Marzolo 8, 315126 Padova, Italy  Padua
Center for Network Medicine, University of Padua, Via F. Marzolo 8, 315126
Padova, Italy  Istituto Nazionale di Fisica Nucleare, Sez. Padova, Italy

###### Abstract

A broad class of systems, including ecological, epidemiological, and
sociological ones, are characterized by populations of individuals assigned to
specific categories, e.g., a chemical species, an opinion or an epidemic
state, that are modeled as compartments. Due to interactions and intrinsic
dynamics, individuals are allowed to change category, leading to
concentrations varying over time with complex behavior, typical of reaction-
diffusion systems. While compartmental modeling provides a powerful framework
for studying the dynamics of such populations and describe the spatiotemporal
evolution of a system, it mostly relies on deterministic mean-field
descriptions to deal with systems with many degrees of freedom. Here, we
propose a method to alleviate some of the limitations of compartmental models
by capitalizing on tools originating from quantum physics to systematically
reduce multi-dimensional systems to an effective one-dimensional
representation. Using this reduced system, we are able to not only investigate
the mean-field dynamics and their critical behavior, but we can additionally
study stochastic representations that capture fundamental features of the
system. We demonstrate the validity of our formalism by studying the critical
behavior of models widely adopted to study epidemic, ecological and economic
systems.

##  I Introduction

Complex systems consist of a multitude of interconnected and interacting
entities that form extensive networks, with typical examples ranging from
chemical to ecological and social ones [1, 2, 3, 4, 5, 6, 7, 8, 9]. To analyze
their complex behavior, often characterized by collective dynamics, phase
transitions and other emergent phenomena, compartmental modeling is usually
adopted. In this framework, units of the system are categorized or stratified
into distinct species, states or categories, referred to as compartments, and
the concentration of units in each compartment is studied over time while
taking into account the change of species due to either unit-unit interactions
or spontaneous transitions into other compartments.

Compartmentalization facilitates the study of dynamic behaviors and phase
transitions in some order parameters as a function of control parameters that
directly depend on modeling assumptions and observables. For instance, in the
case of an epidemic outbreak it is possible to study how the steady-state
density of infected individuals depends on the transmission rate of a pathogen
by means of social contacts and the recovery rate from the corresponding
disease [10, 11, 12] (see [13] for an extensive review). One of the most
successful frameworks to study the behavior of such epidemic systems, as well
as of the ones consisting of multiple interconnected compartments, is based on
the study of the next-generation matrix, an indispensable tool for analyzing
the conditions required to favor the emergence of active states where
epidemics dynamics is sustained and infection has the potential to scale at
system level [14, 15].

Despite its successes, compartmental modeling comes with some limitations. For
instance, the aforementioned next-generation matrix formalism, while robust,
offers a limited scope in addressing the multitude of degrees of freedom
present in large-scale interconnected systems. In fact, the presence of this
constraint requires the introduction of simplifying assumptions and
approximations, often in the form of deterministic mean-field equations, to
deal with the complexity inherent in these systems. Nevertheless, these
approximations can be still very powerful: for example, in the case of complex
networks with an heterogeneous degree distribution it is possible to consider
classes of a specific degree and then study the dynamics of the overall system
by means of the dynamics of such classes. This approach leads to heterogeneous
mean-field equations that allows one to model the collective behavior of a
system more accurately than simple homogeneous mean-field equations [16, 17,
18].

In fact, this approach shares commonalities with the study of the broad class
of reaction-diffusion processes to describe chemical reaction networks [19],
where units are chemical species undergoing transitions because of chemical
reactions or spontaneous transformations. This parallel is striking and it is
therefore natural to wonder if, and under which conditions, it is possible to
analyze population dynamics of a variety of complex systems by means of a
formalism that overcomes the inherent limitations of the current ones.

Chemical reaction networks (broadly speaking) have been extensively studied by
means of the Doi-Peliti formalism [20, 21], offering a rich methodology for
analyzing both stochastic and deterministic dynamics beyond pure mean-field
approaches. However, while chemical reaction networks are subject to strict
constraints, such as the mass balance of chemical elements [22, 23], this is
not always the case for the general systems we are interested into in this
work, where compartments can also be of heterogeneous nature. Therefore, we
anticipate that it will be necessary to adapt the Doi-Peliti formalism to deal
with the variety of cases of interest for network and complexity science.
Accordingly, in the following, we will refer in general to chemical reaction
networks (CRNs) to indicate such a broad spectrum of systems, from chemical to
ecological and social ones.

The dynamics of a CRN can be completely described by a master equation and the
Doi-Peliti formalism recalls the second quantization used in quantum field
theory to apply the methods for studying quantum systems for the analysis of
stochastic classical systems. One of these methods is the path integral
representation [24], that is particularly useful to go beyond the
deterministic mean-field limit of a model, providing insights on its
stochastic behavior. In the last years, attempts to extend this formalism to
systems with many degrees of freedom have been proposed [25, 26, 24, 27],
although providing general results is still a challenging endeavour.

Building upon these works, we develop a framework that describes multi-
compartment systems by reducing them to a one-dimensional effective model,
showing how to map a CRN to a Hamiltonian problem. This mapping allows us to
characterize the critical behavior, if any, of the CRN through an elegant
geometrical interpretation. Furthermore, irrespectively of the presence, or
lack thereof, of phase transitions in the stochastic process of interest, our
formalism provides a way to approximate its steady-state with high accuracy
even for systems where existing approaches cannot be readily used. The idea of
building on methods originating from quantum physics to study complex
networks, or viceversa, has proven successful to stretch the boundaries of
statistical physics [28, 29, 30, 31, 32, 33, 34, 35, 36]. Similarly, we
foresee broad applications of our dimensionality reduction approach to study
the critical behaviour and the population dynamics of multi-compartment
processes.

The article is organized as follows. We first provide the theoretical basis of
our work in Sec. II, where we introduce the terminology, discuss the Doi-
Peliti’s path integral formalism for the study of CRNs, and lay out the
mapping from CRNs to Hamiltonian systems. In Sec. III, we present the
geometrical interpretation of phase transitions by means of the analysis of
phase portraits, and detail how to classify them. In Sec. IV, we present the
method for dimensionality reduction to an effective one-dimensional model for
the dynamics. Throughout the article, we accompany the theoretical results
with applications to stochastic models of different nature and increasing
complexity, thus validating our framework and discussing its limitations. We
close the article with our conclusions in Sec. V.

##  II Theoretical background

###  II.1 Stochastic compartmental models

Many natural, biological and socio-technical phenomena can be modelled by
means of stochastic processes. Through a probabilistic representation, one
focuses on a macroscopic quantity $n$ whose dynamics is probabilistically
described by a distribution $p(n,t)$ that depends on the features of the
system. Some examples include the disease spreading, where $n$ represents the
number of infected individuals, or ecological systems where this variable
describes species abundance. From a mathematical point of view, $n$ is the
output of a discrete random variable, whereas $p(n,t)$ is the probability of
finding that random variable with value $n$ at time $t$. Generally, we can
take $n$ as a discrete quantity. We will focus on homogeneous Markovian
systems only, therefore the dynamical behaviour of the systems is governed by
a master equation

$\partial_{t}p(n,t)=\sum_{m}H_{nm}p(m,t),$ (1)

where the sum is evaluated over all possible states, and $H_{nm}$ elements of
a matrix that is defined by the system interactions. This matrix can be
written in terms of the transition rates $\omega(n|m)$, that represent the
probability per unit of time of transitioning from the state $m$ to the state
$n$, that is

$H_{nm}=\omega(n|m)-\delta_{n,m}\sum_{m^{\prime}}\omega(m^{\prime}|n),$ (2)

where $\delta_{nm}$ is the Kronecker delta. This compartmental representation
is very flexible, and the evolution of the system narrows down to keeping
trace of the number of elements that belong to each state-compartment.

Now let $M$ be the number of compartments necessary to fully characterize a
system, which we denote $\\{I_{i}\\}_{i=1,\ldots,M}$. In this case,
$\textit{{n}}=(n_{1},n_{2},\ldots,n_{M})$ is the (multi-dimensional) variable
of interest, and $n_{i}$ is the population number in the $i$-th compartment.
For a generic system, $M$ can be arbitrarily large, and the interaction
network very complex. For the sake of simplicity, in this work we limit
ourselves to the case of well-mixed reaction networks, given by the set of
transitions

$\sum_{i=1}^{M}A_{i}\,I_{i}\rightarrow\sum_{i=1}^{M}B_{i}\,I_{i}.$ (3)

where $\\{A_{i},B_{i}\\}_{i=1,\ldots,M}$ represent the stoichiometric
coefficients [27]. Note that these are mass-action reactions, meaning that the
rates are proportional to the population numbers.

Next we will provide methods to study systems fulfilling these transitions. We
will show that our proposed formalism will simplify the analysis of these
systems by focusing only on the effective dynamics of a single compartment.
Before introducing our formalism, let us first present the path integral
formulation and its relation to classical Hamiltonian systems.

###  II.2 The Doi-Peliti formalism and the path integral representation

A useful tool to handle CRNs is the Doi-Peliti (DP) formalism [20, 21]. This
formalism describes the dynamics of a population using creation and
annihilation operators, denoted by $a^{\dagger}$ and $a$ respectively. Thus,
analogous with second quantization, an element of the system can be created or
destroyed following the prescriptions of the ME. This approach has been
extensively covered in the literature (see, e.g., Ref. [24]). While there are
several ways one can introduce the main components of the formalism, we adhere
to the standard definition and develop the usual path integral representation
[21, 24, 27, 37].

We proceed by noting that we use Dirac’s braket notation to describe the state
of a system, i.e. $\ket{\textit{{n}}}=\ket{n_{1},n_{2},\ldots,n_{M}}$. The
corresponding bra is given by the normalization condition

$\bra{\textbf{{n}}}\ket{\textbf{{n}}^{\prime}}=\prod_{i=1}^{M}n_{i}!\delta_{n_{i}n_{i}^{\prime}}.$
(4)

With this normalization condition, the creating and annihilation operators are
defined by

$\displaystyle
a_{i}^{\dagger}\ket{\textit{{n}}}=\ket{\textit{{n}}+\mathbb{I}_{i}}\,,\qquad
a_{i}\ket{\textit{{n}}}=n_{i}\ket{\textit{{n}}-\mathbb{I}_{i}},$ (5)

where $\mathbb{I}_{i}$ is the canonical vector with the $i$-th component equal
to 1 and all others zero. Thus, within the context of compartmental models,
$a_{i}^{\dagger}$ and $a_{i}$ respectively create and annihilate an element
belonging to compartment $I_{i}$ respectively, and follow the usual bosonic
commutation relation. To arrive at a path integral representation, let us
introduce the generating function

$\ket{G}=\sum_{\textit{{n}}}p(\textit{{n}},t)\ket{\textit{{n}}},$ (6)

and, by using Eq. 4, the probability distribution of attaining a specific
state at time $t$ reads

$\displaystyle
p(\textit{{n}},t)=\frac{\bra{\textit{{n}}}\ket{G}}{\prod_{i}n_{i}!}.$

Thus, $\ket{G}$ fully determines the statistical behaviour of the system. By
employing the ME (Eq. 1), we obtain a “Schrödinger-like” equation within the
Doi-Peliti representation, which reads

$\partial_{t}\ket{G}=\mathcal{H}\ket{G},$ (7)

where

$\mathcal{H}=\sum_{\bm{n},\bm{m}}\frac{1}{\prod_{i}m_{i}!}\ket{\bm{n}}H_{nm}\bra{\bm{m}}$
(8)

is the _Doi-Peliti Hamiltonian_ (DPH). For a CRN, the DPH can be written as a
polynomial function of the creation and annihilation operators (Appendix D).

As a basis for the Fock space, we work with

$\ket{\textit{{n}}}=\prod_{i=1}^{M}x_{i}^{n_{i}},$ (9)

where ${\bf\it x}=(x_{1},\ldots,x_{M})\in\mathbb{R}^{M}$. Hence, using the
normalization condition, we have that

$\displaystyle\bra{\bm{n}}=\prod\left(\partialderivative{{\,}}{x_{i}}\right)^{n_{i}}\biggr{\rvert}_{{\bf\it
x}=0}$ (10)

Under these bases, the path-integral representation for the generating
function is given by [24]

$\ket{G}=\int_{0}^{t)}\mathcal{D}[\textbf{{x}},\textbf{{q}}]\,e^{-\mathcal{S}[\textbf{{x}},\mathrm{i}\mkern
1.0mu\textbf{{q}}]+\sum_{j}{n_{0j}}\log{x_{0j}}},$ (11)

where $\mathrm{i}\mkern 1.0mu$ the imaginary unit, $\textbf{{n}}_{0}$ is the
initial state of the stochastic process,
$\mathcal{D}[\textbf{{x}},\textbf{{q}}]$ a suitable measure [24, 37], and
$\mathcal{S}$ is the action of the system defined as (Appendix A)

$\mathcal{S}[\bm{x},\mathrm{i}\mkern
1.0mu\bm{q}]=-\int_{0}^{t}\mathop{}\\!\mathrm{d}{\tau}\left[\sum_{j}\mathrm{i}\mkern
1.0muq_{j}(\tau)\partial_{\tau}x_{j}(\tau)+\mathcal{H}(\bm{x}(\tau),\mathrm{i}\mkern
1.0mu\bm{q}(\tau))\right]$ (12)

Within the path integral representation, the operators $a_{j}^{\dagger}$ and
$a_{j}$ are replaced by two fields $x_{j}(\tau)$ and $\mathrm{i}\mkern
1.0muq_{j}(\tau)$. They are linked to the random variable through the formal
relation $\mathrm{i}\mkern 1.0muq_{j}x_{j}=n_{j}$, which comes from the
definition of the number operators $\hat{N}_{j}=a_{j}^{\dagger}a_{j}$
($\hat{N}_{j}\ket{\textit{{n}}}=n_{j}\ket{\textit{{n}}}$). Note that handling
Eq. 11 is generally difficult, but in the case of chemical reactions, we will
show that we can an approximation similar to the Wentzel-Kramers-Brillouin
(WKB) in spirit, reducing the study of a stochastic process to a Hamiltonian
problem.

###  II.3  Stationary paths and Hamiltonian systems

The path integral representation is formally exact, but to obtain explicit
results one has to resort to approximations. Recently, it has been shown that
the WKB approximation for the study of stochastic systems is a good candidate
[26, 38, 39, 40]. In Appendix B we prove this approximation can be applied to
CRNs as those in Eq. 3, thus greatly simplifying the path integral
representation. Using this approximation, the generating function can be
approximated by $\ket{G}\approx\exp(-\mathcal{S})$ where $\mathcal{S}$ here is
the stationary action. Thus, the principal contributors inside the integration
of Eq. 6 are the solutions $(\widetilde{{\bf\it x}},\widetilde{{\bf\it q}})$
of the Hamilton-Jacobi equations with Hamiltonian $-\mathcal{H}$;

$\displaystyle\partialderivative{\widetilde{{\bf\it
x}}}{t}=-\partialderivative{\mathcal{H}}{\widetilde{{\bf\it
q}}},\qquad\partialderivative{\widetilde{{\bf\it
q}}}{t}=\partialderivative{\mathcal{H}}{\widetilde{{\bf\it x}}}.$ (13)

While this system of equations might seem complicated, it represents an
alternative description to the ME, and, as we will see, valuable information
about the dynamical properties of the stochastic process can be extracted from
it.

As we work with chemical reactions where the WKB approximation holds, we will
consider only stationary trajectories denoted with x and q. In Appendix C, we
show that the stationary trajectories can be taken as real positive functions
satisfying the relation $x_{j}q_{j}=n_{j}$.

Note that the trajectory $x_{j}=1$ is a solution of Eq. 13. It plays an
important role from a statistical point of view. Indeed, by our definitions,
the average value of $n_{j}$ is given by

$\langle
n_{j}\rangle=\partialderivative{\ket{G}}{x_{j}}\biggr{\rvert}_{{\bf\it
x}=1}=q_{j},$ (14)

where we have used the Hamilton-Jacobi relations. Thus, in the phase space,
the deterministic behaviour (i.e. the mean-field behavior) of the system is
contained on the mean-field manifold ${\bf\it x}=1$. As we will illustrate,
this manifold is key in the analysis of the phase transitions on the phase
portrait of the system.

##  III Phase portraits and phase transitions

So far we have used the Doi-Peliti formalism to map the dynamics of a CRN to a
Hamiltonian system through two vectors that encode the dynamics of the process
x and $\mathrm{i}\mkern 1.0mu\textit{{q}}$. To provide a clearer physical
interpretation of our variables, we use the Cole-Hopf (CH) transformation [41]

$\displaystyle x_{j}=e^{\mathrm{i}\mkern 1.0mu\theta_{j}},\qquad
q_{j}=\eta_{j}e^{-\mathrm{i}\mkern 1.0mu\theta_{j}},$ (15)

with $\eta_{j},\theta_{j}\in\mathbb{R}$. In doing so, we replace the field
$(\bm{n},\mathrm{i}\bm{q})$ with $(\bm{\eta},\mathrm{i}\bm{\theta})$, where
now $\eta_{j}$ represents the population number $n_{j}$. With these variables,
the action $\mathcal{S}$ becomes

$\displaystyle\mathcal{S}$
$\displaystyle=\int_{t_{0}}^{t}\mathop{}\\!\mathrm{d}{\tau}\,\left[\sum_{j}\mathrm{i}\mkern
1.0mu\theta_{j}\partial_{\tau}\eta_{j}-\mathcal{H}\left(e^{\mathrm{i}\mkern
1.0mu\bm{\theta}},\bm{\eta}e^{-\mathrm{i}\bm{\theta}}\right)\right]$
$\displaystyle\qquad\qquad\qquad\qquad+\mathrm{i}\mkern
1.0mu\bm{\theta}(t_{0})\bm{\eta}(t_{0})-\mathrm{i}\mkern
1.0mu\bm{\theta}(t)\bm{\eta}(t),$

and the system of Eq. 13 reads

$\displaystyle\begin{split}\partialderivative{{\bf\it\eta}}{t}&=-\mathrm{i}\mkern
1.0mu\partialderivative{{\bf\it\theta}}\mathcal{H}\left(e^{\mathrm{i}\mkern
1.0mu{\bf\it\theta}},{\bf\it\eta}e^{-\mathrm{i}\mkern
1.0mu{\bf\it\theta}}\right)\\\ \mathrm{i}\mkern
1.0mu\partialderivative{{\bf\it\theta}}{t}&=\;-\partialderivative{{\bf\it\eta}}\mathcal{H}\left(e^{\mathrm{i}\mkern
1.0mu{\bf\it\theta}},{\bf\it\eta}e^{-\mathrm{i}\mkern
1.0mu{\bf\it\theta}}\right)\end{split}$ (16)

Note that mean-field manifold $\bm{x}=1$ ($\bm{\theta}=0$) is always a
solution of both Eqs. 13 and 16. Thus, it is a universal feature of CRNs, as
discussed in [42]. A more formal proof of this property is provided in
Appendix D.

As it has been done for $(\bm{x},\mathrm{i}\mkern 1.0mu\bm{q})$, in the
stationary limit we replace $(\bm{\eta},\mathrm{i}\mkern 1.0mu\bm{\theta})$
with the (real) solutions of (16). With this change of coordinates, the mean-
field dynamics is recovered on $\bm{\theta}=0$. Thus, the physical
interpretation of $\bm{\theta}$ relates to how much the dynamical behaviour of
the system is close to the deterministic one. In other words, $\bm{\theta}$
expresses the role of the stochastic fluctuations _beyond_ the mean-field
approximation. Therefore, by studying the phase portrait of (16), we can
investigate the time evolution through its attractors. Moreover, we are able
to probe the role of fluctuations, which shall be discussed in Section III.2.

Note that while presented transformation offers a clearer interpretation of
the fields representing the creation and the annihilation operators, it
implies technical difficulties as well, as shall be ilustrated below. In the
following sections, we will alternate between the $(\bm{x},\bm{q})$
representation and the Cole-Hopf $({\bf\it\eta},{\bf\it\theta})$ one, showing
that the former tends to be more suitable for analytic computations, while the
latter is more convenient for graphical visualizations of the phase portraits.

###  III.1 The Doi-Peliti formalism and the characterization of phase
transitions

Following ideas put forward in [25], we build on the formalism presented so
far to develop a mathematical framework for the analysis of phase transitions.
In particular, our approach overcomes the limitations present when considering
one-compartment models only, thus offering a tool for the analysis of the
critical and stochastic behavior of general multidimensional compartmental
model. We consider systems that converge to a steady state for $t\to\infty$.
These systems are typically out of equilibrium, hence they cannot be treated
with the tools of equilibrium statistical mechanics, such as the Landau theory
[43]. The non-equilibrium properties, however, can be easily addressed by
leveraging the Doi-Peliti mapping to the Hamiltonian representation.

We will describe our system by means of a non-equilibrium order parameter
$\rho$, which denotes the asymptotic density of a specific compartment. We
focus on cases in which the initial size of a compartment of a certain state
$I_{1}$ is large. That is, $n_{1}\approx N$, and hence

$\rho=1-\lim_{t\to\infty}\frac{\langle n_{1}\rangle}{N}.$ (17)

Within the context of epidemiology, this means that the number of susceptible
individuals dominates, and the study of the dynamics consists in observing how
the system responds by inserting a small amount of infected individuals.

As we are interested in the large time behaviour, Eq. 7 hints towards focusing
our attention on the zero-energy sets $\mathcal{H}=0$. Within the WKB
approximation, we are able to study the topology of these sets to obtain
information about the system at large times. Here, by topology we mean the
structure of the trajectories on the phase portrait related to $\mathcal{H}$.
The fixed points of the Hamiltonian system on the manifold $\mathcal{H}=0$
play an important role since they are related to the stable and the unstable
manifold that determine the system’s evolution [44]. In particular, the
dynamics converges at the hyperbolic fixed point on the stable manifold, and
diverges from it on the unstable one. Therefore, when modifying the control
parameters, such as the coefficients of transition rates of the web of
interactions among compartments, the configuration of the fixed points
changes. As a result, the dynamics of the system will vary with its
attractors, suggesting the presence of phase transitions. Moreover, as stated
above, $\bm{x}=1$ is the stationary trajectory for which $\mathcal{H}=0$, and
it represents the mean-field dynamics. By studying the fixed points on
$\bm{x}=1$, we can recover the mean-field phase transitions, and, by analysing
the topology of the stable and unstable manifold, we are able to characterize
these transitions [25]. If new fixed points appear when the control parameters
are modified, then the dynamics of the system radically change, leading to
first-order phase transition. Instead, if we observe that the number of fixed
points does not vary but their positions change, then the asymptotic
configuration continuously varies with respect to the control parameters and
thus a second-order transition is present.

After these considerations, we define the following strategy for studying
possible phase transitions for a multi-compartment system:

  1. 1.

find the Doi-Peliti Hamiltonian of the system;

  2. 2.

look for the fixed points on the zero-energy sets;

  3. 3.

look for conditions at which the determinant of the system Jacobian is null;

  4. 4.

find parameter relations that characterize the phase transition;

When looking for the zero-energy sets, one must take all possible constraints
into account, such as conservation of the total number of elements. In doing
so, the system will be reduced to a size $2M-2m$, where $m$ is the number of
constraints (see below). The determinant of the Jacobian is null when the
Jacobian has at least two null eigenvalues. And finally, by varying the
(control) parameters near the critical point(s), we can determine whether
there is a first- or second-order phase transition, as we described above.

Regarding the constraints, we represent these by an equation of the form
$\bm{C}=\mathcal{C}(\bm{x},\bm{q})$, where $\bm{C}$ is now an $m$-dimensional
vector. Assuming that the implicit function theorem holds, we carry out a
function $Q$ such that, for $m$ components of $\bm{q}$, we have:

$(q_{1},\ldots,q_{m})=Q(q_{m+1},\ldots,q_{M},x_{1},\ldots,x_{M},\bm{C}).$ (18)

Furthermore, constraints represent first integrals of the Hamiltonian system.
By substituting Eq. 18 into the DPH, we obtain a Hamiltonian that is
independent of these variables, therefore the time derivatives of the
corresponding conjugate coordinates $(x_{1},\ldots,x_{m})$ are all zero, and
thus we can let $x_{i}=1$ for all $i=1,\ldots,m$. Each constraint reduces the
dimension of the system by $2$, projecting the DPH onto a new, reduced
Hamiltonian $\mathcal{H}^{\star}$ that depends on $2M-2m$ variables. Under
these premises, the system reads

$\displaystyle\partial_{x_{j}}\mathcal{H}^{\star}=0,\qquad\partial_{q_{j}}\mathcal{H}^{\star}=0$
(19)

with the constraints

$\displaystyle x_{j}=1,\qquad\det(\mathbf{J})=0,$ (20)

where J is the Jacobian of $\mathcal{H}^{\star}$ in terms of $\bm{x}$ and
$\bm{q}$. By treating the coefficients of the transition rates as variables,
the system of Eq. 19 gives us the critical points for possible phase
transitions. In the same way, we can build this system by using the Cole-Hopf
coordinates.

In the following section we shall use the Cole-Hopf variables and show how to
derive the Langevin dynamics starting from a generic DPH, and the accompanying
Fokker-Planck equation for the corresponding probability distribution.

###  III.2 Fokker-Planck and Langevin representations

The Cole-Hopf variables allow us to represent the fluctuations around the
mean-field dynamics by means of $\bm{\theta}$. We have seen that, from the
first derivative of $\mathcal{H}^{\star}$ with respect to $\bm{\theta}$
evaluated at zero, we get the dynamics of the reduced system. Now, we show
that, from its second derivative, we can infer the behavior of the system due
to fluctuations, as done in [27, 24].

Let us consider, for now, a single-compartment system and take the reduced
Hamiltonian $\mathcal{H}^{\star}$ in the Cole-Hopf representation. When the
population number $n$ is sufficiently large, it can be approximated by a
continuous variable. In doing so, we can re-write the ME of Eq. 1 using a
Kramers-Moyal expansion [24],

$\displaystyle\partial_{t}p(n,t)$ $\displaystyle=\sum_{m\geq
1}\frac{(-1)^{m}}{m!}\partialderivative{{}^{m}}{n}\big{[}W_{m}(n)p(n,t)\big{]}$
$\displaystyle=\mathcal{H}_{KM}(n,\partial_{n})p(n,t),$ (21)

with

$W_{m}(n)=\int_{\mathbb{R}}\mathop{}\\!\mathrm{d}{y}\,y^{m}\omega(n+y,n),$
(22)

and where $\mathcal{H}_{KM}$ the Kramers-Moyal operator. In this limit, a link
between the Doi-Peliti Hamiltonian and $\mathcal{H}_{KM}$ exists. Using the
definition of $\mathcal{H}$, and considering Eq. (8), one can prove that for
one-dimensional systems the following holds;

$\mathcal{H}^{\dagger}_{KM}(\eta,\partial_{\eta})=\mathcal{H}(\eta,\theta)$
(23)

where $\mathcal{H}^{\dagger}_{KM}$ is the adjoint operator of
$\mathcal{H}_{KM}$. Under this mapping, we have $\eta=n$ and
$\theta=\partial_{n}$.

The same reasoning can be applied to systems with more degrees of freedom by
using the dimension reduction presented above. If we consider the average of
$n_{1}=\eta$ as the order parameter, by expanding $\mathcal{H}^{\star}$ up to
the second-order of $\theta$ around $\theta=0$, we then obtain a Fokker-Planck
equation [45]

$\displaystyle\partial_{t}p(\eta,t)$
$\displaystyle=-\partial_{\eta}\big{[}\partial_{\theta}\mathcal{H}^{\star}_{0}p(\eta,t)\big{]}+\partial^{2}_{\eta}\big{[}\tfrac{1}{2}\partial^{2}_{\theta}\mathcal{H}^{\star}_{0}p(\eta,t)\big{]},$
(24)

where, with some slight abuse of notation, we have
$\partial^{k}_{\theta}\mathcal{H}^{\star}_{0}=(\partial^{k}_{\theta}\mathcal{H}^{\star})|_{\theta=0}$
and $p(\eta,t)$ the conditional probability distribution constrained by the
stationary conditions. Hence, by expanding the reduced Hamiltonian, we can
write an effective one-dimensional Fokker-Planck equation for the multi-
compartment stochastic process, that makes further calculations considerably
easier. From the Fokker-Planck equation, we can readily obtain the Langevin
equation of the problem [45]

$\displaystyle\derivative{\eta}{t}$
$\displaystyle=\partial_{\theta}\mathcal{H}^{\star}_{0}+\sqrt{\partial^{2}_{\theta}\mathcal{H}^{\star}_{0}}\xi(t)$
(25)

with $\xi(t)$ Gaussian white noise. Formal steady-state solutions of Eq. 23
are known to be [46]

$\lim_{t\to\infty}p(\eta,t)=\mathcal{Z}\,e^{-V(\eta)},$ (26)

with the potential

$\displaystyle
V(\eta)=\int\mathop{}\\!\mathrm{d}{\eta}\;\frac{\partial_{\eta}(\partial^{2}_{\theta}\mathcal{H}^{\star}_{0})-2\partial_{\theta}\mathcal{H}^{\star}_{0}}{\partial^{2}_{\theta}\mathcal{H}^{\star}_{0}},$
(27)

and $\mathcal{Z}$ the normalization factor. Here it is clear that that
$\partial_{\theta}^{2}\mathcal{H}^{\star}$ is related to the system’s
fluctuations and that, by employing the Doi-Peliti formalism, we can gather
together both the deterministic and the stochastic dynamics.

Eqs. 19 and 25 constitute the grounds on which we will develop a framework
adapted for the systematic study of multi-dimensional systems. Yet, we shall
first show how to combine the results reported so far to handle models with
two compartments. In Section III.3, we treat a modified $SIS$ model in which
we inserted higher order interactions. We observe how the Doi-Peliti
representation allows us to gain information on the critical behavior by
unveiling tricritical points. In Sec. III.4, we analyse phase transitions
depending on the system initial configuration. Afterwards, we will then extend
this approach to multi-dimensional problems by presenting a novel dimensional
reduction method (Sec. IV.1), which offer insights even in cases where other
methods fail, as the well-known next-generation matrix [15]. Similarly, for
illustrative purposes, we will showcase the implementation of this in models
with interesting critical behavior: a modified $SIRS$ model (Sec. IV.2), two
Lotka-Volterra models (Secs. IV.3 and IV.3.2) and an exploiter-resource model
IV.4 proposed by us. In addition, in Appendix F we report further examples
applied to several other contexts, such as tax evasion models. In Table 1 we
summarise all the models treated in this paper with the main results.

Table 1: Models studied with the Doi-Peliti representation, in order of
appearance in the article. For each model, we report the number of compartment
and constraints, the compartment related to the order parameter, with its
average taken at the stationary limit $t\to\infty$, and the type of phase
transaction identified with our approach. The Modified $SIS$ models are
treated in Sec. III.3, the Modified $SIRS$ model in Sec. IV.2, the Lotka-
Volterra model in Sec. IV.3 and the Resource-Exploiters model in Sec. IV.4. In
the Appendix, we address the $SEIRS$ model in Sec. F.1, the $SEIR$ model with
latent categories $E_{1,\,2}$ in Sec. F.2, the Tax evasion model in Sec. F.3
and the Levin model in Sec. F.4.

###  III.3 Application to epidemiology: A modified $SIS$ model

Let us now analyse a modified Susceptible-Infected-Susceptible (SIS) model to
illustrate our method. The SIS model is an epidemiological model that
describes the interactions between elements of a population that are in either
of two compartments; the susceptible $S$ and the infected $I$ compartments.
Upon contact with an infected individual, a susceptible can can become
infected with rate $\beta$. Infected individuals recover from the disease with
rate $\gamma$. The elementary reactions thus read

$\displaystyle S+I$ $\displaystyle\xrightarrow{\beta/N}2I,$ (28)
$\displaystyle I$ $\displaystyle\xrightarrow{\gamma}S$

This system has been studied extensively and it is known that there is a
continuous phase transition at $\beta=\gamma$ [13]. For $\beta>\gamma$, in the
limit $t\to\infty$ the disease prevails in the population, impacting a
macroscopic fraction of individuals, and so $\rho=1-\langle n_{I}\rangle/N>0$.
For $\beta<\gamma$, the disease will eventually disappear, and the system
converges to the inactive (absorbing) phase. In general, for such
epidemiological models, the initial condition is taken $n_{S}\approx N$, but
obviously $n_{S}\leq N$.

To showcase the potential of this framework with respect to more traditional
approaches, we include an additional reaction that underlies richer critical
scenarios. More specifically, we consider an additive elementary reaction,

$\displaystyle S+2I\xrightarrow{\alpha/N^{2}}3I,$ (29)

that corresponds to an additional transition rate,

$\omega(\bm{n}-\mathbb{I}_{S}+\mathbb{I}_{I},\bm{n})=\alpha\frac{n_{S}n_{I}(n_{I}-1)}{N^{2}}.$
(30)

This is a three-body-like interaction that may interpreted as a process that
helps the diffusion of the disease, mediated by the parameter $\alpha$.

Taken all the rules together, the Doi-Peliti Hamiltonian reads

$\displaystyle\begin{split}\mathcal{H}&=\frac{\alpha}{N^{2}}(a_{I}^{\dagger})^{2}\left(a^{\dagger}_{I}-a_{S}^{\dagger}\right)a_{S}a_{I}^{2}\\\
&\quad+\frac{\beta}{N}(a_{I}^{\dagger}-a_{S}^{\dagger})a_{I}^{\dagger}a_{I}a_{S}\\\
&\quad+\gamma\left(a^{\dagger}_{S}-a^{\dagger}_{I}\right)a_{I},\end{split}$
(31)

and the corresponding mean-field trajectories read

$\displaystyle\partialderivative{\langle n_{S}\rangle}{t}$
$\displaystyle=\gamma\langle n_{I}\rangle-\frac{\beta}{N}\langle
n_{S}\rangle\langle n_{I}\rangle-\frac{\alpha}{N^{2}}\langle
n_{S}\rangle\langle n_{I}\rangle^{2},$ (32)
$\displaystyle\partialderivative{\langle n_{I}\rangle}{t}$
$\displaystyle=-\gamma\langle n_{I}\rangle+\frac{\beta}{N}\langle
n_{S}\rangle\langle n_{I}\rangle+\frac{\alpha}{N^{2}}\langle
n_{S}\rangle\langle n_{I}\rangle^{2}.$ (33)

Here, we the total population is conserved, i.e. $N=n_{S}+n_{I}$. Therefore,
by focusing on the number of susceptible individuals, we define the
constraints $C=N$ and

$q_{I}=Q(q_{S},x_{S},x_{I},N)=N-x_{I}q_{I}.$ (34)

After applying the Cole-Hopf transformation the reduced DPH reads

$\mathcal{H}^{*}=(1-e^{\theta})(N-\eta)\left(\frac{\alpha}{N^{2}}\eta
e^{-\theta}(N-\eta)+\frac{\beta}{N}\eta e^{-\theta}+\gamma\right),$ (35)

where $\eta=\langle n_{S}\rangle$. The zero-energy lines are given by

$\displaystyle\begin{split}\theta=0,&\qquad\eta=N,\\\
\Theta(\eta,\beta,\gamma,\alpha)&=\log\left(-\frac{\eta(\eta\alpha-N(\alpha+\beta))}{N^{2}\gamma}\right).\end{split}$
(36)

Within the mean-field approximation, we consider only the fixed points
$F=(\eta,\theta)$ with $\theta=0$, that is,

$\displaystyle F_{1}$ $\displaystyle=(N,0),$ $\displaystyle F_{\pm}$
$\displaystyle=\left(N\frac{\alpha+\beta\pm\sqrt{\alpha^{2}+2\alpha\beta-4\alpha\gamma+\beta^{2}}}{2\alpha},0\right).$

By solving the system of Eq. 19 and imposing $\alpha,\,\beta,\,\gamma\geq 0$,
we find two different critical conditions,

$\displaystyle\beta=\gamma,\qquad\beta=-\alpha+2\sqrt{\alpha\gamma},$ (37)

and we thus expect to find two distinct phase transitions. We can distinguish
these two critical behaviors, as shown in Fig. 1, where one transition occurs
for $\alpha<\gamma$ and another for $\alpha>\gamma$. First, we note that the
condition $\beta=\gamma$ is a critical point for any $\alpha$. Indeed, as can
be seen in Figs. 1(a) and 1(d), for $\gamma>\beta$, $F_{1}$ is the only
attractor of the systems (provided that $\langle n_{S}\rangle\leq N$). When
$\beta=\gamma$, the attractor changes. This can also be observed when
considering a deterministic approach by means of linearizing the mean-field
equations about the fixed point. When doing so, the dynamics around $F_{1}$
when $\theta=0$ reads,

$\displaystyle\partial_{t}\eta$
$\displaystyle=(\eta-N)\left(\partial^{2}_{\eta,\theta}\mathcal{H}^{\star}\right)\big{\rvert}_{F_{1}}=(\beta-\gamma)(\eta-N).$
(38)

Hence, the asymptotic state varies at $\beta=\gamma$. To fully characterize
the critical behavior, the remaining critical condition (Eq. 37) needs to be
studied as well. We note that, at $\beta\geq-\alpha+2\sqrt{\alpha\gamma}$, two
new fixed points $F_{\pm}$ appear, suggesting the presence of an abrupt
transition. When $\beta=-\alpha+2\sqrt{\alpha\gamma}$, a non-trivial zero-
energy line intersects tangentially the mean-field manifold $\theta=0$ at
$F_{+}=F_{-}$. For $\gamma>\alpha$ the tangent point is above the absorbing
state $F_{1}$, otherwise it stays below it. Fixing $\alpha$ and $\gamma$, with
$\gamma>\alpha$, and varying $\beta$, we observe that $F_{-}=F_{1}$ at
$\beta=\gamma$ (Fig. 1(b)). By further increasing $\beta$, the intersection
$F_{-}$ smoothly moves below the absorbing state becoming the new asymptotic
state (Fig. 1(c)), resulting in a second-order phase transition. On the other
hand, for $\gamma<\alpha$ we have $F_{+}=F_{1}$ at the critical condition
(Fig. 1(e)), and, for $\beta>\gamma$, $F_{+}$ is above the absorbing state
(Fig. 1(f)). In this case, $F_{-}$ becomes attractive, so the dynamics jumps
to $F_{-}$ and the transition is of first-order.

To summarize, our system is characterized by the following critical behavior:

$\displaystyle\beta=\gamma$ $\displaystyle\qquad\textrm{critical condition},$
$\displaystyle\beta=\alpha=\gamma$ $\displaystyle\qquad\textrm{tricritcal
point},$ $\displaystyle\alpha>\gamma$ $\displaystyle\qquad\textrm{first-order
phase transition},$ $\displaystyle\alpha<\gamma$
$\displaystyle\qquad\textrm{second-order phase transition}.$

This analysis exemplifies how one should approach the study of critical
behavior of CRNs within the context of the Doi-Peliti formalism.

(a) $\beta<\gamma,\,\alpha<\gamma$

(b) $\beta=\gamma,\,\alpha<\gamma$

(c) $\beta>\gamma,\,\alpha<\gamma$

(d) $\beta<\gamma,\,\alpha>\gamma$

(e) $\beta=\gamma,\,\alpha>\gamma$

(f) $\beta<\gamma,\,\alpha>\gamma$

Figure 1: Phase portrait for the modified $SIS$ model by using the Cole-Hopf
variables. The orange lines represent the mean-field manifold, while the blue
ones are the other non zero-energy lines. The fixed points are $F_{1}$ and
$F_{\pm}$, Red dots are the attractors of the deterministic dynamics. In
Panels 1(a),1(c) and 1(b) we show a second-order phase transition: at the
critical condition $\beta=\gamma$, $F_{-}$ overlaps with $F_{1}$ becoming the
new attractor. In the lower Panels 1(d),1(e) and 1(f), the first-order
transition is reported: at $\beta=\gamma$, $F_{1}=F_{+}$ and the attractor
changes abruptly.

###  III.4 The role of the initial conditions

In general, CRNs can exhibit highly non-trivial phase diagrams with many
different attractors — especially when the number of compartments grows. Under
certain conditions, this may impact the procedure introduced in this work.
Indeed, this methodology does not provide the critical point certainly, but
indicates conditions in which variations of the topology of the phase space
lead to possible changes in the asymptotic dynamical behaviour of the system.
Since we are analysing systems that are out of equilibrium, the steady states
can depend on the initial conditions, and the dynamical evolution of the
system may be independent of the topological configuration variations
identified by the formalism. In these cases, the possible phase transitions
depend on the initial state and, in general, it is impossible to identify an
analytical condition which provides the critical point and characterizes the
phase transition.

(a) $\alpha<4\gamma$

(b) $\alpha=4\gamma$

(c) $\tilde{\alpha}=\alpha>4\gamma$

(d) $\alpha>\tilde{\alpha}$

(e)

Figure 2: Characterization of the critical behavior of the two-compartment
system with a three-body interaction (modified $SIS$ model of Sec. III.3).
From panels (a) to (d), the orange lines represent the mean-filed manifold,
while the blue one are the other zero-energy lines. Red dots are the
attractors of the deterministic dynamics. $F_{0}$ represent the initial
condition. In (e), we display the order parameter as a function of the
parameter controlling the strength three-body interaction, obtained from
simulations. For the numerical experiments, we have considered $\gamma=0.2$
and we have normalized the population setting $N=1$. We observe that for
$n_{S}(0)>1/2$ the phase transition occurs for $\alpha>4\gamma$. Here
$\tilde{\alpha}_{1}=0.953$ and $\tilde{\alpha}_{2}=1.25$ as predicted by Eq.
(40). Moreover, for all initial conditions, the active phase is described by
the same curve: this curve is given by the value of $F_{-}$ as a function of
$\alpha$.

To illustrate this point, let us consider again the modified $SIS$ model, but
now with $\beta=0$. This leaves only the three-body interaction as the sole
infection mechanism. By repeating the analysis, we find for this configuration
three fixed points along the mean-field line, which are

$\displaystyle F_{1}=(N,0)\,,\quad
F_{\pm}=\left(N\left(\frac{1}{2}\pm\sqrt{\frac{\alpha-4\gamma}{4\alpha}}\right),0\right).$

After rescaling, i.e. $N=1$, and solving Eq. 19, the critical condition reads

$F_{-}=F_{+}=\left(\tfrac{1}{2},0\right)\,,\qquad\alpha=4\gamma$ (39)

At this condition, new fixed points appear (Fig. 2). According to our
expectations, there is a possible first-order transition. However, when
assuming $n_{S}\approx 1$, a phase transition is not observed when numerically
integrating the dynamical system. When varying the rates, we can determine the
attractive or repulsive character of the fixed points. For $\alpha<4\gamma$,
$F_{1}$ is the only attractor (Fig. 2(a)). For $\alpha>4\gamma$, if $n_{S}(0)$
is greater than the fixed value $\tfrac{1}{2}$, the attractor does not vary,
and we can conclude there is no phase transition. In fact, this occurs only
when there is an overlap between a hyperbolic fixed point and the initial
configuration of the system as shown in Figs. 2(c) and 2(d). This condition is
given by a rate $\widetilde{\alpha}$ that can found by solving
$F_{+}=(n_{S}(0),0)$,

$\tilde{\alpha}(n_{S}(0))=\frac{\gamma}{(1-n_{S}(0))n_{S}(0)}.$ (40)

In the case $n_{S}(0)\approx 1$, we have $\widetilde{\alpha}\to\infty$, so the
system does not display an effective phase transition. On the other hand, if
$n_{S}(0)<\tfrac{1}{2}$, then the variation of the eigenspace does affect the
system and we have an effective phase transition described by the condition of
Eq. 39. The existence of this phase transition can be verified by numerical
integration for different initial conditions (Fig. 2(e)).

The very same arguments can be extended to the case with several compartments.
However, the solution becomes more complicated because the inherent difficulty
to systematically identify which eigenspace a generic initial configuration
belongs to. In this case, an alternative solution might be to employ numerical
investigations to verify whether the phase transition takes place or not. We,
however, provide yet another alternative approach that is related to reducing
the dimensions of the system, which we shall now expand on.

##  IV Dimension reduction

###  IV.1 Theoretical results

Dynamical models with many components are common and their characterization as
Hamiltonian systems can be complicated. Moreover, only in the case of single-
compartment systems we are able to find a graphical visualization of the level
sets with the tools presented so far. To find a strategy to extend the Doi-
Peliti representation to multi-compartment models, we start by noting that the
presence of first integrals helps simplify the system. We employ this property
to build our own dimension reduction method, which is based on using the
deterministic equations as first integrals to freeze the degrees of freedom
that are not inherent in the order parameter. From this, we are able to obtain
an effective Hamiltonian $\mathcal{H}^{\star}$ that captures the fundamental
features of the system.

Let us consider a general DPH for a chemical reaction process. Take the
deterministic trajectories

$\displaystyle\begin{split}\partial_{t}\langle n_{1}\rangle&=f_{1}(\langle
n_{1}\rangle,\langle n_{2}\rangle,\ldots),\\\ \partial_{t}\langle
n_{2}\rangle&=f_{2}(\langle n_{1}\rangle,\langle n_{2}\rangle,\ldots),\\\
\vdots&\\\ \partial_{t}\langle n_{d}\rangle&=f_{M}(\langle
n_{1}\rangle,\langle n_{2}\rangle,\ldots),\end{split}$ (41)

where $f_{i}$ are functions of occupation numbers obtained by the solutions of
the Hamiltonian problem setting $x_{j}=1$ and $q_{j}=\langle n_{j}\rangle$.
Since we are assuming that an asymptotically steady state is reached when
$t\to\infty$, this set of equations admits a solution for $\partial_{t}\langle
n_{j}\rangle=0$ for all $j$. As before, we take $I_{1}$ as the starting
compartment with $\langle n_{1}\rangle\approx N$ at $t=0$. For large $t$, we
rewrite the previous system as

$\displaystyle\begin{split}\partial_{t}\langle n_{1}\rangle&=f_{1}(\langle
n_{1}\rangle,\langle n_{2}\rangle,\ldots),\\\ \partial_{t}\langle
n_{j}\rangle&=0,\end{split}$ (42)

meaning that $f_{i}(\langle n_{1}\rangle,\langle n_{2}\rangle,\ldots)=0$ for
all $j>1$ — i.e. for all compartments but the first. In doing so, we are
imposing constraints on the average values of the occupation numbers, from
which we derive the set of first integrals. Fixing these constraints is
equivalent to seeing how all degrees of freedom converge at equilibrium points
as a function of a single compartment. We emphasise that the constraints of
Eq. 42 regards only the average values of the mean-field approximation. Making
an analogy with statistical mechanics, we are setting a constraint on the
expectation value of a quantity, similarly to the canonical ensemble
constraining the average value of the system energy. Therefore, the physical
interpretation of this dimension reduction is as follows: In the large time
limit, the averages of occupation numbers converge to a hyperbolic fixed point
in an $M$-dimensional space following the constraints imposed by Eqs. (42).

Further physical insights can be gained by studying what occurs near this
limit. Let $x_{1}$ and $q_{1}$ be the only variables related to the order
parameter, while all the other degrees of freedom are frozen. From an analytic
point of view, by using Eqs. (42) we can rewrite $q_{j}$ for all $j>1$ as
functions of $q_{1}$ and $x_{1}$. While this is not always possible, we show
that we can add a suitable transition rate to have the Eqs. (42) be solvable
without altering the critical environment (see Appendix E). In other words, we
can define a parametrization $\Lambda$ such that

$\displaystyle\Lambda(x_{1},q_{1})$ $\displaystyle=\left\\{x_{1},x_{\setminus
1}=1,q_{1},q_{\setminus 1}=Q_{\setminus 1}(x_{1},q_{1})\right\\},$ (43)

where $x_{\setminus 1}=1$ denotes that for all $j>1$ we have $x_{j}=1$, and
similar for $q_{\setminus 1}$ we have $q_{j}=Q_{j}(x_{1},q_{1})$ with $Q_{j}$
as defined in Eq. 18. Applying this parametrization to the DPH, we obtain an
effective Hamiltonian $\mathcal{H}^{\star}=\mathcal{H}\circ\Lambda$ that only
depends on $x_{1}$ and $q_{1}$. However, note that the solutions of Eq. 42 can
yield different parametrizations. These may represent different states of the
system, where the corresponding reduced Hamiltonians have different hyperbolic
fixed points. To take the correct $\mathcal{H}^{\star}$, we can investigate
the stability of the fixed points of the full problem, and choose the
parametrization such that $\mathcal{H}^{\star}$ has the same attractor as the
complete system.

To be more specific, we define a domain of $\Lambda$ in the parameter space.
The points at which two different parametrizations provide the same reduced
dynamics are generally critical points (see also Appendix F.4). Moreover, for
many systems the conservation of the total population number hold, hence
reducing the dimension of the problem by one. So, by applying the reduction
procedure, we substitute one mean-field equation by Eq. 34. As before, we have
an extra degree of freedom since we can derive different $\Lambda$ depending
on the first integrals that we are choosing. However, that does not change the
deterministic description. In fact, despite the multiplicity of Hamiltonians
$\mathcal{H}^{\star}$, if the corresponding $\Lambda$ satisfies the condition
discussed above, the fixed points and the asymptotic behaviour of the system
do not vary. Otherwise, from a stochastic point of view, we have different
Langevin equations (see Eq. 25).

Our approach binds the system by imposing certain degrees of freedom fixed in
their steady state. This means that the stochastic information that we get is
related to the conditional probability $p(n_{1}|n_{2},..,n_{N-M})$, which we
are taking as an approximation for $p(n_{1})$ in the stationary limit.
Different parametrizations lead to different constraints, and thus to
different conditional probabilities. In this case we can write several
Langevin equations that will be only qualitative estimates of the stochastic
properties of the system, such as the stationary probability distribution. In
Appendix F.2 we discuss an example in more detail.

As a final remark on the parametrization $\Lambda$, we note that it is
associated to the ($x,q$)-coordinates. Indeed, it is defined from the mean-
field trajectories that are functions of the variables $q$ only. If we were to
use the Cole-Hopf variables of Eq. 72 instead, we would get a transformation
in terms of $\eta$, which is linked to the previous coordinates by
$\eta_{j}=x_{j}q_{j}$. Considering that we are imposing $x_{j}=1$ for $j>1$,
the two representations are equivalent only if we substitute $\langle
n_{1}\rangle$ in Eq. 13 with $x_{1}q_{1}$ instead of $q_{1}$. This would,
however, be incorrect, as we would go beyond the mean-field approximation and
additionally constrain fluctuations. Therefore, we must define $\Lambda$ by
using the $(x,q)$-coordinates, and only after doing this we can apply the
Cole-Hopf transformation and obtain a physical interpretation of the fields.
We want to emphasize that, while the reduction holds only at the large time
limit, it remains useful for the study of the system near the critical
condition. Similar to the work in Ref. [25], our goal is to create an analogy
with Landau theory, which, in equilibrium statistical mechanics, offers a way
to analyse the critical behaviour of a system [47]. Defining the effective
Landau free energy $\mathcal{L}$, one can study the phase transitions by
finding the minima of $\mathcal{L}$. By employing the Doi-Peliti formalism, we
can extend this approach for non-equilibrium systems, by introducing
$\mathcal{L}$ by $\mathcal{H}^{\star}$ and, instead of observing the minima of
the former, focusing on the fixed points of the latter instead.

We would like to note that our formalism does not pretend to give an exact
characterization of the system’s stochastic behavior. Instead, our approach
qualitatively describes this with an effective Langevin equation that we
obtain within a single framework. Within this framework, we derive both the
deterministic limit and investigate critical points. In the absence of this
ambiguity, our dimension reduction scheme proves to be a useful way to
describe the stochastic behavior of a variety of models, which will now be
applied to study CRNs in the remainder of this paper.

Before this, we would like to mention that deterministic studies of the
critical behaviour, such as the one we propose in Eq. 19, can additionally be
done with other approaches. For example, by obtaining a one-dimensional
effective mean-field equation using our method, we can define an effective
potential, as is done in Ref. [48], and investigate the phase transition by
studying its minima. However, our approach offers an alternative to
geometrically attack a broad class of problems. It does not only allow us to
go beyond the mean-field description, but also allows for using path-integral
representations to complex physics problems.

###  IV.2 Application to cascading failures: A modified $SIRS$ model

Let us begin displaying the strength of our proposed framework by considering
a modified version of the $SIRS$ model, which now, in contrast with the
previously analysed $SIS$ model, contains an additional compartment [49]. We
assume that the three compartments correspond to active, $A$, damaged
(failed), $F$, or inactive, $I$, states. The elementary reactions read

$\displaystyle I$ $\displaystyle\xrightarrow{\sigma}A,\quad$ $\displaystyle
A+F$ $\displaystyle\xrightarrow{\beta/N}2F,$ (44) $\displaystyle F$
$\displaystyle\xrightarrow{\gamma}I,\quad$ $\displaystyle A+2F$
$\displaystyle\xrightarrow{\alpha/N^{2}}3F,$

where $\sigma$ is the recovery rate of failed units, $\gamma$ the transition
rate of failed units to inactive ones, and $\alpha$ and $\beta$ denote the
“infection” reaction rates corresponding to three-body long-range failures
related to cascading failures, and the standard failure belonging to local
spread of failures. Hence, the model differs from a standard $SIRS$-model as
we include a three-body interaction111Note that this is similar to the three-
body interaction we introduced in the modified $SIS$-model that was discussed
in Section IV.2. with rate $\alpha/N^{2}$. This reaction can be interpreted as
the effect of a field mediating the higher-order effect and is the all-to-all
(fully-connected) limit of a model that, when studied in infrastructural or
social interconnected systems, is used to describe the non-local unfolding of
cascading failures. These cascading effects are observed, for instance, in
electrical power grids [50] or in communication networks [51].

We apply on this model the methodological pipeline defined above. The Doi-
Peliti Hamiltonian reads

$\displaystyle\begin{split}\mathcal{H}&=\frac{\beta}{N}a_{F}^{\dagger}(a_{F}^{\dagger}-a_{A}^{\dagger})a_{F}a_{A}\\\
&\quad+\frac{\alpha}{N^{2}}(a_{F}^{\dagger})^{2}(a^{\dagger}_{F}-a^{\dagger}_{A})a_{F}^{2}a_{A}\\\
&\quad+\gamma(a_{I}^{\dagger}-a_{F}^{\dagger})a_{F}+\sigma(a_{A}^{\dagger}-a_{I}^{\dagger})a_{I},\end{split}$
(45)

and the mean-field equations are given by

$\displaystyle\partial_{t}\langle n_{A}\rangle$
$\displaystyle=-\frac{\beta}{N}\langle n_{A}\rangle\langle
n_{F}\rangle-\frac{\alpha}{N^{2}}\langle n_{A}\rangle\langle
n_{F}\rangle^{2}+\sigma\langle n_{I}\rangle,$ (46)
$\displaystyle\partial_{t}\langle n_{F}\rangle$
$\displaystyle=\frac{\beta}{N}\langle n_{A}\rangle\langle
n_{F}\rangle+\frac{\alpha}{N^{2}}\langle n_{A}\rangle\langle
n_{F}\rangle^{2}-\gamma\langle n_{F}\rangle,$ (47)
$\displaystyle\partial_{t}\langle n_{R}\rangle$ $\displaystyle=\gamma\langle
n_{F}\rangle-\sigma\langle n_{I}\rangle.$ (48)

As done before, we set $\partial_{t}\langle n_{R}\rangle=0$, and by including
also the conservation of total population $N$, which we re-scale to $N=1$, we
obtain the transformations

$\displaystyle q_{F}=\frac{1-x_{A}q_{A}}{\gamma+\sigma}\sigma\,\,,\qquad
q_{R}=\frac{1-x_{A}q_{A}}{\gamma+\sigma}\gamma.$

The Cole-Hopf Hamiltonian related to the reduced problem is

$\displaystyle\mathcal{H}^{\star}$
$\displaystyle=-(\gamma+\sigma)^{-2}\biggr{\\{}(\eta-1)e^{-\theta}\left(e^{\theta}-1\right)\sigma\left[\alpha(\eta-1)\eta\sigma\right.$
$\displaystyle\left.\qquad\qquad\quad-\beta\eta(\gamma+\sigma)+\gamma
e^{\theta}(\gamma+\sigma)\right]\biggr{\\}},$ (49)

and the zero-energy lines are

$\displaystyle\eta$ $\displaystyle=1,\qquad\theta=0,$

and

$\displaystyle\theta$ $\displaystyle=\Theta(\eta,\beta,\gamma,\alpha,\sigma)$
$\displaystyle=\log\left(-\frac{\eta(\alpha\eta\sigma-\alpha\sigma-\beta\gamma-\beta\sigma)}{\gamma(\gamma+\sigma)}\right).$

These curves show the same structure as the ones of the modified $SIS$ model
(see Section IV.2), thus we recover the same topology reported in Fig. 1, and
the same phase portraits shown in Fig. A3. Therefore, we have three fixed
points on $\theta=0$: $F_{1}=(1,0)$ and $F_{\pm}$. We expect to $\gamma=\beta$
to be a critical point, whose nature is conditioned by the three-body
interaction mediated by $\alpha$. For $\alpha$ larger than a certain threshold
$\widetilde{\alpha}$, the system will display a first-order, abrupt
transition, while otherwise a second-order appears. By studying $\Theta$, we
note that the sign of the eigenvalues along the mean-field manifold changes at
$\beta=\gamma$ since $\Theta(\eta,\beta,\beta,\alpha,\sigma)=0$. In the range
$\beta>\gamma$, the fixed points $F_{\pm}$ become the new attractors. Let us
focus now on the attractor with $\eta\leq 1$. If at $\beta=\gamma$ we have
$F_{+}=F_{1}$, then the new attractor in the active phase is $F_{-}$, and so
there is an abrupt transition. Again, this happens when
$\alpha>\widetilde{\alpha}$. When $\alpha<\widetilde{\alpha}$, at the critical
point we have $F_{-}=F_{1}$, and the new attractor at $\beta>\gamma$ is still
$F_{-}$. In this case, the position of the attractor varies smoothly, and so
there is a second-order phase transition.

Further insights on the critical behavior can be obtained by investigation of
$\Theta$. Indeed, the tricritical condition, i.e., $\gamma=\beta$ and
$\alpha=\widetilde{\alpha}$ such that $F_{1}=F_{\pm}$, corresponds to the
solution of the following system

$\displaystyle\Theta(\eta,\beta,\gamma,\alpha,\sigma)$
$\displaystyle=\partialderivative{\eta}\Theta(\eta,\beta,\gamma,\alpha,\sigma)=0,$
$\displaystyle\eta$ $\displaystyle=1,$ $\displaystyle\beta$
$\displaystyle=\gamma$

We find that the tricritical point in the parameter space is given by

$\displaystyle\begin{split}\beta&=\gamma,\\\
\widetilde{\alpha}&=\gamma(\gamma+\sigma)/\sigma.\end{split}$ (50)

These results are consistent with the deterministic description of the system.
Furthermore, following the strategy discussed in Section III.2, we can
investigate the behavior beyond the mean-field limit. For this purpose, we
analyse the statistical fluctuations of the state of the system. Using the
expression of Eq. 49 for the reduced Hamiltonian $\mathcal{H}^{\star}$, the
reduced Langevin equation of Eq. 25 reads

$\displaystyle\mathop{}\\!\frac{\mathrm{d}\eta}{\mathrm{d}t}$
$\displaystyle=\left(\frac{\sigma(N-\eta)\left(\gamma
N^{2}(\gamma+\sigma)-N\eta(\alpha\sigma+\beta(\gamma+\sigma))+\alpha\sigma\eta^{2}\right)}{N^{2}(\gamma+\sigma)^{2}}\right)$
$\displaystyle\qquad\qquad\qquad+\left(\frac{\sigma(N-\eta)\left(\gamma
N^{2}(\gamma+\sigma)+N\eta(\alpha\sigma+\beta(\gamma+\sigma))-\alpha\sigma\eta^{2}\right)}{N^{2}(\gamma+\sigma)^{2}}\right)^{\frac{1}{2}}\cdot\xi.$
(51)

We compare between the order parameter for the reduced dynamics and for the
full dynamics carried out by the Kramers-Moyal expansion, and illustrate this
comparison in Fig. 3. We observe that, for $t$ large, the reduced mean-field
dynamics converges to the complete one, and the fluctuations display similar
statistical properties in both cases. Indeed, in Fig. 3 we give a comparison
between the stationary probability distribution function (PDF) of the full
system and the one derived from the reduced dynamics. In the case of reduced
dynamics, its PDF is given by Eq. 26. These figures indicate a good agreement
between our approximations and the real distribution, displaying the efficacy
of our approach.

Figure 3: Top: Numerical simulation for the modified $SIRS$ model. We have
considered $N=10^{3}$, $\beta=0.8$, $\gamma=0.4$, $\alpha=0.15$ and
$\sigma=1$. The dashed lines represent the deterministic dynamics for both the
reduced and the full system. The solid lines are the stochastic trajectories
generated by Eq. (25) applied to $\mathcal{H}$ and to $\mathcal{H}^{\star}$.
Bottom: Probability distribution function at the stationary limit of the
modified $SIRS$ model. The blue histogram represents the PDF obtained by the
long time behaviour of numerical simulations realized considering the reduced
dynamics. The orange one is the distribution related to the full system. The
solid blue line depicts the analytic reduced PDF given by Eq. (26). Relevant
parameters are $N=10^{3}$, $\beta=0.8$, $\gamma=0.4$, $\alpha=0.15$ and
$\sigma=1$.

This example highlights that we can analyse the critical environment of a
system by studying a simple function $\Theta$. This demonstrates one of the
advantages of our approach. Given a multidimensional system of which we want
to characterize some of its critical features, we can systematically follow
the steps of our formalism to reduce all the complexity to an analysis of some
parametric functions. That said, our approach is not limited to providing
insights into the critical properties, but also offers quantitative
assessments on the properties of the steady state, such as the value of the
order parameter and its distribution. In the next sections, we shall provide a
more detailed explanation of this.

###  IV.3 Application to ecology: the generalized Lotka-Volterra model

Let us now turn our attention to ecological systems. Within this context, one
most often considers the generalized Lotka-Volterra model, which describes the
dynamics of many interacting species. Recently, this model has been used to
investigate biodiversity, with the aim of providing an answer to the so-called
‘diversity-stability’-paradox, as theoretical approaches indicate that
diversity should hamper stability [52, 53, 54], yet observational studies find
the exact opposite must be true [55]. More recently, dynamical mean-field
theory has been applied to the generalized Lotka-Volterra model with great
success [56, 57], and there have been successful attempts at describing the
Lotka-Volterra model as a chemical reaction network [27]. Yet, to the best of
our knowledge, reduced dynamics (stochastic) Lotka-Volterra models have not
been put forward, but yet may provide valuable insights in (critical) behavior
of complex ecological models.

To this end, let us study the stochastic Lotka-Volterra model. This model can
be described by the following elementary reactions for each of the species
$I_{i}$,

$\displaystyle I_{i}$ $\displaystyle\xrightarrow{r_{i}}2I_{i},\quad$
$\displaystyle I_{i}$ $\displaystyle\xrightarrow{\mu_{i}}\emptyset,$ (52)
$\displaystyle\emptyset$ $\displaystyle\xrightarrow{\lambda_{i}}I_{i},\quad$
$\displaystyle I_{i}+I_{j}$ $\displaystyle\xrightarrow{b_{ij}}I_{j},$

where $r_{i}$, $\mu_{i}$ and $\lambda_{i}$ are, respectively, the growth,
death and migration rate of species $i$, and $b_{ij}$ is the interaction
coefficient that captures the effect of species $j$ on species $i$ (see, e.g.,
Ref. [27], for more details). The DPH of the Lotka-Volterra model reads

$\displaystyle\mathcal{H}$
$\displaystyle=\sum_{i}r_{i}a_{i}^{\dagger}(a_{i}^{\dagger}-1)a_{i}+\sum_{i}\lambda_{i}(a_{i}^{\dagger}-1)$
$\displaystyle\quad+\sum_{i}\mu_{i}(1-a_{i}^{\dagger})a_{i}+\sum_{i,j}b_{ij}a^{\dagger}_{i}(1-a_{j}^{\dagger})a_{i}a_{j}.$
(53)

####  IV.3.1 A simple two-species Lotka-Volterra system

As an example, we consider a simple two-species system consisting of native
and invasive species, denoted with $N$ and $I$, that compete. The elementary
reactions for this example read

$\displaystyle N$ $\displaystyle\xrightarrow{r}2N,$ $\displaystyle N$
$\displaystyle\xrightarrow{\mu}\emptyset,$ $\displaystyle 2N$
$\displaystyle\xrightarrow{b_{N}}N,$ $\displaystyle N+I$
$\displaystyle\xrightarrow{b}I,$ (54) $\displaystyle I$
$\displaystyle\xrightarrow{r}2I,$ $\displaystyle I$
$\displaystyle\xrightarrow{\mu}\emptyset,$ $\displaystyle 2I$
$\displaystyle\xrightarrow{b_{I}}I,$ $\displaystyle\emptyset$
$\displaystyle\xrightarrow{\lambda}I,$

where $b_{I}$ and $b_{N}$ are the self-regulation term corresponding to an
effective carrying capacity, $b$ is the effect of the invasive species on the
native one, and $\lambda$ is the migration rate (invasion rate) of the
invasive species. In this case, the DPH reads

$\displaystyle\mathcal{H}$
$\displaystyle=r\left(a^{\dagger}_{I}(a^{\dagger}_{I}-1)a_{I}+a^{\dagger}_{I}(a_{I}^{\dagger}-1)a_{I}\right)+\lambda(a^{\dagger}_{I}-1)$
(55)
$\displaystyle\quad+\mu\left((1-a^{\dagger}_{I})a_{I}+(1-a^{\dagger}_{N})a_{N}\right)$
$\displaystyle\quad+ba^{\dagger}_{I}(1-a_{N}^{\dagger})a_{N}a_{I}$
$\displaystyle\quad+b_{I}a^{\dagger}_{I}(1-a^{\dagger}_{I})a_{I}^{2}+b_{N}a^{\dagger}_{N}(1-a^{\dagger}_{N})a_{N}^{2}$

In order to proceed, we consider $b_{N},b_{I},b\ll 1$ to work in the range of
validity of our theory, i.e., the chemical reaction follow the scaling rule of
Sec. II.3. Moreover, we assume that $b_{I}>b$, that is, the carrying capacity
of the invasive species is not so large to dominate the native species for any
migration rate $\lambda>0$. We choose the asymptotic population of the native
species $N$ as the order parameter and derive the reduced Hamiltonian using
Eq. (42). We solve the system using the invasion rate $\lambda$ as a free
parameter, and we obtain a critical value

$\displaystyle\widetilde{\lambda}:=\frac{(r-\mu)^{2}(b_{I}-b)}{b^{2}}.$ (56)

The critical points for the reduced model are

$\displaystyle F_{1}$ $\displaystyle=(0,0),$ $\displaystyle F_{2}$
$\displaystyle=\left(\frac{r-\mu}{b_{N}}-\frac{b\left(\sqrt{(r-\mu)^{2}+4b_{I}\lambda}-\mu+r\right)}{2b_{I}b_{N}},0\right).$

For $\lambda>\widetilde{\lambda}$ the native species is driven to extinction
by the invasive one (Panel 4(c)), and $n_{B}\rightarrow 0$ for
$t\rightarrow\infty$. When $\lambda<\widetilde{\lambda}$ we find that the two
species can coexist (that is $n_{N}>0$) as shown in Panel 4(a). For lower
effective growth rates $r-\mu$, asymptotic abundances of the native species
reduce, indicating that slow growing species are more sensitive to invasion.
By careful investigation of the phase portrait (Fig. 4), it becomes clear that
at the critical point $\lambda=\widetilde{\lambda}$ the two fixed points
overlap and the native species is driven to extinction (see Fig. 4(b)). As the
equilibrium state changes smoothly, the extinction phase transition is of
second order. We report a numerical simulation in Fig. 4(d).

(a) $\lambda<\tilde{\lambda}$

(b) $\lambda=\tilde{\lambda}$

(c) $\lambda>\tilde{\lambda}$

(d)

Figure 4: Reduced dynamics of the Lotka-Volterra model. For
$\lambda<\tilde{\lambda}$ the attractor is $F_{2}$ with $n_{N}>0$. At the
critical condition we have the superposition of the fixed points
$F_{1}=F_{2}$, and the native species goes to the extinction for larger
$\lambda$. As confirmed by the numerical integration in Fig. (4(d)), this is a
second order phase transition. We have considered $b=0.75$, $b=0.7$,
$r=1/2500$, $r_{A}=1/2000$ and $r_{B}=1/3000$. With these values we have
$\tilde{\lambda}=1.563$.

####  IV.3.2 The Lotka-Volterra system with 8 species

Consider now a stochastic Lotka-Volterra model with a larger number of species
to illustrate how the dimension reduction can provide a very good
approximation for the asymptotic distributions. We take an 8-species $I_{i}$
model governed by the reactions of Eq. 52. We let $r_{i}=\lambda_{i}=1/2$ and
$\mu_{i}=0$, and specify the interaction matrix $\mathbf{B}$, with elements
$b_{ij}$, as

$\mathbf{B}=\frac{1}{1000}\begin{pmatrix}1&1&\frac{1}{2}&0&\frac{3}{10}&0&0&0\\\
\frac{1}{15}&\frac{11}{10}&0&0&\frac{1}{17}&0&0&0\\\
0&0&\frac{4}{5}&0&0&1&0&0\\\ 0&0&0&1&0&0&0&\frac{1}{9}\\\
0&0&0&0&\frac{11}{10}&0&0&0\\\ 1&0&0&0&0&1&0&0\\\
0&0&0&0&0&0&\frac{9}{10}&0\\\
\frac{1}{6}&0&0&0&0&0&\frac{1}{5}&1\end{pmatrix}.$ (57)

While we do not explore possible phase transitions here (for that, see [54,
56], among others), we will focus here on an analysis of the reduced Langevin
dynamics instead. The DPH for this system is given by Eq. 53. Without loss of
generality, we will use the first species as the order parameter. A Langevin
equation as in Eq. 25 can then be obtained, for which the stationary
distribution of the focal species can be derived by numerical integration
(Fig. 5). The results of our reduction procedure are in excellent agreement
with numerical integration of the full system, highlighting again the strength
of our proposed method.

Figure 5: Numerical asymptotic abundance distribution of the species $I_{1}$
of a system described by the generalized Lotka-Volterra model with $8$
species. Orange bars represent the normalized probability distribution
obtained by numerical integration of the one-dimensional reduced system, while
blue bars are results from numerical integration of the full system. We have
considered $r=\lambda=1/2$ and $\mathbf{B}$ given by Eq. 57.

###  IV.4 Resource-exploiter model

In this last section we treat an ecological model which describes the
interactions between resources $R$, exploiters $E$ and pollutants $X$. The
rationale behind this model is that there are resources that can be used and
transformed into something else (e.g., goods, infrastructure, etc.) because of
interactions with a population of exploiters. The process of transformation of
resources comes at some cost: it produces pollutants that in part can be re-
used to produce other resources. A typical example of the last process is the
production of CO2 due to human activities and the partial re-absorption of
this gas by plants to produce new resources (e.g., oxygen, plant-matter that
feeds herbivores, agricultural uses, etc.).

By denoting with $n_{R},\,n_{E}$ and $n_{X}$ the corresponding occupation
numbers, we can model such a system by assuming the following transition
rates:

$\displaystyle\omega$ $\displaystyle(\bm{n}+\mathbb{I}_{R}|\bm{n})=\alpha
n_{R}\,,$ $\displaystyle\omega$
$\displaystyle(\bm{n}+\mathbb{I}_{E},\bm{n})=\mu n_{E}\,,$
$\displaystyle\omega$
$\displaystyle(\bm{n}-\mathbb{I}_{R}|\bm{n})=n_{X}\xi+\frac{\alpha}{k_{R}}n_{R}(n_{R}-1)\,,$
$\displaystyle\omega$
$\displaystyle(\bm{n}-\mathbb{I}_{E}|\bm{n})=n_{E}\delta+\frac{\mu}{k_{E}}n_{E}(n_{E}-1)\,,$
$\displaystyle\omega$
$\displaystyle(\bm{n}-\mathbb{I}_{E}+\mathbb{I}_{X}|\bm{n})=\beta
n_{E}n_{R}\,,$ $\displaystyle\omega$
$\displaystyle(\bm{n}+\mathbb{I}_{R}-\mathbb{I}_{X}|\bm{n})=\beta\sigma
n_{R}n_{E}+\psi n_{X}n_{R}\,,$ $\displaystyle\omega$
$\displaystyle(\bm{n}+\mathbb{I}_{E}-\mathbb{I}_{X}|\bm{n})=\beta\gamma
n_{E}n_{R}\,.$

To reflect the constraints of empirical systems, we are considering a case
with resources limited by the capacity $k_{R}$. Exploiters consume resources
to their advantage while producing pollutants and this interaction is here
represented by a contact process: accordingly, the rate is proportional to
$n_{E}n_{R}$. By exploiting the resources $R$, $E$ increases, yet their
abundance is also limited a carrying capacity $k_{E}$. The carrying capacities
prevent unbounded population growth. Exploiters are further subjected to
natural birth and death processes with rates $\mu$ and $\delta$ respectively.
Resources grow with rate $\alpha$ (e.g., this generation can be due to the
injection of energy into the system, such as the one reaching the Earth which
is due to the sun) and can be recovered from $X$ by means of two contact
interactions: $\beta\sigma n_{E}n_{R}$, thanks to the exploiters, and $\psi
n_{X}n_{R}$. Finally, the pollutants reduce the number of resources by means
of the rate $\xi$.

(a) $\delta>\tilde{\delta}$

(b) $\delta=\tilde{\delta}$

(c) $\delta<\tilde{\delta}$

(d)

Figure 6: Phase transition for the resource-exploiter model. Figs. 6(a), 6(b)
and 6(b) show the phase portraits of the reduced system, where the orange line
is the mean-field manifold, the blue ones are the zero-energy lines of
$\mathcal{H}^{\star(A)}$ while the red curves are related to
$\mathcal{H}^{\star(B)}$. Red dots represent the fixed points on the mean-
field manifold. We observe the superposition of the attractos of
$\mathcal{H}^{\star(A)}$ and $\mathcal{H}^{\star(B)}$ at
$\delta=\tilde{\delta}$. The solid lines represent the effective phase
portrait, while the dashed ones are related to the alternative
parameterization that does not represent the system. In Panel. 6(d) we report
a numerical simulation of the model with $\beta=9\cdot 10^{-6}$,
$\gamma=8\cdot 10^{-2}$, $\sigma=0.1$, $\epsilon=0.9$, $\psi=10^{-5}$,
$\mu=0.1$, $\alpha=0.3$, $k_{R}=7\cdot 10^{5}$ and $k_{E}=8.5\cdot 10^{5}$.
The critical point is given by $\tilde{\delta}=0.604$.

We note that this system is hard to analyse using available tools, such as the
next-generation matrix. The reason is that the pollutants $X$, which one might
interpret as the “infected” counterpart as in standard epidemiological models,
can be generated by processes that are independent of their abundance (see the
system definition, and Eq. 61). In epidemiological terms, the “disease-free”
equilibrium is not stable, as the infected compartment can grow even when its
abundance has reached zero. This contradicts the main assumptions underlying
the next-generation matrix [15, 58]. By contrast, while our formalism does not
provide an interpretation of epidemic-related measures, such as the basic
reproduction number, it can be applied to systems where the main assumptions
underlying the next-generation matrix do not hold.

We proceed by deriving DPH for this system, which reads

$\displaystyle\mathcal{H}$ $\displaystyle=$ $\displaystyle\alpha
a_{R}^{\dagger}\left(a_{R}^{\dagger}-1\right)a_{R}+\mu
a_{E}^{\dagger}\left(a_{E}^{\dagger}-1\right)a_{E}$ (58) $\displaystyle+\xi
a_{X}^{\dagger}a_{X}\left((a_{R}^{\dagger})^{-1}-1\right)+\frac{\alpha}{k_{R}}a_{R}^{\dagger}\left(1-a_{R}^{\dagger}\right)a_{R}^{2}$
$\displaystyle+\frac{\mu}{k_{E}}a_{E}^{\dagger}\left(1-a_{E}^{\dagger}\right)a_{E}^{2}+\delta(1-a_{E}^{\dagger})a_{E}$
$\displaystyle+\beta
a_{E}^{\dagger}\left(a_{X}^{\dagger}-a_{R}^{\dagger}\right)a_{R}a_{E}$
$\displaystyle+\beta\sigma
a_{R}^{\dagger}\left(a_{R}^{\dagger}(a_{X}^{\dagger})^{-1}-1\right)a_{R}a_{E}^{\dagger}a_{E}$
$\displaystyle+\psi
a_{R}^{\dagger}\left(a_{R}^{\dagger}-a_{X}^{\dagger}\right)a_{R}a_{X}$
$\displaystyle+\gamma\beta
a_{R}^{\dagger}\left(a_{E}^{\dagger}(a_{X}^{\dagger})^{-1}-1\right)a_{E}a_{R}^{\dagger}a_{R}$

Using this expression, we obtain the following mean-field equations:

$\displaystyle\begin{split}\partial_{t}\langle n_{R}\rangle&=\alpha\langle
n_{R}\rangle\left(1-\frac{\langle n_{R}\rangle}{k_{R}}\right)-\beta\langle
n_{E}\rangle\langle n_{R}\rangle\\\ &\quad-\xi\langle
n_{X}\rangle+\beta\sigma\langle n_{R}\rangle\langle n_{E}\rangle+\psi\langle
n_{R}\rangle\langle n_{X}\rangle\end{split}$ (59)
$\displaystyle\partial_{t}\langle n_{E}\rangle$ $\displaystyle=\mu\langle
n_{E}\rangle\left(1-\frac{\langle
n_{E}\rangle}{k_{E}}\right)+\beta\gamma\langle n_{R}\rangle\langle
n_{E}\rangle-\delta\langle n_{E}\rangle$ (60)
$\displaystyle\partial_{t}\langle n_{X}\rangle$
$\displaystyle=(1-\sigma-\gamma)\beta\langle n_{R}\rangle\langle
n_{E}\rangle-\psi\langle n_{R}\rangle\langle n_{X}\rangle$ (61)

The WKB approximation is valid when the following conditions are satisfied:
(1) $k_{R},\,k_{E}\gg 1$, meaning that we consider relatively large
populations and that we are able to apply a system size expansion, and (2)
$\beta,\,\psi\ll 1$, in order to rescale the contact processes with respect to
the size of the system. We consider $\rho=1-\lim_{t\to\infty}\langle
n_{R}(t)\rangle/k_{R}$ as the order parameter. By letting $\partial_{t}\langle
n_{E}\rangle=\partial_{t}\langle n_{X}\rangle=0$ we obtain two different
parametrizations, $A$ and $B$, that read

$\displaystyle\Lambda_{A}=\left(x_{R},x_{E}=1,x_{X}=1,q_{R},q_{E}=0,q_{X}=0\right),$
$\displaystyle\Lambda_{B}=\Big{(}x_{R},1,1,q_{R},$
$\displaystyle\,\frac{k_{E}}{\mu}(q_{R}\beta\gamma-\delta+\mu),\frac{k_{E}\beta}{\mu\psi}(q_{R}\beta\gamma-\delta+\mu)(1-\gamma-\sigma)\Big{)}.$

By solving Eq. 19, we identify $\widetilde{\delta}=\mu+\gamma\beta k_{R}$ as a
critical point. By using $\Lambda_{A}$, the Cole-Hopf reduced Hamiltonian
reads

$\displaystyle\mathcal{H}^{\star}_{A}=\frac{\alpha\eta
e^{-\theta}\left(e^{\theta}-1\right)\left(e^{\theta}k_{R}-\eta\right)}{k_{R}},$
(62)

or, by employing $\Lambda_{B}$ instead, we find

$\displaystyle\begin{split}\mathcal{H}^{\star}_{B}=\frac{e^{-2\theta}\left(e^{\theta}-1\right)}{k_{R}\mu\phi}&\Big{(}\eta
e^{2\theta}k_{R}\phi(\alpha\mu+\beta(\gamma-1)k_{E}(\delta-\mu))+\beta^{2}\gamma\eta
k_{R}k_{E}(\epsilon(\gamma+\sigma-1)-\eta\phi)\\\
&\qquad\quad-e^{\theta}\big{(}\alpha\eta^{2}\mu\phi+\beta
k_{R}k_{E}(\eta\phi(\beta(\gamma-1)\gamma\eta-\delta+\mu)+\epsilon(\gamma+\sigma-1)(\delta-\mu))\big{)}\Big{)}\end{split}$
(63)

By studying the eigenvalues of the deterministic system, we find that $\langle
n_{R}\rangle=k_{R}$ is the attractor for $\delta>\widetilde{\delta}$ (Fig.
6(a)), otherwise the attractor is one of the mean-field fixed point of
$\mathcal{H}^{\star}_{B}$ as reported in Fig. 6(c). Following the discussion
of Section IV, $\mathcal{H}^{\star}_{A}$ is well defined for
$\delta>\widetilde{\delta}$, while $\mathcal{H}^{\star}_{B}$ represents the
system for $\delta<\widetilde{\delta}$. The energy levels of both reduced
Hamiltonians are shown in Fig. 6. From Figs. 6(a), 6(b) and 6(c), we note that
the attractors exchange their position in a continuous manner at the critical
point, which denotes a second-order transition. This is confirmed by numerical
simulations (Fig. 6(d)).

Since the DPH contains many parameters, it is difficult to explore all
possible ranges. Indeed, the system exhibits (numerical) instabilities, for
example when $\delta$ is very small. When Eqs. 59, 60 and 61 converge to an
equilibrium point (as is the case in Fig. 6), and thus satisfies the stability
condition necessary for our approach, we see instead that we are able to fully
characterize the critical behaviour using our method.

##  V Conclusions

In this work we have presented a unifying framework that shows how to
systematically apply the Doi-Peliti formalism to a broad class of
compartmental models. Several studies have proven how path integral approaches
can offer advantages when describing stochastic systems, characterizing phase
transitions and giving information beyond the deterministic limit. Yet,
dealing with systems with many compartments is generally complicated,
especially regarding the visualization of the geometry of phase portraits to
identify the nature of the phase transitions. Here we have derived a novel
strategy that allows us to unify these descriptions thanks to the introduction
of a dimension reduction scheme. In practice, through our method we gain
information about the nature of the phase transitions, and analyze both
qualitatively and quantitatively the stochastic dynamics. We do so through a
transparent physical interpretation: since the observable of interest, the
order parameter, corresponds to a single compartment, one can freeze all other
degrees of freedom and study only the compartment corresponding to the
observable.

We have illustrated the usefulness of this methodology by applying it to
several multi-compartmental stochastic processes of increasing complexity and
of different nature, including epidemiological, ecological and sociotechnical
ones. Starting from the simplest case of two-compartment systems, we have
explored the role that three-body interactions, recently reported to
significantly impact the behavior of many empirical systems [59], play on the
critical aspects of epidemiological-inspired models. This opens the door to
systematically analyze systems whose interactions are beyond pairwise. For
systems with more than two components, our dimension reduction proposal offers
very good approximations for the stationary probability density function of
the order parameter, and successfully captures their critical properties.

To derive our dimension reduction scheme, we have used some approximations
that do not impact the predictions of the critical behavior, yet they do have
an effect on the transient dynamics. Indeed, many methods of dimensionality
reduction aim to describe not only the stationary limit, but also to provide
information about the whole dynamics, which is however often at the expense of
needing more assumptions. Hence, there is a compromise between the amount of
information one can analytically retrieve and the generality of the model one
can actually study. Our method positions itself in a place in which it can be
generally applied to a range of systems that is wider than what traditional
methods, such as those based on time scale separation [60, 61], offer, but it
misses the possibility of an accurate prediction of the transient phase.

Offering insights on the transient dynamics from the point of view developed
here remains as a future research perspective. Additionally, our method is
most useful in scenarios in which only one degree of freedom is effectively
operating in the long-time behaviour, so it not yet clear how to reduce the
dimensionality in a case with coupled degrees of freedom. Another aspect that
we believe is worth further consideration is related to the parametrization
$\Lambda$ and the reduced Hamiltonian,
$\mathcal{H}^{\star}=\Lambda\circ\mathcal{H}$. More specifically, the
difficulties are associated with the definition of the validity domain and the
potential degeneration of $\Lambda$. A remaining open question, which we
foresee having a great impact due to the wealth of potential applications, is
the generalization of this dimension reduction scheme to structured, networked
populations [62].

All in all, we have exploited the combination of the path-integral
representation of stochastic processes with a dimension reduction scheme to
better understand the critical properties and the stationary behavior of
coupled compartmental stochastic processes. We hope that our formalism will
provide new perspectives and physical insights into, and beyond, the
characterization of phase transitions and fluctuation dynamics multi-
compartment complex systems.

## Acknowledgments

O.A. acknowledges financial support from the Spanish Ministry of Universities
through the Recovery, Transformation and Resilience Plan funded by the
European Union (Next Generation EU) and the University of the Balearic
Islands. J.N. acknowledges financial support from the Human Frontier Science
Program Organization (HFSP Ref. RGY0064/2022). M.D.D. acknowledges partial
financial support from the MUR - PNC (DD n. 1511 30-09-2022) Project no.
PNC0000002, DigitAl lifelong pRevEntion (DARE); from MUR funding within the
PRIN 2022 PNRR (DD n. 1214 31-07-2023) Project no. P2022A889F and from MUR
funding within the FIS (DD n. 1219 31-07-2023) Project no. FIS00000158. T.S.
acknowledges financial support from the University of Padua (PRD-BIRD 2022).

##  Appendix A Derivation of the path integral representation

Here we show the derivation of the path integral representation for the
generating function $\ket{G}$ by using the Doi-Peliti formalism as done in
[24]. Let us take the simplest one-compartment case $\ket{\bm{n}}=\ket{n}$
with the basis (9). We take also $\mathcal{H}$ in its normal ordered form. By
writing the time dependency explicitly $\ket{G}=\ket{G(t)}_{x}$ and by
splitting $t$ in infinitesimal intervals $t=M\epsilon$ with $\epsilon\to 0$
and $M\to 0$, we have

$\displaystyle\ket{G(t)}_{x}=e^{\epsilon\mathcal{H}(a^{\dagger},a)}\ket{G(t-\epsilon)}_{x}+O(\epsilon)$
$\displaystyle=\int_{\mathbb{R}}\mathop{}\\!\mathrm{d}{x_{M-1}}e^{\epsilon\mathcal{H}(x,\partial_{x})}\delta(x-x_{M-1})\ket{G(t-\epsilon)}_{x}+(\epsilon^{2})$
$\displaystyle=\int_{\mathbb{R}^{2}}\frac{\mathop{}\\!\mathrm{d}{x_{N-1}}\mathop{}\\!\mathrm{d}{q_{N-1}}}{2\pi}\,e^{\epsilon\mathcal{H}(x_{N},\mathrm{i}q_{N-1})}$
$\displaystyle\qquad\qquad\times
e^{-\mathrm{i}q_{N-1}(x_{N-1}-x_{N})}\ket{G(t_{N-1})}_{x_{N-1}}+O(\epsilon^{2})$
$\displaystyle=\int_{\mathbb{R}^{2N}}\prod_{j={0}}^{N-1}\frac{\mathop{}\\!\mathrm{d}{x_{j}}\mathop{}\\!\mathrm{d}{q_{j}}}{2\pi}\,e^{\epsilon\sum_{j=0}^{M-1}\Big{(}\mathcal{H}(x_{j+1},\mathrm{i}q_{j})-\mathrm{i}q_{j}\frac{x_{j}-x_{j+1}}{\epsilon}\Big{)}}$
$\displaystyle\qquad\qquad\qquad\qquad\times\ket{G(t_{0})}_{x_{0}}+O(\epsilon^{2})$
$\displaystyle=\int_{[t_{0}}^{t)}\mathcal{D}[x,q]\,e^{-\mathcal{S}[x,\mathrm{i}q]}\ket{G(t_{0})}_{x_{0}}$

where in the last equation we have taken the continuous limit by introducing
the measure

$\int_{[t_{0}}^{t)}\mathcal{D}[x,q]=\lim_{N\to\infty}\int_{\mathbb{R}^{2N}}\prod_{i={0}}^{N-1}\frac{\mathop{}\\!\mathrm{d}{x_{i}}\mathop{}\\!\mathrm{d}{q_{i}}}{2\pi},$

and $\mathcal{S}$ is the action defined by Eq. (12). For the multidimensional
case the derivation is trivially the same.

##  Appendix B WKB approximation for CRNs

Here, we follow the reasoning of [26] and demonstrate that similar arguments
can be applied for extend the WKB to the case of chemical reaction networks
(3) as well. Let us denote by $\omega_{k}$ the transition rate related to the
reaction between $k$ reagents, and by $N$ the total population in the system.
For illustrative purposes, let us consider a simple contact process

$\displaystyle I_{j}+I_{j}\xrightarrow{\beta}I_{j}.$

The transition rate is proportional to the number of elements, $n_{j}$, times
the probability of finding another element in the heterogeneous network,
$(n_{i}-1)/N$. Thus, assuming $n_{j},N\gg 0$ and $n_{j}\ll N$, we have
$\omega_{2}(\textit{{n}}-\mathbb{I}_{j}|\textit{{n}})=\beta
n_{j}(n_{j}-1)/N\propto n_{j}^{2}/N$. Generalizing the same reasoning for a
contact process between $k$ elements, the transition rate is
$\omega_{k}\propto n_{j}^{k}/N^{k-1}$. Since the Doi-Peliti Hamiltonian is a
polynomial function of $x_{j}$ and $\mathrm{i}q_{j}$ , we conclude that the
term related to the highest power of $\mathrm{i}q_{j}$ inside $\mathcal{H}$ is
proportional to
$(\mathrm{i}q_{h})^{k}N^{-k+1}\,(=N^{-k+1}\partial^{k}_{x_{j}})$. Furthermore,
with these rates we can write the DPH and the action (12) as extensive
quantities with respect to $N$ as done in [26]. If we define the concentration
variables $\eta_{j}=n_{j}/N$, then the transition rate reads

$\omega_{k}\approx\beta\eta_{j}^{k}N.$

By means of Eq. 2 and Eq. 8, we see that

$H_{n,m}=Nh_{m,m}(\eta)\,,\qquad\mathcal{S}=Ns(\eta),$

with $\\{h_{n,m}\\}_{n,m}$ and $s$ being functions of the concentrations only.
Using the ansatz $\ket{G}=\exp(-Ns)$, we obtain

$\frac{\partial_{x_{j}}^{k}\exp(-Ns)}{N^{k-1}}=N(-\partial_{x_{j}}s)^{k}+O(1)\approx
N(-\partial_{x_{j}}s)^{k}.$

Inserting the same ansatz in Eq. (7), the relation
$\mathrm{i}q_{j}=-\partial_{x_{j}}\mathcal{S}$ holds for any $j$ and we find
the Hamilton-Jacobi equation

$\partial_{t}\mathcal{S}=-\mathcal{H}\left(\textit{{x}},-\partial_{\textit{{x}}}\mathcal{S}\right),$
(64)

with
$\partial_{\textit{{x}}}\mathcal{S}=(\partial_{x_{1}}\mathcal{S},\partial_{x_{2}}\mathcal{S},\ldots)$.
Therefore, for CRNs we can always take the stationary action defined by the
solutions of Eq.(13) and $\ket{G}=\exp(-\mathcal{S})$.

##  Appendix C A basis for the stationary paths

Now we report the basis for the WKB approximation presented in [24].
Considering the one-dimensional case and denoting with $\tilde{x}(t)$ and
$\tilde{q}(t)$ the solutions of the Hamilton’s equations (13), we take the
following basis

$\displaystyle\ket{\tilde{n}}=\big{(}\Delta
x+\tilde{x}(t)\big{)}^{n}e^{-\tilde{q}(t)(\Delta x+\tilde{x}(t))}$
$\displaystyle\bra{\tilde{n}}=\big{(}\partial_{\Delta
x}+\tilde{q}(t)\big{)}^{n},$

where $\Delta x$ and $\mathrm{i}\Delta q$ are fluctuations of $x(t)$ and
$\mathrm{i}q(t)$ around the stationary trajectories. The creation and
annihilation operators act according to

$\displaystyle
a\ket{\tilde{n}}=n\ket{\tilde{n}-1}=(\partial_{x}+\tilde{q}(t))\ket{\tilde{n}}$
(65) $\displaystyle
a^{\dagger}\ket{\tilde{n}}=\ket{\tilde{n}+1}=(x+\tilde{x}(t))\ket{\tilde{n}}.$
(66)

With this basis, in the derivation of the path integral representation we
still have to include the contribution resulting from the time dependency of
$\ket{\tilde{n}}$. For the sake of simplicity, we write
$\tilde{x}(t)=\tilde{x}$ and $\tilde{q}(t)=\tilde{q}$, and by taking the time
derivative of $\ket{\tilde{n}}$ we have

$\displaystyle\partial_{t}\ket{\tilde{n}}$
$\displaystyle=\Big{[}n\partial_{t}\tilde{x}(\Delta x+\tilde{x}^{n-1})$
$\displaystyle\quad\,-(\Delta x+\tilde{x})^{n}(\partial_{t}\tilde{q}(\Delta
x+\tilde{x})+\tilde{q}\partial_{t}\tilde{x})\Big{]}e^{-\tilde{q}(\Delta
x-\tilde{x})}$
$\displaystyle=n\partial_{t}\tilde{x}\ket{\tilde{n}-1}-\left(\partial_{t}\tilde{q}(\Delta
x+\tilde{x})+\tilde{q}\partial_{t}\tilde{x}\right)\ket{\tilde{n}}$
$\displaystyle=\left(\partial_{t}\tilde{x}\,a-\partial_{t}\tilde{q}\,a^{\dagger}-\tilde{q}\partial_{t}\tilde{x}\right)\ket{\tilde{n}}$
$\displaystyle=:\mathcal{B}\ket{\tilde{n}}$ (67)

where $\mathcal{B}$ is the time evolution operator of the basis. Within this
representation, the ME (1) becomes

$\displaystyle\partial_{t}\ket{G}$
$\displaystyle=\left(\mathcal{H}+\mathcal{B}\right)\ket{G}.$

Now, by replacing $(x,\,\mathrm{i}q)$ by $(\Delta x,\,\mathrm{i}\Delta q)$,
and by substituting $\mathcal{H}$ with $\mathcal{H}+\mathcal{B}$ in the steps
presented of Appendix A, we obtain the path integral representation with this
new basis. In detail, the action reads

$\displaystyle-$ $\displaystyle\mathcal{S}[\tilde{x}+\Delta
x,\tilde{q}+\mathrm{i}\Delta q]$
$\displaystyle=\int_{t_{0}}^{t}\mathop{}\\!\mathrm{d}{\tau}\,\mathrm{i}\Delta
q\partial_{\tau}\Delta x+\mathcal{H}(\tilde{x}+\Delta
x,\tilde{q}+\mathrm{i}\Delta q)$ (68)
$\displaystyle\qquad\qquad+\mathcal{B}(\tilde{x}+\Delta
x,\tilde{q}+\mathrm{i}\Delta q)$
$\displaystyle=\int_{t_{0}}^{t}\mathop{}\\!\mathrm{d}{\tau}\,\Big{[}\mathrm{i}\Delta
q\partial_{\tau}\Delta x+\mathcal{H}(\tilde{x},\tilde{q})+\Delta
x\partial_{\Delta x}\mathcal{H}(\tilde{x},\tilde{q})$
$\displaystyle\,\,\qquad+\mathrm{i}\Delta q\partial_{\mathrm{i}\Delta
q}\mathcal{H}(\tilde{x},\tilde{q})+\Delta\mathcal{H}+\tilde{q}\partial_{\tau}\tilde{x}+\mathrm{i}\Delta
q\partial_{\tau}\tilde{x}$ $\displaystyle\qquad\qquad\quad\quad-\Delta
x\partial_{\tau}\tilde{q}-\tilde{x}\partial_{\tau}\tilde{q}-\tilde{q}\partial_{\tau}\tilde{x}\Big{]}$
$\displaystyle=\int_{t_{0}}^{t}\mathop{}\\!\mathrm{d}{\tau}\,\Big{[}\mathrm{i}\Delta
q\partial_{\tau}\Delta x+\mathcal{H}(\tilde{x},\tilde{q})+\Delta
x\partial_{\tau}\tilde{q}-\mathrm{i}\Delta q\partial_{\tau}\tilde{x}$
$\displaystyle\qquad\qquad\quad\quad+\Delta\mathcal{H}+\mathrm{i}\Delta
q\partial_{\tau}\tilde{x}-\Delta
x\partial_{\tau}\tilde{q}-\tilde{x}\partial_{\tau}\tilde{q}\Big{]}$
$\displaystyle=\int_{t_{0}}^{t}\mathop{}\\!\mathrm{d}{\tau}\,\Big{[}\mathcal{H}(\tilde{x},\tilde{q})+\tilde{q}\partial_{\tau}\tilde{x}\Big{]}-\tilde{x}(t)\tilde{q}(t)+\tilde{x}_{0}\tilde{q}_{0}$
$\displaystyle\qquad\quad\qquad\quad+\int_{t_{0}}^{t}\mathop{}\\!\mathrm{d}{\tau}\,\Big{[}\mathrm{i}\Delta
q\partial_{\tau}\Delta x+\Delta\mathcal{H}\Big{]}$
$\displaystyle=-\tilde{\mathcal{S}}-\tilde{x}(t)\tilde{q}(t)+\tilde{x}_{0}\tilde{q}_{0}$
$\displaystyle\qquad\qquad\,\,+\int_{t_{0}}^{t}\mathop{}\\!\mathrm{d}{\tau}\,\Big{[}\mathrm{i}\Delta
q\partial_{\tau}\Delta x+\Delta\mathcal{H}\Big{]}$ (69)

where in the second equality we have expanded $\mathcal{H}(\tilde{x}+\Delta
x,\tilde{q}+\mathrm{i}\Delta q)$ around the stationary trajectories, and
$\Delta\mathcal{H}$ includes all the contributions of order $O(\Delta
x^{2},\Delta q^{2},\Delta x\Delta q)$. Here $\tilde{\mathcal{S}}$ satisfies
the Hamilton-Jacobi equation (64).

Let us now rewrite the generating function in terms of the basis
$\ket{\tilde{n}}$,

$\displaystyle\ket{G(t)}_{x}=e^{\tilde{x}(t)\tilde{q}(t)}\sum_{n}p(n,t)\ket{\tilde{n}}\Big{|}_{\Delta
x=0}.$

Thus, the path integral expression is

$\displaystyle\ket{G(t)}_{x}=$ $\displaystyle
e^{-\tilde{\mathcal{S}}}\int_{[t_{0}}^{t)}\mathcal{D}[\Delta x,\Delta
q]\,e^{\int_{t_{0}}^{t}\mathop{}\\!\mathrm{d}{\tau}\,\big{(}\mathrm{i}\Delta
q\partial_{\tau}\Delta
x+\Delta\mathcal{H}\big{)}}\ket{n_{0}}\Big{|}_{\genfrac{}{}{0.0pt}{}{\Delta
x(t)=0}{\tilde{x}(t)=x}}$

If we assume fixed initial condition, then
$\ket{n}_{0}=\ket{\tilde{n}_{0}}e^{\tilde{x}_{0}\tilde{q}_{0}}=e^{n_{0}\log{\tilde{x}_{0}}}$.
Within the WKB approximation, the fluctuating part does not contribute and we
have $\ket{G}=\exp\big({-\tilde{\mathcal{S}}}\big{missing})$ where $x$ and
$\mathrm{i}q$ are replaced by the real trajectories $\tilde{x}$ and
$\tilde{q}$.

##  Appendix D Doi-Peliti Hamiltonian for Complex Reaction Networks

In this section, we present the generic structure of a DPH in the case of
chemical reaction networks, providing a proof of the relation (23). For the
systems treated in this work, we can distinguish five fundamental processes
and, for each of them, we associate a term of the DPH by using definitions (5)
and (8). These read

$\displaystyle\emptyset\to A$ $\displaystyle a_{A}^{\dagger}-1$ $\displaystyle
A\to 2A\qquad$ $\displaystyle a_{A}^{\dagger}(a_{A}^{\dagger}-1)a_{A}$
$\displaystyle A\to\emptyset\qquad$ $\displaystyle(1-a_{A}^{\dagger})a_{A}$
$\displaystyle A\to B\qquad$
$\displaystyle(a_{B}^{\dagger}-a_{A}^{\dagger})a_{A}$ $\displaystyle A+B\to
A+C\qquad$ $\displaystyle
a_{A}^{\dagger}(a_{C}^{\dagger}-a_{B}^{\dagger})a_{A}a_{B}.$

These derivations can easily be generalized by more complicated interactions
such as the ones introduced in Sec. IV.2. For a CRN, generally we can build
the DPH as a summation of the elements presented above.

Knowing the structure of the DPH, we can now prove Eq. (23). For a two-
compartment model, the Kramers-Moyal expansion of the ME reads

$\displaystyle\partial_{t}p(\bm{n},t)=\sum_{\begin{subarray}{c}m_{1},m_{2}\\\
m_{1}+m_{2}>0\end{subarray}}\frac{1}{m_{1}!m_{2}!}\partial^{m_{1}}_{y_{1}}\partial^{m_{2}}_{y_{2}}W^{m_{1},m_{2}}(\bm{n})p(\bm{n},t),$

where

$W^{m_{1},m_{2}}(\bm{n})=\int\mathop{}\\!\mathrm{d}{y_{1}}\mathop{}\\!\mathrm{d}{y_{2}}\,\omega(\bm{n}+y_{1}+y_{2}|\bm{n})y_{1}^{m_{1}}y_{2}^{m_{2}}.$

Taking a two-body interaction

$pA+qB\xrightarrow{\gamma}(p+k)A+(q-l)B,$

the transition rate can be written as

$\displaystyle\omega(\bm{n}+y_{1}+y_{2}|\bm{n})=\gamma
n_{A}^{p}n_{B}^{q}\delta(y_{1}+k)\delta(y_{2}+l).$ (70)

Hence

$W^{m_{1},m_{2}}(\bm{n})=\gamma n_{A}^{p}n_{B}^{q}k^{m_{1}}(-l)^{m_{2}},$

and we find an expression for $\mathcal{H}^{\dagger}_{KM}$

$\displaystyle\mathcal{H}^{\dagger}_{KM}(\bm{\eta},\mathrm{i}\bm{\theta})$
$\displaystyle=\sum_{\begin{subarray}{c}m_{1},m_{2}\\\
m_{1}+m_{2}>0\end{subarray}}\frac{(\mathrm{i}\theta_{1}k)^{m_{1}}(-\mathrm{i}\theta_{2}l)^{m_{2}}}{m_{1}!m_{2}!}\gamma
n_{A}^{p}n_{B}^{q}$ $\displaystyle=\gamma
n_{A}^{p}n_{B}^{q}\left(e^{\mathrm{i}\theta_{1}k-\mathrm{i}\theta_{2}l}-1\right)$
(71)

On the other hand, the DPH for this reaction is

$\displaystyle\mathcal{H}$
$\displaystyle=\gamma(a^{\dagger}_{A})^{p}(a^{\dagger}_{B})^{q-l}\left((a^{\dagger}_{A})^{k}-(a_{B}^{\dagger})^{l}\right)a_{A}^{p}a_{B}^{q}$
$\displaystyle=\gamma\left((a^{\dagger}_{A})^{k}(a_{B}^{\dagger})^{-l}-1\right)(a^{\dagger}_{A})^{p}a_{A}^{p}(a^{\dagger})^{q}_{B}a_{B}^{q}$
(72)

By applying the Cole-Hopf transformation, we substitute
$a^{\dagger}_{A,B}=e^{\mathrm{i}\theta_{A,B}}$ and
$a_{A,B}=e^{-\mathrm{i}\theta_{A,B}}\eta_{A,B}$ into Eq. (72) obtaining the
same expression of the Kramers-Moyal expansion (71). These arguments can be
extended to a more complex and general interactions

$\sum_{i=1}^{M}p_{i}\,I_{i}\xrightarrow{\gamma}\sum_{i=1}^{M}(p_{i}+k_{i})\,I_{i},$

with $M$ representing the number of reagents/compartments $\\{I_{i}\\}$ and
$k_{i}$ can be any integer. In this case we define the vector
$\bm{k}=(k_{1},k_{2},\ldots,k_{M})$ and the transition rate is given by

$\omega(\bm{n}+\bm{y}|\bm{n})=f(\bm{n})\prod_{j=1}^{M}\delta(y_{j}-k_{j}),$

where $f$ a suitable function of the occupation numbers. Thus

$\displaystyle W^{m_{1},m_{2},\ldots,m_{M}}(\bm{n})$
$\displaystyle=\int\mathop{}\\!\mathrm{d}{\bm{y}}\left(\prod_{j=1}^{M}y_{j}^{m_{j}}\right)\omega(\bm{n}+\bm{y}|\bm{n})f(\bm{n})$
$\displaystyle=f(\bm{n})\prod_{j=1}^{M}k_{j}^{m_{j}},$

and the adjoint of the Kramers-Moyal operator is

$\displaystyle\mathcal{H}^{\dagger}_{KM}(\bm{\eta},\mathrm{i}\bm{\theta})$
$\displaystyle=\sum_{\begin{subarray}{c}m_{1},m_{2},\ldots,m_{M}\\\
\sum_{j}m_{j}>0\end{subarray}}\prod_{j=1}^{M}\frac{\left(\mathrm{i}\theta_{j}k_{j}\right)^{m_{j}}}{m_{j}!}f(\bm{\eta})$
$\displaystyle=\left(\prod_{j=1}^{M}e^{\mathrm{i}\theta_{j}k_{j}}-1\right)f(\bm{\eta}).$
(73)

In terms of the Doi-Peliti representation, the corresponding DPH is

$\mathcal{H}=\gamma\sum_{i=1}^{M}\left((a^{\dagger}_{i})^{k_{i}}-1\right)f(a^{\dagger}a),$
(74)

which can also be written as a function of the Cole-Hopf variables,

$\mathcal{H}(\bm{\eta},\mathrm{i}\bm{\theta})=\gamma\left(\prod_{j=1}^{M}e^{\mathrm{i}\theta_{j}k_{j}}-1\right)f(\bm{\eta}).$

This is equal to the KM operator (73). Eq. (74) shows that the mean-field
manifold $\bm{q}=0$ ($\bm{\theta}=0$) is always a zero of $\mathcal{H}$.

By means of the relation (23), we can show the link with the Langevin equation
(25). The procedure for deriving Langevin dynamics from $\mathcal{H}_{KM}$ has
been discusses in several works (e.g. [24] and [27]) and, for completeness, we
sketch the derivation in the following.

In the simple one-dimensional case, we can develop a path-integral
representation of the probability distribution starting from the Kramers-Moyal
expansion (21). Following the same steps reported in Appendix A, we write

$p(n,t)=\int_{[0}^{t)}\mathcal{D}[\eta,\theta]\,e^{-\mathcal{S}^{\dagger}_{KM}[\eta,\mathrm{i}\theta]}\delta(n-\eta(t)),$

where the action is defined by

$\mathcal{S}^{\dagger}_{KM}=\int_{0}^{t}\mathop{}\\!\mathrm{d}{\tau}\,\mathrm{i}\theta(\tau)\partial_{\tau}\eta(\tau)-\mathcal{H}^{\dagger}_{KM}(\eta(\tau),\mathrm{i}\theta(\tau)).$

As discussed in the main text, we know that the deterministic limit is
recovered for $\theta\to 0$. To obtain the Langevin expression, we expand
$\mathcal{S}^{\dagger}$ around $\theta=0$ up to the second order and integrate
with respect $\theta(\tau)$,

$\displaystyle
p(n,t)=\int_{[0}^{t)}\mathcal{D}[\eta]\,e^{-\int_{0}^{t}\frac{(\partial_{\tau}\eta-\partial_{\theta}\mathcal{H}^{\dagger}_{KM}|_{\theta=0})^{2}}{2\partial^{2}_{\theta}\mathcal{H}_{KM}^{\dagger}|_{\theta=0}}}\delta(n-\eta(t)),$

where $\mathcal{D}[\eta]$ is the re-normalized measure. The latter expression
is nothing but the Feynman-Kac formula [24]

$p(n,t)=\langle\delta(n-\eta(t))\rangle_{\eta},$

where $\langle\cdot\rangle_{\eta}$ is the average over the trajectories
defined by

$\mathop{}\\!\frac{\mathrm{d}\eta(t)}{\mathrm{d}t}=\partial_{\theta}\mathcal{H}^{\dagger}_{KM}\big{|}_{\theta=0}+\sqrt{\partial^{2}_{\theta}\mathcal{H}^{\dagger}_{KM}\big{|}_{\theta=0}}\cdot\xi(t)$

where $\xi(t)$ is the standard white noise. Including the relation (23), we
recover Eq. (25). The generalization to the multidimensional case is
straightforward.

##  Appendix E Closure transition rates

Here we discuss some limits of our formalism. In the main text, we have
considered models such as the $SIS$ or $SIRS$ that have a closed diagrammatic
structure, and systems where each compartment is subject to at least two
interactions. Indeed, as we will show, these features provide the necessary
conditions to apply the Doi-Peliti mapping.

Let us take a chemical reaction model with $M$ compartments
$\\{I_{1},I_{2},\ldots,I_{M}\\}$ and where only spontaneous transitions and
contact interactions occur. We assume $I_{M}$ is connected to other
compartments by a single spontaneous transition
$I_{j}\xrightarrow{\gamma}I_{M}$. Thus, $I_{M}$ is a termination of the model
diagram. Moreover, we consider to work with fixed $N$ and, as for the $SIS$
model, we take $I_{1}\approx N$. We choose the order parameter $\rho$ as in
Eq. (17). Now, we focus on $I_{1}$ and $I_{M}$ to understand if and when a
phase transition occurs. We rewrite the diagram of the system as

(a)

(b)

Figure A1: $SIR$-like diagram for the effective open process (a) and
$SIRS$-like diagram for the effective closed process (b).

where we have replaced the $M-2$ compartments $I_{2},\ldots,I_{M-1}$ by an
effective one $I_{\text{eff}}$. The effective transition rates
$\omega_{I_{1}}$ and $\omega_{I_{\text{eff}}}$ are functions of all
populations except $n_{M}$. Thus, the system is active when $I_{\text{eff}}$
becomes asymptotically populated and we can write

$\partial_{t}\langle
n_{I_{\text{eff}}}\rangle(t)=f(\\{\beta\\},I_{1},I_{\text{eff}}),$ (75)

with $\\{\beta\\}$ representing the transition coefficients of the model, and
$f$ is a suitable function depending on the system we are considering. The
problem is that the mean-field equations include

$\partial_{t}\langle n_{M}\rangle=\gamma\langle n_{j}\rangle,$ (76)

from which, at the stationary limit, $\langle n_{j}\rangle=0$. The presence of
null terms like this leads to a degeneration of the system (19). Thus
$\text{det}(\textbf{J})=0$ for any fixed point and we cannot use the strategy
presented in the main text to analyse the critical behaviour. To overcome this
problem, we introduce an effective transition rate that closes the model
diagram connecting $I_{M}$ directly to $I_{1}$ by means of a spontaneous
transition $\omega(\bm{n}-\mathbb{I}_{M}+\mathbb{I}_{I_{1}}|\bm{n})$. In doing
so we obtain a $SIRS$-like model, as depicted in 1(b).

From now on we will call open a system described by diagrams like Fig. 1(a),
and closed a system represented by diagrams like Fig. 1(b). By following a
classical treatment [13], one can prove that both the $SIR$ and $SIRS$ models,
have the same critical condition. Indeed, both of them are active when the
number of elements in state $I$ grows during the early stage of the process.
For $t\approx 0$ we have $n_{S}\approx N$, and from the mean-field equations
it is easy to find the critical condition

$\partial_{t}\langle n_{I}\rangle=\frac{\beta}{N}\langle n_{I}\rangle\langle
n_{S}\rangle-\gamma\langle n_{I}\rangle\approx(\beta-\gamma)\langle
n_{I}\rangle$ (77)

Therefore, the epidemic spreads when initially the population of infected
individuals manages to grow fast enough, that is, when $\beta>\gamma$. Since
Eq. (77) holds for both $SIR$ and $SIRS$, we have that the critical condition
does not change.

Focusing back on a general process, we can repeat the same steps to find that

$\partial_{t}\langle
n_{I_{\text{eff}}}\rangle(t)=f(\\{\beta\\},I_{1},I_{\text{eff}}),$

where $\\{\beta\\}$ are the transition coefficients of the $SIR$-like process,
and $f$ is a suitable function depending on the system we are considering. As
for the $SIR$ and $SIRS$ models, by assuming $I_{1}\approx N$ for $t\approx
0$, the critical behavior is defined by the system

$\Big{\\{}\begin{array}[]{l}\partial_{t}\langle
n_{I_{\text{eff}}}\rangle=f(\\{\beta\\},I_{\text{eff}})\\\ \partial_{t}\langle
n_{M}\rangle=\epsilon\langle n_{j}\rangle\end{array},$ (78)

where we are assuming that $I_{j}\in I_{\text{eff}}$ is the only compartment
linked to the termination $I_{M}$ by a spontaneous transition with per capita
rate $\epsilon>0$. As we will show below, these arguments can be generalized
for system with several terminations. As for the $SIRS$ model, by adding a re-
filling reaction
$\omega(\bm{n}-\mathbb{I}_{M}+\mathbb{I}_{I_{1}},\bm{n})=\sigma I_{M}$, the
system (78) becomes

$\Big{\\{}\begin{array}[]{l}\partial_{t}\langle
n_{I_{\text{eff}}}\rangle=f(\\{\beta\\},I_{\text{eff}})\\\ \partial_{t}\langle
I_{M}\rangle=\epsilon\langle I_{j}\rangle-\sigma\langle
I_{M}\rangle\end{array}.$

Since the phase transition occurs when the $I_{\text{eff}}$ becomes active,
the critical condition is independent of $I_{M}$ and it is determined by the
first equation. Thus, we conclude that the critical condition for such system
is the same of the corresponding closed diagram. We can generalize these
arguments to more complex systems by imposing only two conditions

  1. 1.

each extremity of the open diagram is linked to the main process by a
spontaneous transition;

  2. 2.

the closed diagram is given by adding spontaneous transitions from the
extremity to the initial empty state.

The first condition ensures the inner part of the process is independent of
the population numbers of the extremities 222For example, with pair
interactions, $f(\\{\beta\\})$ would depend on the population of the ends..
The second condition guarantees to have the same diagrammatic structure of the
$SIRS$ model.

##  Appendix F Further examples of models studied through the dimension
reduction

With the aim of showing the depth and breadth of our formalism, we discuss in
this Appendix section its application to a variety of models not treated in
the main text.

###  F.1 The standard $SEIRS$ model

We start with the usual $SEIRS$ model, where the spread of a disease is
assumed to occur by mean by of 4 compartments. Individuals can be in one of
the susceptible $S$, incubation $E$, infected $I$ and recovered $R$ states.
The model is defined by the following transition rates:

$\displaystyle S+I$ $\displaystyle\xrightarrow{\beta/N}I+E,$ $\displaystyle E$
$\displaystyle\xrightarrow{\lambda}I,$ $\displaystyle I$
$\displaystyle\xrightarrow{\gamma}R,$ $\displaystyle R$
$\displaystyle\xrightarrow{\sigma}S.$

In this case we are considering a system with a fixed total population $N$,
and we consider $n_{S}(0)\approx N$. Therefore, according to the notation
introduced above, $S=I_{1}$. The Doi-Peliti Hamiltonian for this model is

$\displaystyle\mathcal{H}$
$\displaystyle=\frac{\beta}{N}(a_{E}^{\dagger}-a_{S}^{\dagger})a_{S}a_{I}^{\dagger}a_{I}+\lambda(a_{I}^{\dagger}-a_{E}^{\dagger})a_{E}+$
$\displaystyle\qquad\qquad+\gamma(a_{R}^{\dagger}-a_{I}^{\dagger})a_{I}+\sigma(a_{S}^{\dagger}-a_{R}^{\dagger})a_{R},$

and the mean-field trajectories are

$\displaystyle\partial_{t}\langle n_{S}\rangle=-\frac{\beta}{N}\langle
n_{S}\rangle\langle n_{I}\rangle+\sigma\langle n_{R}\rangle,$
$\displaystyle\partial_{t}\langle n_{E}\rangle=\frac{\beta}{N}\langle
n_{S}\rangle\langle n_{I}\rangle-\lambda\langle n_{E}\rangle,$ (79)
$\displaystyle\partial_{t}\langle n_{I}\rangle=\lambda\langle
n_{E}\rangle-\gamma\langle n_{I}\rangle,$ $\displaystyle\partial_{t}\langle
n_{R}\rangle=\gamma\langle n_{I}\rangle-\sigma\langle n_{R}\rangle.$

By setting the last two equations equal to zero, and by taking into account
the conservation of total (normalized) population
$1=x_{S}q_{S}+q_{E}+q_{I}+q_{R}$, we find the following parameterization

$\displaystyle x_{E}=x_{I}=x_{R}=1\,\,,\quad
q_{E}=-\frac{\beta\sigma(x_{S}q_{S}-1)q_{S}}{\gamma\lambda+\lambda\sigma+\beta\sigma
q_{S}},$ $\displaystyle
q_{I}=\frac{\lambda\sigma(x_{S}q_{S}-1)}{\gamma\lambda+\lambda\sigma+\beta\sigma
q_{S}}\,\,,\quad
q_{R}=\frac{\lambda\gamma(x_{S}q_{S}-1)}{\gamma\lambda+\lambda\sigma+\beta\sigma
q_{S}}.$

By putting these expressions into the DPH, we obtain the reduced Hamiltonian
in the CH representation

$\mathcal{H}^{\star}=\frac{(1-\eta)e^{-\theta}\left(e^{\theta}-1\right)\lambda\sigma\left(\gamma
e^{\theta}-\beta\eta\right)}{\gamma(\lambda+\sigma)+\lambda\sigma}.$

Three zero-energy lines are found, namely,

$\eta=1\,,\quad\theta=0\,,\quad\theta=\Theta(\eta)=\log{\left(\eta\beta/\gamma\right)}.$

The absorbing phase, where the disease gets extincted, corresponds to $(1,0)$
on the plane $(\eta,\theta)$. By plotting these lines, we observe in Fig. (A2)
the well-known phase transition for the $SEIRS$ model at $\beta=\gamma$. Since
no new fixed points appear, we identify it as a second-order transition, as
expected. We can choose different ways to define the parametriazation, that is
we have several possible $\Lambda$ with different reduced DPHs. In this
particular case, the zero energy phase portraits of the different Hamiltoninas
are identical. In the next section we will show an example in which the
degeneration of $\Lambda$ leads to to have different phase portraits.

(a) $\beta<\gamma$

(b) $\beta=\gamma$

(c) $\beta>\gamma$

Figure A2: Identifying phase transitions in the $SEIRS$ model. Panels 2(a),
2(b) and 2(c) show the phase portrait of the reduced system by varying
$\beta/\gamma$. The orange curves represent the deterministic trajectories,
while the blue ones are the non-trivial energy lines. Black and red dots
represent the fixed points on the mean-field manifold. In particular the red
ones are the attractors of the dynamics. At the critical point,
$\beta=\gamma$, the attractor changes smootly and we observe a second-order
phase transition.

###  F.2 $SEIR$ model with two latent categories

Now we study the epidemiological model discussed in [15] by means of the next-
generation matrix method [14]. This model is an extension of the $SEIR$, where
two latent categories $E_{1}$ and $E_{2}$ are considered. The transition rates
are

$\displaystyle S,\,E_{1,2}$ $\displaystyle,\,I,\,R\xrightarrow{b}\emptyset,$
$\displaystyle\emptyset$ $\displaystyle\xrightarrow{bN}S,$ $\displaystyle S+I$
$\displaystyle\xrightarrow{\beta/N}E_{1,2}+I,$ $\displaystyle E_{1,2}$
$\displaystyle\xrightarrow{\nu_{1,2}}I,$ $\displaystyle I$
$\displaystyle\xrightarrow{\gamma}R,$

where $p_{1}=p$ and $p_{2}=1-p$. Its DPH reads

$\displaystyle\mathcal{H}=bN(a_{S}^{\dagger}-1)+b\Big{[}(1-a_{S}^{\dagger})a_{S}+(1-a_{E_{1}}^{\dagger})a_{E_{1}}+(1-a_{E_{2}}^{\dagger})a_{E_{2}}+(1-a_{I}^{\dagger})a_{I}+(1-a_{R}^{\dagger})a_{R}\Big{]}$
$\displaystyle\,\,+p\beta\frac{(a_{E_{1}}^{\dagger}-a_{S}^{\dagger})a^{\dagger}_{I}a_{I}a_{S}}{N}+\beta(1-p)\frac{(a_{E_{2}}^{\dagger}-a_{S}^{\dagger})a_{I}^{\dagger}a_{I}a_{S}}{N}+\nu_{1}(a_{I}^{\dagger}-a_{E_{1}}^{\dagger})a_{E_{1}}+\nu_{2}(a_{I}^{\dagger}-a_{E_{2}}^{\dagger})a_{E_{2}}+\gamma(a_{R}^{\dagger}-a_{I}^{\dagger})a_{I}$
(80)

By means of the Hamilton’s equations related to (80), the differential
equations for the deterministic trajectories are

$\displaystyle\partial_{t}\langle n_{S}\rangle=-\beta\frac{\langle
n_{S}\rangle\langle n_{I}\rangle}{N}+Nb-\langle n_{S}\rangle b$
$\displaystyle\partial_{t}\langle n_{E_{1,2}}\rangle=p_{1,2}\beta\frac{\langle
n_{S}\rangle\langle n_{I}\rangle}{N}-(b+\nu_{1,2})\langle n_{E_{1,2}}\rangle$
$\displaystyle\partial_{t}\langle n_{I}\rangle=-(\gamma+b)\langle
n_{I}\rangle+\nu_{1}\langle n_{E_{1}}\rangle+\nu_{2}\langle n_{E_{2}}\rangle$
$\displaystyle\partial_{t}\langle n_{R}\rangle=\gamma\langle
n_{I}\rangle-b\langle n_{R}\rangle$

By setting the last three equations equal to 0, and taking into account
conservation of number of particles $N$, we find the following transformations
(renormalized with respect to $N$):

$\displaystyle q_{E_{1}}$
$\displaystyle=-\frac{b(x_{S}q_{S}-1)\left[(\gamma+b)(b+\nu_{2})+\beta\nu_{2}(p-1)q_{S}\right]}{(\gamma+b)(b+\nu_{1})(b+\nu_{2})-\beta
b(p-1)q_{S}(\nu_{1}-\nu_{2})}$ $\displaystyle q_{E_{2}}$
$\displaystyle=\frac{\beta
b\nu_{1}(p-1)q_{S}(x_{S}q_{S}-1)}{(\gamma+b)(b+\nu_{1})(b+\nu_{2})-\beta
b(p-1)q_{S}(\nu_{1}-\nu_{2})}$ $\displaystyle q_{I}$
$\displaystyle=-\frac{b\nu_{1}(b+\nu_{2})(x_{S}q_{S}-1)}{(\gamma+b)(b+\nu_{1})(b+\nu_{2})-\beta
b(p-1)q_{S}(\nu_{1}-\nu_{2})}$ $\displaystyle q_{R}$
$\displaystyle=-\frac{\gamma\nu_{1}(b+\nu_{2})(x_{S}q_{S}-1)}{(\gamma+b)(b+\nu_{1})(b+\nu_{2})-\beta
b(p-1)q_{S}(\nu_{1}-\nu_{2})}$

These equations define the map $\Lambda$ from which we derive the reduced
Hamiltonian

$\displaystyle\mathcal{H}^{\star}=b\Bigg{[}\eta\left(e^{-\theta}-1\right)+e^{\theta}-1+$
$\displaystyle\qquad-\frac{\beta\eta\left(e^{\theta}-1\right)\nu_{1}(\eta-1)(b+\nu_{2})}{\beta
b(p-1)\eta(\nu_{1}-\nu_{2})-e^{\theta}(\gamma+b)(b+\nu_{1})(b+\nu_{2})}\Bigg{]}$

By solving Eqs. (19), we see that two different hyperbolic points overlap for

$\tilde{\beta}=\frac{(\gamma+b)(b+\nu_{1})(b+\nu_{2})}{b\nu_{2}+\nu_{1}\nu_{2}+b\nu_{1}p-b\nu_{2}p},$
(81)

which is the same result obtained by the next-generation matrix [15]. The
phase portraits are reported in Fig. (A3). From Figs. 3(a), 3(b) and 3(c), we
see that at the critical condition no new hyperbolic fixed points appear,
therefore the phase transition is of the second-order as verified by the
numerical simulations (Fig. 3(d)).

(a) $\beta<\tilde{\beta}$

(b) $\beta=\tilde{\beta}$

(c) $\beta>\tilde{\beta}$

(d)

Figure A3: Phase transition for the $SEIR$ with two latent categories. Panels
3(a), 3(b) and 3(c) show the phase portrait of the reduced system, where the
blue line is the mean-field manifold and the orange ones are the non-trivial
zero-energy levels. The dashed lines represent the different
parameterizations. Black and red dots represent the fixed points on the mean-
field manifold. In particular, the red ones are the attractors of the
dynamics. Panel 3(d) shows the dependence of the order parameter as a function
of $\beta$, from numerical simulation. We have considered: $b=0.2$,
$\gamma=0.3$, $p=0.9$, $\nu_{1}=0.4$ and $\nu_{2}=0.1$. Therefore, the phase
transition occurs at $\tilde{\beta}=0.7895$.

As discussed in the main text, since we have the population conservation as
first integral, there is an arbitrariness in the derivation of $\Lambda$. In
this case, in addition to the parametrization presented before, we have
studied two others $\Lambda_{2}$ and $\Lambda_{3}$. Choosing as fist integrals
the conservation of $N$, the second, the fourth and the fifth deterministic
equations we obtain $\Lambda_{2}$, while by taking the $N$ conservation with
the second, the third and the fifth equations we have $\Lambda_{3}$. The phase
portrait of the corresponding reduced Hamiltonians is represented in Fig.
(A3). We have analyzed the long-time behavior of the different Langevin
dynamics (25). For each equation, we have performed a simulation in the active
phase and compared the stationary distribution for $n_{S}$. The results are
shown in Fig. (A4). Distributions show differences as expected, however they
have the identical average and provide a qualitative description of the
stochastic environment of the system. Further studies are needed to research a
criterion for selecting only one parameterization over the others.

Figure A4: Numerical asymptotic distribution of the susceptible individuals
for the $SEIR$ model with two latent categories. The different colors
represent the three different parametrizations: $\Lambda$, $\Lambda_{2}$ and
$\Lambda_{3}$, We have considered $N=10^{4}$, $\beta=1.11$, $b=0.2$,
$\nu_{1}=0.4$, $\nu_{2}=0.1$, $p=0.9$ and $\gamma=0.3$. In this case, the
critical point is $\tilde{\beta}=0.7895$.

###  F.3 Tax evasion model

We analyse with our formalism a simple model for tax evasion that undergoes
first-order phase transitions [63]. The models includes 3 compartments: honest
individuals $H$, susceptible ones $S$ and evaders $E$. The dynamics are
described by the following elementary reactions,

$\displaystyle E+H$ $\displaystyle\xrightarrow{\lambda/N}E+S,$ $\displaystyle
S$ $\displaystyle\xrightarrow{\alpha}E,$ $\displaystyle S+H$
$\displaystyle\xrightarrow{\epsilon/N}2H,$ $\displaystyle H+E$
$\displaystyle\xrightarrow{\delta/N}2H,$ $\displaystyle E$
$\displaystyle\xrightarrow{\beta}H.$

The DPH is given by

$\displaystyle\mathcal{H}=\frac{\epsilon}{N}a_{H}^{\dagger}(a_{H}^{\dagger}-a_{S}^{\dagger})a_{H}a_{S}+\frac{\delta}{N}a_{H}^{\dagger}(a_{H}^{\dagger}-a_{E}^{\dagger})a_{H}a_{E}+$
$\displaystyle\frac{\lambda}{N}a_{E}^{\dagger}(a_{S}^{\dagger}-a_{H}^{\dagger})a_{H}a_{E}+\alpha(a_{E}^{\dagger}-a_{S}^{\dagger})a_{S}+\beta(a_{H}^{\dagger}-a_{E}^{\dagger})a_{E},$

and the mean-field equations read

$\displaystyle\partial_{t}\langle n_{H}\rangle$
$\displaystyle=\frac{\epsilon}{N}\langle n_{H}\rangle\langle
n_{S}\rangle+\beta\langle n_{E}\rangle+\frac{\delta}{N}\langle
n_{H}\rangle\langle n_{E}\rangle$
$\displaystyle\qquad\qquad\qquad-\frac{\lambda}{N}\langle n_{H}\rangle\langle
n_{E}\rangle$ $\displaystyle\partial_{t}\langle n_{S}\rangle$
$\displaystyle=-\alpha\langle n_{S}\rangle-\frac{\epsilon}{N}\langle
n_{H}\rangle\langle n_{S}\rangle+\frac{\lambda}{N}\langle n_{H}\rangle\langle
n_{E}\rangle$ $\displaystyle\partial_{t}\langle n_{E}\rangle$
$\displaystyle=\alpha\langle n_{S}\rangle-\beta\langle
n_{E}\rangle-\frac{\delta}{N}\langle n_{E}\rangle\langle n_{H}\rangle.$

In this model, $N$ is fixed and we consider $\rho=1-\lim_{t\to\infty}\langle
n_{H}\rangle/N$ as the order parameter. After re-scaling the system size to
$N=1$, we take as first integrals the last mean-field equation and the
population conservation, obtaining the parametrization

$q_{S}=-\frac{\lambda
q_{H}(x_{H}q_{H}-1)}{\alpha+q_{H}(\lambda+\epsilon)}\,\,,\qquad
q_{E}=\frac{(\alpha+q_{S}\epsilon)(1-x_{H}q_{H})}{\alpha+q_{H}(\epsilon+\lambda)}.$

The corresponding reduced Hamiltonian is:

$\displaystyle\mathcal{H}^{\star}$
$\displaystyle=-\frac{(\eta-1)e^{-\theta}\left(e^{\theta}-1\right)}{\alpha
e^{\theta}+\eta(\lambda+\epsilon)}\Big{(}\alpha
e^{2\theta}(\beta+\delta\eta)+$ $\displaystyle\qquad\eta
e^{\theta}(-\alpha\lambda+\beta\epsilon+\delta\eta\epsilon+\eta\lambda\epsilon)-\eta^{2}\lambda\epsilon\Big{)}.$

The zero-energy curves intersecting the mean-field line are

$\eta\equiv
1\,,\qquad\theta=\Theta(\eta)=\log\left(\eta\frac{\sqrt{(-\alpha\lambda+\beta\epsilon+\delta\eta\epsilon+\eta\lambda\epsilon)^{2}+4\alpha\lambda\epsilon(\beta+\delta\eta)}+\alpha\lambda-\beta\epsilon-\delta\eta\epsilon-\eta\lambda\epsilon}{2\alpha(\beta+\delta\eta)}\right).$

We denote by $F_{-}$ and $F_{+}$ the intersection of $\Theta$ with the mean-
field line, where $\eta_{F_{-}}\leq\eta_{F_{+}}$, and $F_{1}=(1,0)$ is the
usual absorbing state (see Fig. A5). By considering $\lambda$ as control
parameter, we find that $\Theta(1)=0$ for $\lambda$ equal to

$\tilde{\lambda}:=\frac{(\alpha+\epsilon)(\beta+\delta)}{\alpha}.$ (82)

For $\lambda<\tilde{\lambda}$, as reported in Figs. 5(a) and 5(d), the system
converges to $F_{1}$ and there are no phase transition. Moreover, at the
critical point, for $\epsilon$ equal to

$\tilde{\epsilon}:=\frac{\alpha\beta}{\delta},$

$\Theta$ is tangent to the mean-field line and $F_{+}=F_{-}=(1,0)$, so we find
a tricritical point. In particular, for $\epsilon<\tilde{\epsilon}$, at the
critical point we find $F_{-}=(1,0)$ (Panel 5(b)). By studying the
eigenvalues, we observe that $F_{-}$ becomes the attractor (Panel 5(c)) and a
second-order phase transition occurs. Otherwise, if
$\epsilon>\tilde{\epsilon}$, at $\lambda=\tilde{\lambda}$ we have
$F_{+}=(1,0)$ as shown in Fig. 5(e) and for $\lambda>\tilde{\lambda}$ the
dynamics jumps to the new attractor $F_{-}$ (Fig. 5(f)). Therefore, there is
an abrupt changing of the asymptotic behavior, that is a first-order
transition. A picture of this behaviour is reported in Fig. A5 and in Fig. A6
we show a numerical simulation for both the first- and second- order phase
transition. We point out that, while the critical value (82) corresponds to
the result predicted in [63], with our method we find a continuous critical
behaviour from low $\epsilon$ that had not been investigated before.

(a) $\lambda<\tilde{\lambda},\,\epsilon<\tilde{\epsilon}$

(b) $\lambda=\tilde{\lambda},\,\epsilon<\tilde{\epsilon}$

(c) $\lambda>\tilde{\lambda},\,\epsilon<\tilde{\epsilon}$

(d) $\lambda<\tilde{\lambda},\,\epsilon>\tilde{\epsilon}$

(e) $\lambda=\tilde{\lambda},\,\epsilon>\tilde{\epsilon}$

(f) $\lambda>\tilde{\lambda},\,\epsilon>\tilde{\epsilon}$

Figure A5: Phase portrait for the taxes evasion model. The orange lines
represent the mean-field manifold, while the blue ones are the other non zero-
energy lines. $F_{1}$ and $F_{\pm}$ are the system fixed points. Red dots are
the attractors of the deterministic dynamics. In the upper Panels 5(a), 5(b)
and 5(c), we report the second-order phase transition behavior: with
$\epsilon<\tilde{\epsilon}$, we see that the intersection $F_{1}$ overlaps
with $F_{1}$ and becomes the new attractor. The asymptotic state changes
smoothly and we have a continuous phase transitions. Otherwise, in the lower
Panels 5(d), 5(e) and 5(f), when $\lambda=\tilde{\lambda}$, we have
$F_{1}=F_{+}$ and the dynamics varies in a discontinuous way converging in
$F_{-}$. So we assist to a first-order phase transitions. Figure A6: Numerical
asymptotic distribution of the honest individuals. The blue line represents
the second order phase transition, while the orange one is the first order
transition. For both the cases we have considered: $\beta=0.2$, $\alpha=0.3$,
$\gamma=0.1$, thus $\tilde{\epsilon}=0.6$. For the first order transition
$\epsilon=0.8$, hence $\tilde{\lambda}=1.1$. For the continuous transition
instead $\epsilon=0.3$ and $\tilde{\lambda}=0.6$.

###  F.4 Levin’s model

Let us now study another model of ecological interest widely treated in the
literature: the Levin’s model. This model is characterized by two types of
processes, namely, one of population proliferation in an empty patch, and one
of extinction. We focus on systems containing two species $A$ and $B$ whose
dynamics are defined by the rates

$\displaystyle\omega(\bm{n}+\mathbb{I}_{A}-\mathbb{I}_{\emptyset},\bm{n})=\frac{c}{N}n_{A}(n_{\emptyset}-D)$
(83)
$\displaystyle\omega(\bm{n}+\mathbb{I}_{B}-\mathbb{I}_{\emptyset},\bm{n})=\frac{\beta}{N^{2}}n_{A}n_{B}(n_{\emptyset}-D)$
$\displaystyle\omega(\bm{n}-\mathbb{I}_{A,B}+\mathbb{I}_{\emptyset},\bm{n})=n_{A,B}\epsilon$

where $n_{\emptyset}$ represents the number of non-occupied patches. Moreover,
we are including $D$ patches that have been damaged and can no longer be
occupied. The population $A$ has an independent growth, while the dynamics of
$B$ depends on $A$. From these rates, we see that $N\geq n_{\emptyset}\geq D$
to have a well defined dynamics.

The Doi-Peliti Hamiltonian reads

$\displaystyle\mathcal{H}$
$\displaystyle=\frac{cD}{N}a_{A}^{\dagger}\left(a_{A}^{\dagger}-a_{\emptyset}^{\dagger}\right)a_{A}a_{\emptyset}-\frac{cD}{N}a_{A}^{\dagger}\left[(a_{\emptyset})^{-1}-1\right]a_{A}+\frac{\beta}{N}\left(a_{B}^{\dagger}-a_{\emptyset}^{\dagger}\right)a_{A}^{\dagger}a_{B}^{\dagger}a_{A}a_{B}a_{\emptyset}$
$\displaystyle\,\qquad\qquad\qquad-\frac{\beta
D}{N}\left[a_{B}^{\dagger}(a^{\dagger}_{\emptyset})^{-1}-1\right]a^{\dagger}_{A}a^{\dagger}_{B}a_{A}a_{B}+\epsilon\left(a_{\emptyset}^{\dagger}-a_{A}^{\dagger}\right)a_{A}+\epsilon\left(a_{\emptyset}^{\dagger}-a_{B}^{\dagger}\right)a_{B},$

and the mean-field equations are

$\displaystyle\partial_{t}\langle n_{A}\rangle=\frac{c}{N}\langle
n_{A}\rangle\left(\langle n_{\emptyset}\rangle-D\right)-\epsilon\langle
n_{A}\rangle$ $\displaystyle\partial\langle
n_{B}\rangle=\frac{\beta}{N^{2}}\langle n_{A}\rangle\langle
n_{B}\rangle\left(\langle n_{\emptyset}\rangle-D\right)-\epsilon\langle
n_{B}\rangle$ $\displaystyle\partial\langle
n_{\emptyset}\rangle=-\partial_{t}\langle n_{A}\rangle-\partial_{t}\langle
n_{B}\rangle$

We choose $\langle n_{A}\rangle/N$ as the order parameter and $D$ as the
control variable. The solution of (19) yields two critical points

A) $\displaystyle\frac{\langle n_{A}\rangle}{N}=0,\,\,\,\frac{\langle
n_{B}\rangle}{N}=0,\,\,\,\frac{\langle
n_{\emptyset}\rangle}{N}=1,\,\,\,\tilde{D}_{1}=N\left(1-\frac{\epsilon}{c}\right),$
B) $\displaystyle\frac{\langle
n_{A}\rangle}{N}=\frac{c}{\beta},\,\,\,\frac{\langle
n_{B}\rangle}{N}=0,\,\,\,\frac{\langle
n_{\emptyset}\rangle}{N}=\left(1-\frac{c}{\beta}\right),$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\,\,\,\tilde{D}_{2}=N\left(1-\frac{c}{\beta}-\frac{\epsilon}{c}\right).$

Thus $\tilde{D}_{1}>\tilde{D}_{2}$, and we expect to observe two phase
transitions. By solving Eq. (18) with the second mean-field equation and the
population conservation as constraints, we obtain two parametrizations
$\Lambda(x_{A},q_{A})=(x_{A},x_{B},x_{\emptyset},q_{A},q_{B},q_{\emptyset})$,
as in Appendix IV.4. They are

$\displaystyle\Lambda^{(A)}$ $\displaystyle=$
$\displaystyle\,\,\left(x_{A},1,1,q_{A},N-D-q_{A}x_{A}-\frac{\epsilon
N^{2}}{q_{A}\beta},D+\frac{\epsilon}{q_{A}\beta}\right),$
$\displaystyle\Lambda^{(B)}$
$\displaystyle=\left(x_{A},1,1,q_{A},0,N-x_{A}q_{A}\right).$

We denote by $\mathcal{H}^{\star(A)}$ and $\mathcal{H}^{\star(B)}$
respectively the reduced Hamiltonians for $\Lambda^{(A)}$ and $\Lambda^{(B)}$,
that read

$\displaystyle\mathcal{H}^{\star(A)}$
$\displaystyle=\frac{e^{-\theta}\left(e^{\theta}-1\right)\epsilon\left(ce^{2\theta}-\beta\eta\right)}{\beta}$
$\displaystyle\mathcal{H}^{\star(B)}$
$\displaystyle=-e^{-\theta}\left(e^{\theta}-1\right)\eta\left(ce^{\theta}(d+\eta-1)+\epsilon\right),$

where we have re-normalized the population with respect to $N$ and we have
introduced $d=D/N$. By analysing their zero-energy lines with the attractors
of the dynamics, we note that the system is described by
$\mathcal{H}^{\star(B)}$ for $D<\tilde{D}_{1}$, and by
$\mathcal{H}^{\star(A)}$ otherwise. A sketch of the phase portraits is shown
in panels 7(a), 7(b) and 7(c).

A numerical simulation for the system is reported in panel 7(d). There, we can
observe 3 regimes, separated by the two phase transitions found above. The
first regime, for $D<\tilde{D}_{2}$, is characterize by the coexistence of
both populations. In the second one, for $\tilde{D}_{2}<D<\tilde{D}_{1}$, only
population $A$ survives. For $D>\tilde{D}_{1}$, the fragmentation is too
strong and there is complete extinction of all species. In both cases we have
a second-order phase transition.

(a) $D<\tilde{D}_{2}$

(b) $\tilde{D}_{2}<D<\tilde{D}_{1}$

(c) $D>\tilde{D}_{1}$

(d)

Figure A7: Phase transition for the Levin’s model. Panels 7(a), 7(b) and 7(c)
show the phase portraits of the reduced system, where the orange lines are the
mean-field manifold, the red lines are the zero-energy line of
$\mathcal{H}^{\star(A)}$ and the blue ones are related to
$\mathcal{H}^{\star(B)}$. Red dots represent the fixed points on the mean-
field manifold. Panel 7(d) is a numerical simulation of the model with
$\beta=1$, $c=0.7$ and $\epsilon=0.1$, hence $\tilde{d}_{1}=0.157$ and
$\tilde{d}_{2}=0.857$.

## References

  * Freeman  _et al._ [2004] L. Freeman _et al._ , The development of social network analysis, A Study in the Sociology of Science 1, 159 (2004).
  * Montoya  _et al._ [2006] J. M. Montoya, S. L. Pimm, and R. V. Solé, Ecological networks and their fragility, Nature 442, 259 (2006).
  * Kossinets and Watts [2006] G. Kossinets and D. J. Watts, Empirical analysis of an evolving social network, Science 311, 88 (2006).
  * Fath  _et al._ [2007] B. D. Fath, U. M. Scharler, R. E. Ulanowicz, and B. Hannon, Ecological network analysis: network construction, Ecological Modelling 208, 49 (2007).
  * Castellano  _et al._ [2009] C. Castellano, S. Fortunato, and V. Loreto, Statistical physics of social dynamics, Reviews of Modern Physics 81, 591 (2009).