# Nonequilibrium Thermodynamics of Non-Ideal Reaction-Diffusion Systems:
Implications for Active Self-Organization

Francesco Avanzini<EMAIL_ADDRESS>Department of Chemical
Sciences, University of Padova, Via F. Marzolo, 1, I-35131 Padova, Italy
Timur Aslyamov<EMAIL_ADDRESS>Department of Physics and Materials
Science, University of Luxembourg, L-1511 Luxembourg City, Luxembourg  Étienne
Fodor<EMAIL_ADDRESS>Department of Physics and Materials Science,
University of Luxembourg, L-1511 Luxembourg City, Luxembourg  Massimiliano
Esposito<EMAIL_ADDRESS>Department of Physics and Materials
Science, University of Luxembourg, L-1511 Luxembourg City, Luxembourg

###### Abstract

We develop a framework describing the dynamics and thermodynamics of open non-
ideal reaction-diffusion systems, which embodies Flory-Huggins theories of
mixtures and chemical reaction network theories. Our theory elucidates the
mechanisms underpinning the emergence of self-organized dissipative structures
in these systems. It evaluates the dissipation needed to sustain and control
them, discriminating the contributions from each reaction and diffusion
process with spatial resolution. It also reveals the role of the reaction
network in powering and shaping these structures. We identify particular
classes of networks in which diffusion processes always equilibrate within the
structures, while dissipation occurs solely due to chemical reactions. The
spatial configurations resulting from these processes can be derived by
minimizing a kinetic potential, contrasting with the minimization of the
thermodynamic free energy in passive systems. This framework opens the way to
investigating the energetic cost of phenomena such as liquid-liquid phase
separation, coacervation, and the formation of biomolecular condensates.

##  I Introduction

Self-organization of molecules in solution can produce complex spatio-temporal
structures that range from self-assemblies to phase separations, encompassing
stationary patterns and oscillations. It results from nonlinear effects which
can be caused by a variety of mechanisms.

One possible mechanism is based on molecular interactions, typical of non-
ideal systems, which can lead to self-organization even at equilibrium. When
no external energy is provided to sustain self-organization, this mechanism is
passive and described by a dynamics exhibiting detailed balance. The
theoretical description for this mechanism is typically provided by Cahn-
Hilliard theory of spinodal decomposition Cahn and Hilliard (1958) and its
variants Bray (1994). It can be further understood from a thermodynamic
perspective through theories like that proposed by Flory and Huggins Huggins
(1941); Flory (1942) as a relaxation towards equilibrium.

Another mechanism involves multimolecular chemical reactions in open systems,
which can also lead to self-organization in the absence of interactions.
Predicted by A. M. Turing Turing (1952) and extensively studied by the
Brussels School Prigogine and Nicolis (1971); Nicolis and Prigogine (1977);
Lefever (2018), its experimental validation took nearly four decades Castets
_et al._ (1990). From a thermodynamic viewpoint, this type of self-
organization is maintained by the continuous consumption of energy Falasco
_et al._ (2018); Avanzini  _et al._ (2019, 2020), making the mechanism active
with dynamics that exhibits broken detailed balance.

Recent years have seen renewed interest in the study of self-organization due
to its crucial role in biological systems, where biomolecular condensates
Banani  _et al._ (2017); Lyon  _et al._ (2021) resulting from liquid-liquid
phase separation from membraneless compartments. These condensates have been
observed to serve various functions, including storage of molecules Klosin
_et al._ (2020), enhancement of reaction rates by concentrating enzymes and
substrates Castellana  _et al._ (2014); Peeples and Rosen (2020), and
maintenance of macromolecule folding states Frottin  _et al._ (2019).

Furthermore, the interplay between passive Berry  _et al._ (2018) and active
Weber  _et al._ (2019); Zwicker (2022) mechanisms in the formation and control
of these condensates has gained significant attention. Indeed, their
nonequilibrium nature is not easily discernible by merely examining the self-
organized spatial structures. Instead, it requires examining both diffusion
and reaction fluxes, and how these fluxes relate to energy consumption. This
critical feature, known as thermodynamic consistency, is often overlooked in
models that only aim to reproduce the concentration profiles of self-organized
structures, as seen in previous works where chemical reactions Demarchi  _et
al._ (2023); Glotzer  _et al._ (1995) or diffusion processes Brauns  _et al._
(2020); Haas and Goldstein (2021) were modeled in a thermodynamically
inconsistent manner.

Understanding the control of self-organized structures by chemical reactions
and the resulting dissipation of energy is crucial, especially in the field of
biology. Previous studies have explored the regulation of such structures
through simple unimolecular reactions Kirschbaum and Zwicker (2021); Bauermann
_et al._ (2022) or examined their stability in multimolecular reactions
Aslyamov  _et al._ (2023), yet a comprehensive framework to analyze the
nonequilibrium aspects of self-organization is still lacking. This work
addresses this gap by highlighting the critical role of the reaction network
topology in defining the relationship between self-organized structures and
their thermodynamics.

We develop a general thermodynamic theory for open non-ideal reaction-
diffusion systems by integrating the thermodynamic theories of both non-ideal
chemical reaction networks Avanzini  _et al._ (2021) and ideal reaction-
diffusion systems Falasco  _et al._ (2018); Avanzini  _et al._ (2019). This
unified approach enables us to investigate the energetics of self-organized
structures and the interplay between passive and active mechanisms. We examine
various classes of reaction networks. For detailed balanced ones (open or
closed), we show that the dynamics always relaxes to equilibrium while
minimizing a suitable thermodynamic potential which plays the role of a
Lyapunov function that depends only on the energetics and not on the kinetics.
In these passive systems, the resulting self-organized structures do not
dissipate. For pseudo-detailed balanced reaction networks, we show that the
dynamics can be mapped into a detailed-balanced dynamics minimizing a kinetic
potential that we identify and that plays the role of a Lyapunov function.
Nevertheless, these systems are active and the resulting self-organized
structures dissipate. For complex balanced reaction networks (characterized by
Arrhenius-like reaction fluxes), such a mapping does not exist, but the
dynamics still minimizes a kinetic potential that we construct and that plays
the role of a Lyapunov function. These systems are also active, and the self-
organized structures they produce also dissipate. Remarkably, even though
these systems are active, diffusion processes equilibrate at steady state in
both pseudo-detailed balanced and complex balanced networks. However, most
active systems do not fall into these classes. Nonetheless, we show that our
theoretical framework can be used to systematically analyze the energetics and
the dissipation produced by the various active processes generating the self-
organized structures.

The paper is structured as follows. We begin in Section II by outlining the
basic setup. Specifically, we define the representations of chemical reactions
in terms of chemical reaction networks and graphs of complexes, which are
crucial for the derivations of our results. In Section III, we introduce the
dynamical description of the concentration fields and explore the concept of
detailed balance in Subsection III.1. In Subsection III.2, we discuss the
conservation laws of chemical reaction networks to identify molecule
fragments, known as moieties, which are transferred between chemical species
through reactions while remaining intact. These moieties will be instrumental
in identifying the nonconservative forces that break detailed balance. We turn
to thermodynamics in Section IV. We first discuss the conditions needed for
the dynamics to be thermodynamically consistent. We then use moieties to
rewrite the second law of thermodynamics in a way that identifies the
thermodynamic forces driving the system out of equilibrium and the
thermodynamic potential acting as a Lyapunov function when forces are switched
off. In Section V, we make use of thermodynamics, along with the
representations of chemical reactions in terms of chemical reaction networks
and graphs of complexes, to derive our main findings for pseudo-detailed
balanced networks and complex balanced networks. We end our study in Section
VI by analyzing the energetics of various reaction networks using numerical
simulations, including a network that is neither pseudo-detailed balanced nor
complex balanced. Conclusions are drawn in Section VII.

##  II Basic Setup

We consider non-ideal reaction-diffusion (RD) systems in solution: $n_{s}$
interacting, reacting, and diffusing chemical species (identified by the
labels $i\in\mathcal{S}$) are mixed together with a non-reacting and abundant
species called solvent which maintains the temperature $T$ and the volume $V$
constant. We further assume that RD systems are embedded within impermeable
walls. The species $\mathcal{S}$ are either chemostatted (labeled $y\in Y$),
if exchanged with external reservoirs called chemostats (labeled $y\in Y$
too), or internal (labeled $x\in X$). We use $n_{y}$ and $n_{x}$ for the
number of chemostatted and internal species, respectively. Due to the
exchanges with the chemostats, RD systems are said to be open. Without
exchanges, RD systems would be closed.

###  II.1 Chemical Reaction Networks and Graphs of Complexes

Chemical species are interconverted via chemical reactions
$\rho\in\mathcal{R}=\\{\pm 1,\pm 2,\dots,\pm n_{r}\\}$ characterized by the
equation

$\small\sum_{x}\nu_{x,\rho}\,Z_{x}+\sum_{y}\nu_{y,\rho}\,Z_{y}\ch{<=>[\rho][-\rho]}\sum_{x}\nu_{x,-\rho}\,Z_{x}+\sum_{y}\nu_{y,-\rho}\,Z_{y}\,.$
(1)

Here, $Z_{i}$ is the chemical symbol of the species $i$, while $\nu_{i,\rho}$
(resp. $\nu_{i,-\rho}$) is its stoichiometric coefficient in reaction $\rho$
(resp. $-\rho$). The set of species $\mathcal{S}$ might include non-reacting
species (labeled for instance nr). All reactions are assumed to be reversible:
for every (forward) reaction $\rho\in\mathcal{R}$, reaction
$-\rho\in\mathcal{R}$ and denotes its backward counterpart.

On the one hand, chemical reactions (1) can be represented as a hypergraph,
known as chemical reaction network (CRN), by mapping species into nodes and
reactions into edges Klamt  _et al._ (2009). Its topology is encoded in the
$n_{s}\times n_{r}$ stoichiometric matrix $\mathbb{S}$ whose entries (for
$\rho>0$),

${S}_{i,\rho}=\nu_{i,-\rho}-\nu_{i,\rho}\,,$ (2)

specify the net variation of molecule number of species $i$ in reaction
$\rho$. By definition, ${S}_{i,-\rho}=-{S}_{i,\rho}$ and
${S}_{\text{nr},\rho}=0$ $\forall\rho$.

On the other hand, the set of chemical reactions (1) can also define the so-
called graph of complexes. The $n_{k}$ complexes (identified by the labels
$k\in\mathcal{K}$) are aggregates of internal species acting as reactants of a
reaction,

$V_{k(\rho)}=\sum_{x}\nu_{x,\rho}\,Z_{x}\,,$ (3)

and the graph of complexes is obtained using the complexes as nodes and
reactions as edges:

$V_{k(\rho)}\ch{<=>[\rho][-\rho]}V_{k(-\rho)}\,.$ (4)

Note that the same complex can be involved in different reactions. For
instance, if reactions have no internal species among their reactants, they
share the same complex (usually represented by the symbol of empty set, i.e.,
$\emptyset$, as done in the example in Subs. VI.4). The topology of the graph
of complexes is encoded in the $n_{k}\times n_{r}$ incident matrix $\partial$
whose entries (for $\rho>0$),

$\partial_{k,\rho}={\mathbb{1}}_{k,k(-\rho)}-{\mathbb{1}}_{k,k(\rho)},$ (5)

specify whether the complex $k$ is a product (positive value) or a reagent
(negative value) of reaction $\rho$ (with ${\mathbb{1}}_{k,k^{\prime}}$ being
the Kronecker delta).

The graph of complexes allows us to divide the set of reactions $\mathcal{R}$
into $2n_{e}$ disjoint lumped reactions identified by the labels
$\varepsilon\equiv(k,k^{\prime})\in\mathcal{E}$, that collect the reactions
interconverting the complex $V_{k}$ into $V_{k^{\prime}}$. Namely,

$\mathcal{R}=\bigcup_{\varepsilon}\mathcal{R}_{\varepsilon}\,,$ (6)

with

$\mathcal{R}_{\varepsilon}\equiv\big{\\{}\rho\in\mathcal{R}\text{ such that
}V_{k(\rho)}=V_{k}\text{ and }V_{k(-\rho)}=V_{k^{\prime}}\big{\\}}\,.$ (7)

Note that, if $\rho\in\mathcal{R}_{\varepsilon}$, then
$-\rho\in\mathcal{R}_{-\varepsilon}$ with $-\varepsilon\equiv(k^{\prime},k)$.
By definition, all reactions $\rho\in\mathcal{R}_{\varepsilon}$ have the same
stoichiometric coefficient (resp. matrix) for each internal species
${\nu}_{x,\rho}$ (resp. ${S}_{x,\rho}$), labeled ${\nu}_{x,\varepsilon}$
(resp. ${S}_{x,\varepsilon}$) hereafter. On the other hand, there might be two
reactions $\rho$ and $\rho^{\prime}$ in $\mathcal{R}_{\varepsilon}$ such that
${S}_{y,\rho}\neq{S}_{y,\rho^{\prime}}$.

We conclude this subsection by introducing a topological quantity, called
deficiency, which essentially quantifies the number of “hidden” cycles of CRNs
(i.e., cycles that have not a graphical representation in the graph of
complexes) and is known to underlie the emergence of complex behavior in ideal
CRNs Feinberg (1987); Craciun  _et al._ (2006); Lubensky (2010); Anderson  _et
al._ (2015); Polettini  _et al._ (2015); Marehalli Srinivas  _et al._ (2023).
A cycle is formally defined as a right-null eigenvector
$\boldsymbol{\phi}=(\dots,\phi_{\rho},\dots)$ of the substoichiometric matrix
for the internal species $\mathbb{S}_{X}$:

$\sum_{\rho>0}{S}_{x,\rho}\phi_{\rho}{}=0.$ (8)

It thus represents a sequence of reactions (by specifying the number of times
each reaction occurs) that upon completion leaves the molecule number of the
internal species unchanged. Some cycles can be represented as loops in the
graph of complexes. Indeed, the substoichiometric matrix for the internal
species $\mathbb{S}_{X}$ can be expressed in terms of the incidence matrix
$\partial$ using the composition matrix $\Gamma$ as

${S}_{x,\rho}=\sum_{k}\Gamma_{x,k}\,\partial_{k,\rho}\,,$ (9)

where

$\Gamma_{x,k(\rho)}=\nu_{x,\rho}$ (10)

specifies the stoichiometric coefficient of each internal species $x$ in each
complex $k$. Hence, from Eqs. (8) and (9), it follows that any right-null
eigenvector of the incidence matrix $\partial$, defining a loop in the graph
of complexes, corresponds to a cycle. However, the converse is not true and
this difference is encoded in the deficiency

$\delta\equiv\text{dim ker}(\mathbb{S}_{X})-\text{dim ker}(\partial)\,,$ (11)

where $\text{dim ker}(\bullet)$ returns the dimension of the kernel of a
matrix, i.e., the number of linearly-independent right-null eigenvectors.

###  II.2 Example

Consider a system composed of one internal species $\\{\ch{X_{1}}\\}$ and
three chemostatted species $\\{\ch{Y_{1}},\ch{Y_{2}},\ch{Y_{3}}\\}$ which are
interconverted by the following chemical reactions:

$\begin{split}\ch{Y_{1}}&\ch{<=>[$1$][$-1$]}\ch{X_{1}}\,,\\\
\ch{X_{1}}+\ch{Y_{2}}&\ch{<=>[$2$][$-2$]}2\ch{X_{1}}\,,\\\
\ch{X_{1}}+\ch{Y_{3}}&\ch{<=>[$3$][$-3$]}2\ch{X_{1}}\,.\end{split}$ (12)

The first pair of reactions $\\{\pm 1\\}$ represents the direct
interconversion of $\ch{Y_{1}}$ into $\ch{X_{1}}$, while the second (resp.
third) pair of reactions $\\{\pm 2\\}$ (resp. $\\{\pm 3\\}$) represents the
interconversion of $\ch{Y_{2}}$ (resp. $\ch{Y_{3}}$) into $\ch{X_{1}}$
promoted by $\ch{X_{1}}$ in an autocatalytic way. In Fig. 1, we illustrate the
CRN and graph of complexes representations introduced in Subs. II.1 for the
chemical reactions in Eq. (12).

Note that the substoichiometric matrix $\mathbb{S}_{X}$ in Fig. 1 admits two
(linearly-independent) cycles, namely,

$\boldsymbol{\phi}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr}}}}\right)$}}\quad\text{ and
}\quad\boldsymbol{\phi}^{\prime}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern
0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr}}}}\right)$}}\,.$ (13)

The former is also a right-null eigenvector of the incidence matrix
corresponding to the loop

$V_{1}\ch{->[2][\end{equation}]}$

V_2 -¿[$-3$][

$]$

V_1 , while the latter is not. Hence, $\delta=1$.

Figure 1: Different representations of the chemical reactions in Eq. (12). (a)
CRN representation in terms of a hypergraph where the species
$\\{\ch{X_{1}},\ch{Y_{1}},\ch{Y_{2}},\ch{Y_{3}}\\}$ are mapped into nodes,
while the reactions $\rho\in\\{1,2,3\\}$ are mapped into edges. The green
double-headed arrows represent the exchange processes of the chemostatted
species. The topological properties of the CRN are encoded in its
stoichiometric matrix $\mathbb{S}$, where the red dashed (resp. green dotted)
line highlights the substoichiometric matrix for the internal (resp.
chemostatted) species $\mathbb{S}_{X}$ (resp. $\mathbb{S}_{Y}$). (b) Graph of
complexes representation where the complexes $\emptyset=\ch{Y_{1}}$,
$V_{1}=\ch{X_{1}}$, and $V_{2}=\ch{2X_{1}}$ are mapped into nodes, while the
reactions $\rho\in\\{1,2,3\\}$ are mapped into edges, and the exchange
processes are not represented. The topological properties of the graph of
complexes are encoded in its incidence matrix $\partial$. The
substoichiometric matrix for the internal species $\mathbb{S}_{X}$ is related
to the incidence matrix $\partial$ by the composition matrix $\Gamma$. The
reactions $\rho=2$ and $\rho=3$ (highlighted by the blue dashed line)
interconvert the same pair of complexes, i.e., $V_{1}$ and $V_{2}$, and thus
correspond to a single lumped reaction. (c) Simplified representation of the
graph of complexes where the lumped reactions
$\varepsilon_{1}=(\emptyset,V_{1})$ and $\varepsilon_{2}=(V_{1},V_{2})$ are
mapped into edges instead of the reactions. Note that all reactions, as well
all lumped reactions, are assumed to be reversible even though only the
forward reactions and lumped reactions are represented.

##  III Dynamics

The abundances of the chemical species $i\in\mathcal{S}$ in every point of
space $\boldsymbol{r}\in V$ are specified by the concentration fields
$\boldsymbol{c}(\boldsymbol{r})=(\dots,c_{i}(\boldsymbol{r}),\dots)$ following
the RD equation:

$\partial_{t}c_{i}(\boldsymbol{r})=-\boldsymbol{\nabla}\cdot\boldsymbol{J}_{i}(\boldsymbol{c}(\boldsymbol{r}))+\sum_{\rho>0}{S}_{i,\rho}\,j_{\rho}(\boldsymbol{c}(\boldsymbol{r}))+I_{i}(\boldsymbol{r})\,,$
(14)

where $I_{i}(\boldsymbol{r})$ is the exchange current of species $i$ with the
corresponding chemostat (and thus $I_{x}(\boldsymbol{r})=0$
$\forall\boldsymbol{r}$); $j_{\rho}(\boldsymbol{c}(\boldsymbol{r}))$ is the
net current of reaction $\rho>0$ given by difference between the forward
$\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))\geq 0$ and backward
$\omega_{-\rho}(\boldsymbol{c}(\boldsymbol{r}))\geq 0$ reaction flux,

$j_{\rho}(\boldsymbol{c}(\boldsymbol{r}))=\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))-\omega_{-\rho}(\boldsymbol{c}(\boldsymbol{r}))\,;$
(15)

$\boldsymbol{J}_{i}(\boldsymbol{c}(\boldsymbol{r}))$ is the diffusion current
of species $i$. For sake of compactness, we omit the time-dependence of any
function (see for instance the concentration fields in Eq. (14)).

Remark. Note that we implicitly assumed that the concentration fields of the
chemical species are large enough that their dynamics becomes deterministic
Falasco and Esposito (2024).

###  III.1 Steady State

The steady-state solutions of the RD equation (14), i.e.,
$\partial_{t}\boldsymbol{c}(\boldsymbol{r})=0$, if they exist, can be
classified as either equilibrium or nonequilibrium steady states. Equilibrium
steady states $\boldsymbol{c}_{\mathrm{eq}}(\boldsymbol{r})$ satisfy

$j_{\rho}(\boldsymbol{c}_{\mathrm{eq}}(\boldsymbol{r}))=0\text{ }\text{
}\text{ }\text{ and }\text{ }\text{ }\text{
}\boldsymbol{J}_{i}(\boldsymbol{c}_{\mathrm{eq}}(\boldsymbol{r}))=0$ (16)

for all $\rho\in\mathcal{R}$, $i\in\mathcal{S}$, and $\boldsymbol{r}\in V$.
When the RD dynamics admits well-defined equilibrium steady states, RD systems
are said to be detailed balanced. Nonequilibrium steady states
$\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ are characterized by

$j_{\rho}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))\neq 0\text{ }\text{
}\text{ }\text{ and/or }\text{ }\text{ }\text{
}\boldsymbol{J}_{i}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))\neq 0\,,$
(17)

for some $\rho\in\mathcal{R}$ and/or $i\in\mathcal{S}$.

If the steady-state concentration fields are constant for every
$\boldsymbol{r}\in V$, they are said to be homogeneous and represented by
$\boldsymbol{c}_{\mathrm{eq}}^{\mathrm{h}}$ and
$\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}$ for equilibrium and nonequilibrium
steady states, respectively. In general, steady-state concentration fields are
$\boldsymbol{r}$-dependent and are said to be inhomogeneous. Also equilibrium
steady states can be inhomogeneous when chemical species interact (as
predicted by the Cahn-Hilliard Cahn and Hilliard (1958) and Flory-Huggins
Huggins (1941); Flory (1942) theories).

###  III.2 Conservation Laws

The conservation laws Polettini and Esposito (2014); Rao and Esposito (2016);
Falasco  _et al._ (2018) of the RD dynamics (14) (labeled $\lambda\in\Lambda$)
are linearly-independent left-null eigenvectors
$\boldsymbol{\ell}^{\lambda}=(\dots,{\ell}^{\lambda}_{i},\dots)$ of the
stoichiometric matrix, namely,

$\sum_{i}{\ell}^{\lambda}_{i}\,{S}_{i,\rho}=0\quad\forall\rho\,.$ (18)

They physically identify fragments of (or entire) molecules, called moieties,
that remain unaltered by the chemical reactions and the diffusion processes:
they are transferred from one species to another via the chemical reactions
and from one point in space to another via the diffusion processes. Indeed,
the moiety abundances, given by

$L^{\lambda}[\boldsymbol{c}]=\sum_{i}{\ell}^{\lambda}_{i}\int_{V}\mathrm{d}\boldsymbol{r}\,c_{i}(\boldsymbol{r})\,,$
(19)

are conserved by the RD dynamics (14) in closed RD systems (i.e., when
$I_{i}(\boldsymbol{r})=0$ $\forall i$):

$\begin{split}\mathrm{d}_{t}L^{\lambda}[\boldsymbol{c}]=&\sum_{\rho>0}\underbrace{\sum_{i}{\ell}^{\lambda}_{i}\,{S}_{i,\rho}}_{=0}\int_{V}\mathrm{d}\boldsymbol{r}\,j_{\rho}(\boldsymbol{c}(\boldsymbol{r}))\\\
&-\sum_{i}{\ell}^{\lambda}_{i}\underbrace{\int_{V}\mathrm{d}\boldsymbol{r}\,\boldsymbol{\nabla}\cdot\boldsymbol{J}_{i}(\boldsymbol{c}(\boldsymbol{r}))}_{=0}=0\,,\end{split}$
(20)

where we used the definition of conservation law in Eq. (18) and the fact that
the volume $V$ has impermeable walls, i.e.,
$\int_{V}\mathrm{d}\boldsymbol{r}\,\boldsymbol{\nabla}\cdot\boldsymbol{J}_{i}(\boldsymbol{c}(\boldsymbol{r}))=0$.
Note that there is always at least the conservation law
$\boldsymbol{\ell}^{m}=(\dots,m_{i},\dots)$ (with $m_{i}$ being the molecular
mass of species $i$) since the total mass is conserved by the chemical
reactions. Furthermore, the set of conservation laws
$\\{\boldsymbol{\ell}^{\lambda}\\}$ is not unique: a linear combination of
conservation laws is still a conservation law. Different sets (i.e., different
representations) identify different moieties.

In open RD systems, when the $Y$ species are chemostatted, conservation laws
are split into unbroken and broken conservation laws, labeled
$\lambda_{u}\in\Lambda_{u}$ and $\lambda_{b}\in\Lambda_{b}$, respectively Rao
and Esposito (2016); Avanzini  _et al._ (2021). The former are the largest
subset of conservation laws that can be written with null entries for the
chemostatted species, i.e., ${\ell}^{\lambda_{u}}_{y}=0$ $\forall y$ and
$\forall\lambda_{u}$, and hence
$\sum_{x}{\ell}^{\lambda_{u}}_{x}{S}_{x,\rho}=0$. They identify moieties that
are carried by the internal species only and, consequently, are not exchanged
with the chemostats. The corresponding abundances are still conserved,

$\mathrm{d}_{t}L^{\lambda_{u}}[\boldsymbol{c}]=\sum_{y}\underbrace{{\ell}^{\lambda_{u}}_{y}}_{=0}\int_{V}\mathrm{d}\boldsymbol{r}\,I_{y}(\boldsymbol{r})=0\,,$
(21)

despite the RD system being open. The latter are the other conservation laws:
$\Lambda_{b}=\Lambda\setminus\Lambda_{u}$. They identify moieties that are
carried by the chemostatted species too and, consequently, are exchanged with
the chemostats. The corresponding abundances are, in general, not conserved:

$\mathrm{d}_{t}L^{\lambda_{b}}[\boldsymbol{c}]=\sum_{y}\underbrace{{\ell}^{\lambda_{b}}_{y}}_{\neq
0}\int_{V}\mathrm{d}\boldsymbol{r}\,I_{y}(\boldsymbol{r})\neq 0\,.$ (22)

Note that open RD systems have always at least one broken conservation law
since the total mass is not conserved by the exchange processes and,
consequently, $\boldsymbol{\ell}^{m}$ is broken.

In Subs. IV.4, we will use the unbroken conservation laws to quantify the free
energy exchanged between RD systems and chemostats. To this aim, it is
important to recognize that chemostatting a species does not always break a
conservation law Falasco  _et al._ (2018); Avanzini  _et al._ (2021). Hence,
the set of chemostatted species $Y$ (and the set of the corresponding
chemostats) can be split into two disjoint subsets.

The potential species (labeled ${y_{p}{}}\in Y_{p}$) are the smallest subset
of chemostatted species such that all conservation laws
$\\{\boldsymbol{\ell}^{\lambda_{b}}\\}$ with $\lambda_{b}\in\Lambda_{b}$ are
broken. Crucially, there is always a representation of the broken conservation
laws such that each corresponding moiety is exchanged with only one potential
chemostat. Indeed, the matrix (see the graphical representation in Fig. 2)
whose entries are $\\{{\ell}^{\lambda_{b}}_{{y_{p}{}}}\\}$ (with
$\lambda_{b}\in\Lambda_{b}$ and ${y_{p}{}}\in Y_{p}$) is square and invertible
Rao and Esposito (2016); Avanzini  _et al._ (2024) and, by labeling
$\\{\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\\}$ the entries of the inverse
matrix, the broken conservation laws can be represented in terms of the
following linear combination as

$\sum_{\lambda_{b}}\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\,\boldsymbol{\ell}^{\lambda_{b}}\,.$
(23)

The abundance of each moiety in this representation

$M_{{y_{p}{}}}[\boldsymbol{c}]=\sum_{i}\sum_{\lambda_{b}}\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\,{\ell}^{\lambda_{b}}_{i}\int_{V}\mathrm{d}\boldsymbol{r}\,c_{i}(\boldsymbol{r})\,,$
(24)

where $i$ runs over all specie $\mathcal{S}$, changes (according to the RD
equation (14)) only because of the exchanges with one specific potential
chemostat:

$\small\mathrm{d}_{t}M_{{y_{p}{}}}[\boldsymbol{c}]=\sum_{{y_{p}{}^{\prime}}}\underbrace{\sum_{\lambda_{b}}\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\,{\ell}^{\lambda_{b}}_{{y_{p}{}^{\prime}}}}_{={\mathbb{1}}_{{y_{p}{}},{y_{p}{}^{\prime}}}}\int_{V}\mathrm{d}\boldsymbol{r}\,I_{{y_{p}{}^{\prime}}}(\boldsymbol{r})=\int_{V}\mathrm{d}\boldsymbol{r}\,I_{{y_{p}{}}}(\boldsymbol{r})$
(25)

if all chemostatted species are potential species, i.e., $Y=Y_{p}$. This
allows us to label each moiety in Eq. (24) using the label of the chemostat it
is exchanged with, i.e., ${y_{p}{}}$, rather than the label of the
corresponding broken conservation law, i.e., $\lambda_{b}$. Note that Eq. (25)
specializes Eq. (22) for open RD systems with only potential chemostatted
species and with the conservation laws expressed according to the
representation in Eq. (23). We emphasize that the representation in Eq. (23)
will be the one used to quantify the free energy exchanged between RD systems
and chemostats by entering the definition of the thermodynamic potential in
Eq. (41).

Figure 2:  Representation of the matrix with entries
$\\{{\ell}^{\lambda_{b}}_{i}\\}$ where $\lambda_{b}$ is the row index and $i$
is the column index (namely, the matrix constructed using the broken
conservation laws $\\{\boldsymbol{\ell}^{\lambda_{b}}\\}$ as rows). The matrix
column (resp. the entries) can be further specialized into those corresponding
to the internal species $x$ (resp. ${\ell}^{\lambda_{b}}_{x}$), force
chemostatted species ${y_{f}{}}$ (resp. ${\ell}^{\lambda_{b}}_{{y_{f}{}}}$),
and potential chemostatted species ${y_{p}{}}$ (resp.
${\ell}^{\lambda_{b}}_{{y_{p}{}}}$). The submatrix with entries
$\\{{\ell}^{\lambda_{b}}_{{y_{p}{}}}\\}$ (red dotted line) is square and
invertible Rao and Esposito (2016); Avanzini  _et al._ (2024) and we label
$\\{\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\\}$ the entries of the inverse
matrix where ${y_{p}{}}$ is the row index and $\lambda_{b}$ is the column
index.

The force species (${y_{f}{}}\in Y_{f}$) are the other chemostatted species:
$Y_{f}=Y\setminus Y_{p}$. The corresponding chemostats exchange moieties that
are already exchanged with the potential chemostats. Indeed, in general, the
abundance of each moiety in the representation given in Eq. (23) changes
according to

$\mathrm{d}_{t}M_{{y_{p}{}}}[\boldsymbol{c}]=\int_{V}\mathrm{d}\boldsymbol{r}\,I_{{y_{p}{}}}(\boldsymbol{r})+\sum_{{y_{f}{}}}\sum_{\lambda_{b}}\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\,{\ell}^{\lambda_{b}}_{{y_{f}{}}}\int_{V}\mathrm{d}\boldsymbol{r}\,I_{{y_{f}{}}}(\boldsymbol{r})\,,$
(26)

using Eqs. (14) and (24) and the splitting $Y=Y_{p}\bigcup Y_{f}$. Hence,
chemostatting the force species (after the potential species were already
chemostatted) does not break any new conservation law, but creates a flux of
the same moiety between different chemostats. The energetic cost of sustaining
these fluxes will enter the definition of the nonconsevative work in Eq. (42).

###  III.3 Example

A RD system with the chemical species and reactions introduced in Subs. II.2
admits only one conservation law (see the stoichiometric matrix in Fig. 1),
i.e.,

$\boldsymbol{\ell}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr}}}}\right)$}}\,.$ (27)

This physically means that all species
$\\{\ch{X_{1}},\ch{Y_{1}},\ch{Y_{2}},\ch{Y_{3}}\\}$ are made of the same
atoms, namely, the same moiety. Therefore,
$L[\boldsymbol{c}]=\int_{V}\mathrm{d}\boldsymbol{r}\,(c_{\ch{X_{1}}}(\boldsymbol{r})+c_{\ch{Y_{1}}}(\boldsymbol{r})+c_{\ch{Y_{2}}}(\boldsymbol{r})+c_{\ch{Y_{3}}}(\boldsymbol{r}))$
is conserved if no species are chemostatted. By chemostatting just one of the
species $\\{\ch{Y_{1}},\ch{Y_{2}},\ch{Y_{3}}\\}$, the conservation law is
broken. Hence, any of the species $\\{\ch{Y_{1}},\ch{Y_{2}},\ch{Y_{3}}\\}$ can
be chosen as a potential species. We choose Y_1. By chemostatting another
species, like for instance Y_2, no conservation law is broken. Hence, Y_2 is a
force species and the same moiety is exchanged between 2 different chemostats
as illustrated in Fig. 3. Note that this implies that the matrix with entries
$\\{{\ell}^{\lambda_{b}}_{{y_{p}{}}}\\}$ has just one entry corresponding to
${\ell}^{\lambda_{b}}_{\ch{Y_{1}}}=1$ and that the conservation law in Eq.
(27) is already written in the representation given in Eq. (23).

Figure 3: Illustration of the sequences of reactions transferring the same
moiety from a potential chemostat to a force chemostat (for a RD system with
the chemical species and reactions introduced in Subs. II.2). (a) When only
the species Y_1 and Y_2 are chemostatted, the moiety can be injected into the
RD system by a chemostat as the species $\ch{Y_{1}}$, transferred via the
sequence of reaction $\\{1,-2\\}$ to the species $\ch{Y_{2}}$, and finally
extracted from the RD system by a different chemostat as the species
$\ch{Y_{2}}$. (b) When only the species Y_1 and Y_3 are chemostatted, the
moiety can be injected into the RD system by a chemostat as the species
$\ch{Y_{1}}$, transferred via the sequence of reaction $\\{1,-3\\}$ to the
species $\ch{Y_{3}}$, and finally extracted from the RD system by a different
chemostat as the species $\ch{Y_{3}}$. Note that when all three species
$\\{\ch{Y_{1}},\ch{Y_{2}},\ch{Y_{3}}\\}$ are chemostatted, the same moiety can
be transferred between each pair of chemostats.

##  IV Thermodynamics

We now develop a nonequilibrium thermodynamic theory for open RD systems of
interacting chemical species undergoing the RD dynamics (14) by combining the
thermodynamics of noninteracting (i.e., ideal) RD systems Falasco  _et al._
(2018); Avanzini  _et al._ (2019) and interacting (i.e., non-ideal) CRNs
Avanzini  _et al._ (2021).

###  IV.1 Nonequilibrium Free Energy

Nonequilibrium thermodynamics presumes that all degrees of freedom other than
concentration fields are equilibrated (e.g., the temperature $T$ and the
volume $V$ are fixed, while the fields responsible for the interactions, e.g.,
electrostatic fields, relax instantaneously to mean-field values depending on
the concentration fields only). In this framework, the nonequilibrium
Helmholtz free energy $F[\boldsymbol{c}]$ of RD systems can be specified by
its equilibrium form, but as a functional of the nonequilibrium concentration
fields. Furthermore, $F[\boldsymbol{c}]$ can always be written as the sum of
an ideal and a non-ideal contribution:

$F[\boldsymbol{c}]=F^{\mathrm{id}}[\boldsymbol{c}]+F^{\mathrm{ni}}[\boldsymbol{c}]\,.$
(28)

Here, the ideal term $F^{\mathrm{id}}[\boldsymbol{c}]$ is the ideal solution
Helmholtz free energy, while the non-ideal term
$F^{\mathrm{ni}}[\boldsymbol{c}]$ accounts for the interactions between the
chemical species.

The free energy contribution of each chemical species is specified by the
chemical potential

$\mu_{i}(\boldsymbol{c}(\boldsymbol{r}))=\frac{\delta
F[\boldsymbol{c}]}{\delta
c_{i}(\boldsymbol{r})}=\mu^{\mathrm{id}}_{i}(c_{i}(\boldsymbol{r}))+RT\ln\gamma_{i}(\boldsymbol{c}(\boldsymbol{r}))\,,$
(29)

where $\mu^{\mathrm{id}}_{i}(c_{i}(\boldsymbol{r}))$ is the ideal chemical
potential

$\mu^{\mathrm{id}}_{i}(c_{i}(\boldsymbol{r}))=\mu_{i}^{\circ}+RT\ln{c_{i}(\boldsymbol{r})}\,,$
(30)

with $\mu_{i}^{\circ}$ the standard chemical potential, $R$ the gas constant,
and $\gamma_{i}(\boldsymbol{c}(\boldsymbol{r}))$ the activity coefficient
accounting for the interactions,

$RT\ln\gamma_{i}(\boldsymbol{c}(\boldsymbol{r}))=\frac{\delta
F^{\mathrm{ni}}[\boldsymbol{c}]}{\delta c_{i}(\boldsymbol{r})}\,.$ (31)

Note that the ideal free energy $F^{\mathrm{id}}[\boldsymbol{c}]$ reads
Falasco  _et al._ (2018); Avanzini  _et al._ (2019)

$F^{\mathrm{id}}[\boldsymbol{c}]=\sum_{i}\int_{V}\mathrm{d}\boldsymbol{r}\,{\big{(}\mu^{\mathrm{id}}_{i}(c_{i}(\boldsymbol{r}))-RT\big{)}c_{i}(\boldsymbol{r})}\,.$
(32)

###  IV.2 Thermodynamic Consistency of the RD Dynamics

We start by introducing the entropy production rate

$\dot{\Sigma}[\boldsymbol{c}]\equiv\dot{\Sigma}_{\mathrm{rct}}[\boldsymbol{c}]+\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}]\,,$
(33)

that quantifies the dissipation resulting from the changes in chemical
potentials caused by the two types of processes occurring in the system,
namely, reactions and diffusion:

$\displaystyle T\dot{\Sigma}_{\mathrm{rct}}[\boldsymbol{c}]$
$\displaystyle\equiv-\sum_{i}\sum_{\rho>0}\int_{V}\mathrm{d}\boldsymbol{r}\,\mu_{i}(\boldsymbol{c}(\boldsymbol{r})){S}_{i,\rho}j_{\rho}(\boldsymbol{c}(\boldsymbol{r}))\,,$
(34a) $\displaystyle T\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}]$
$\displaystyle\equiv-\sum_{i}\int_{V}\mathrm{d}\boldsymbol{r}\,\boldsymbol{\nabla}\mu_{i}(\boldsymbol{c}(\boldsymbol{r}))\cdot\boldsymbol{J}_{i}(\boldsymbol{c}(\boldsymbol{r}))\,.$
(34b)

The central assumption ensuring thermodynamic consistency of the RD dynamics
is that the reaction fluxes $\omega_{\pm\rho}(\boldsymbol{c}(\boldsymbol{r}))$
satisfy the local detailed balance condition:

$\small
RT\ln\frac{\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))}{\omega_{-\rho}(\boldsymbol{c}(\boldsymbol{r}))}=-\Big{(}\sum_{x}\mu_{x}(\boldsymbol{c}(\boldsymbol{r})){S}_{x,\rho}+\sum_{y}\mu_{y}(\boldsymbol{c}(\boldsymbol{r})){S}_{y,\rho}\Big{)}\,.$
(35)

We further assume that the diffusion current
$\boldsymbol{J}_{i}(\boldsymbol{c}(\boldsymbol{r}))$ of each species is
treated within linear regime and expressed as

$\boldsymbol{J}_{i}(\boldsymbol{c}(\boldsymbol{r}))=-\sum_{j}O_{i,j}(\boldsymbol{c}(\boldsymbol{r}))\boldsymbol{\nabla}\mu_{j}(\boldsymbol{c}(\boldsymbol{r}))\,,$
(36)

where $\\{O_{i,j}(\boldsymbol{c}(\boldsymbol{r}))\\}$ are the entries of the
positive-definite Onsager matrix $\mathbb{O}(\boldsymbol{c}(\boldsymbol{r}))$
de Groot and Mazur (1984). Using these two assumptions in Eqs. (35) and (36),
we see that the two entropy production rates in Eq. (34) can be rewritten in a
manifestly non-negative form

$\displaystyle T\dot{\Sigma}_{\mathrm{rct}}[\boldsymbol{c}]$
$\displaystyle=\int_{V}\mathrm{d}\boldsymbol{r}\underbrace{RT\sum_{\rho>0}j_{\rho}(\boldsymbol{c}(\boldsymbol{r}))\ln\frac{\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))}{\omega_{-\rho}(\boldsymbol{c}(\boldsymbol{r}))}}_{\equiv
T\dot{\sigma}_{\mathrm{rct}}(\boldsymbol{r})\geq 0}\geq 0\,,$ (37a)
$\displaystyle T\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}]$
$\displaystyle=\int_{V}\mathrm{d}\boldsymbol{r}\underbrace{\sum_{i,j}\boldsymbol{\nabla}\mu_{i}(\boldsymbol{c}(\boldsymbol{r}))\cdot
O_{i,j}(\boldsymbol{c}(\boldsymbol{r}))\boldsymbol{\nabla}\mu_{j}(\boldsymbol{c}(\boldsymbol{r}))}_{\equiv
T\dot{\sigma}_{\mathrm{dff}}(\boldsymbol{r})\geq 0}\geq 0\,,$ (37b)

where the equality sign is achieved only at equilibrium (16).

The change in free energy (28) in closed RD systems, using Eq. (14) (when
$I_{i}(\boldsymbol{r})=0$ $\forall i$) and integrating by parts, gives the
second law of thermodynamics for closed RD systems

$\mathrm{d}_{t}F[\boldsymbol{c}]=-T\dot{\Sigma}[\boldsymbol{c}]\leq 0\,.$ (38)

With the additional assumption that the free energy in Eq. (28) is lower
bounded, $F[\boldsymbol{c}]$ becomes a thermodynamic potential acting as a
Lyapunov function of the RD dynamics (14). This ensures thermodynamic
consistency, namely, closed RD systems always relax towards an equilibrium
steady state.

Remark. In Eqs. (37a) and (37b) we introduced the densities of reaction
$\dot{\sigma}_{\mathrm{rct}}(\boldsymbol{r})$ and diffusion
$\dot{\sigma}_{\mathrm{dff}}(\boldsymbol{r})$ entropy production rate,
respectively. They are non-negative for every $\boldsymbol{r}\in V$ and vanish
only at equilibrium (16). The density of reaction entropy production rate can
be also split into the non-negative contribution of each reaction
$\dot{\sigma}_{\rho}(\boldsymbol{r})\geq 0$:
$\dot{\sigma}_{\mathrm{rct}}(\boldsymbol{r})=\sum_{\rho>0}\dot{\sigma}_{\rho}(\boldsymbol{r})$.

Remark. The expressions of entropy production rates in Eqs. (37a) and (37b)
together with the local detailed balance condition (35) and
$\mathbb{O}(\boldsymbol{c}(\boldsymbol{r}))$ being a positive-definite matrix
imply that the thermodynamic equilibrium conditions,
$\dot{\Sigma}_{\mathrm{rct}}[\boldsymbol{c}_{\mathrm{eq}}]=0$ and
$\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}_{\mathrm{eq}}]=0$, are equivalent
to

$\displaystyle\sum_{i}\mu_{i}(\boldsymbol{c}_{\mathrm{eq}}(\boldsymbol{r})){S}_{i,\rho}=0\,,$
(39a)
$\displaystyle\boldsymbol{\nabla}\mu_{i}(\boldsymbol{c}_{\mathrm{eq}}(\boldsymbol{r}))=0\,,$
(39b)

for all $\boldsymbol{r}\in V$, for all $\rho\in\mathcal{R}$, and for all
$i\in\mathcal{S}$. The equilibrium conditions in Eq. (39) can also be derived
by minimizing the free energy (28) as done in App. A.

###  IV.3 Chemostats

Thermodynamically, chemostats control the chemical potentials
$\mu_{y}(\boldsymbol{c}(\boldsymbol{r}))$ of the chemostatted species in every
point of space $\boldsymbol{r}\in V$, which can thus be treated as (in
general, $\boldsymbol{r}$\- and time-dependent) control parameters, i.e,
$\\{\mu_{y}(\boldsymbol{r})\\}$. In ideal RD systems, where the chemical
potentials are given in Eq. (30), chemostats control the concentration fields
of the chemostatted species too:
$c_{y}(\boldsymbol{r})=\exp((\mu_{y}(\boldsymbol{r})-\mu_{y}^{\circ})/RT)$.
The exchange current $I_{y}(\boldsymbol{r})$ in Eq. (14) represents the
physical mechanism implementing this control. In non-ideal RD systems, where
the chemical potentials are given in Eq. (29), controlling the chemical
potentials of the chemostatted species do not univocally determine their
concentrations because of the interactions Avanzini  _et al._ (2021).

###  IV.4 Second Law for Open RD Systems

Using the expression for the entropy production rate in Eq. (33), the second
law of thermodynamics for open RD systems can now be rewritten as

$T\dot{\Sigma}[\boldsymbol{c}]=-\mathrm{d}_{t}\mathcal{F}[\boldsymbol{c}]+\dot{W}_{\mathrm{nc}}[\boldsymbol{c}]+\dot{W}_{\mathrm{driv}}[\boldsymbol{c}]\geq
0\,.$ (40)

This can be checked by direct substitution of the different terms on the
right-hand side that we now define together with the RD equation (14).

The semigrand Helmholtz free energy $\mathcal{F}[\boldsymbol{c}]$ is the
proper thermodynamic potential of open RD systems. It reads

$\mathcal{F}[\boldsymbol{c}]=F[\boldsymbol{c}]-\sum_{{y_{p}{}}}\mu_{{y_{p}{}}}^{\mathrm{ref}}\,M_{{y_{p}{}}}[\boldsymbol{c}]$
(41)

and is obtained, in analogy to equilibrium thermodynamics when passing from
the canonical to the grand canonical ensembles, from the Helmholtz free energy
(28) by removing the energetic contribution of the matter exchanged with the
chemostats. The latter accounts for the abundances of the moieties
$\\{M_{{y_{p}{}}}[\boldsymbol{c}]\\}$, given in Eq. (24), times a reference
value of chemical potentials of the potential species
$\\{\mu_{{y_{p}{}}}^{\mathrm{ref}}\\}$, controlled by the corresponding
chemostats. If the potential chemostats impose different values of the
chemical potentials $\\{\mu_{{y_{p}{}}}(\boldsymbol{r})\\}$ in different
$\boldsymbol{r}\in V$, then the reference chemical potential
$\\{\mu_{{y_{p}{}}}^{\mathrm{ref}}\\}$ can be chosen arbitrarily among these
values. This is equivalent to choosing a reference equilibrium steady state
$\boldsymbol{c}_{\mathrm{eq}}(\boldsymbol{r})$ satisfying
$\mathcal{F}[\boldsymbol{c}_{\mathrm{eq}}]\leq\mathcal{F}[\boldsymbol{c}]$ as
discussed in App. B.

The nonconservative work rate

$\small\dot{W}_{\mathrm{nc}}[\boldsymbol{c}]=\int_{V}\mathrm{d}\boldsymbol{r}\,\underbrace{\sum_{y}f^{\mathrm{nc}}_{y}(\boldsymbol{r})I_{y}(\boldsymbol{r})}_{\equiv\dot{w}_{\mathrm{nc}}(\boldsymbol{r})}\\\
$ (42)

quantifies the energetic cost for sustaining fluxes of the same moiety between
chemostats with different chemical potentials by means of the nonconservative
forces

$f^{\mathrm{nc}}_{y}(\boldsymbol{r})\equiv\mu_{y}(\boldsymbol{r})-\sum_{{y_{p}{}},\lambda_{b}}\mu_{{y_{p}{}}}^{\mathrm{ref}}\,\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\,{\ell}^{\lambda_{b}}_{y}\,.$
(43)

These forces emerge in two cases. When
$\mu_{{y_{p}{}}}(\boldsymbol{r})\neq\mu_{{y_{p}{}}}^{\mathrm{ref}}$ (because
$\sum_{\lambda_{b}}\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\,{\ell}^{\lambda_{b}}_{{y_{p}{}^{\prime}}}={\mathbb{1}}_{{y_{p}{}},{y_{p}{}^{\prime}}}$),
the forces $f^{\mathrm{nc}}_{{y_{p}{}}}(\boldsymbol{r})$ correspond to
diffusive forces that can be present even in the absence of chemical reactions
Falasco  _et al._ (2018). They arise because the potential chemostats
${y_{p}{}}$ impose values of the chemical potential that are different from
the reference equilibrium ones $\mu_{{y_{p}{}}}^{\mathrm{ref}}$ (see App. B).
When
$\mu_{{y_{f}{}}}(\boldsymbol{r})\neq\sum_{{y_{p}{}},\lambda_{b}}\mu_{{y_{p}{}}}^{\mathrm{ref}}\,\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\,{\ell}^{\lambda_{b}}_{{y_{f}{}}}$,
the forces $f^{\mathrm{nc}}_{{y_{f}{}}}(\boldsymbol{r})$ correspond to
chemical forces that can also be present in homogeneous systems Avanzini  _et
al._ (2021). They arise because force chemostats ${y_{f}{}}$ impose chemical
potentials that are different from the one they would have if all reactions
were equilibrated
$\sum_{{y_{p}{}},\lambda_{b}}\mu_{{y_{p}{}}}^{\mathrm{ref}}\,\overline{\ell}{}_{\lambda_{b}}^{{y_{p}{}}}\,{\ell}^{\lambda_{b}}_{{y_{f}{}}}$
(see App. B).

The driving work rate reads

$\dot{W}_{\mathrm{driv}}[\boldsymbol{c}]=-\sum_{{y_{p}{}}}\big{(}\mathrm{d}_{t}\mu_{{y_{p}{}}}^{\mathrm{ref}}\big{)}\,M_{{y_{p}{}}}[\boldsymbol{c}]$
(44)

and quantifies the energetic cost of the time-dependent changes of the
reference chemical potentials $\\{\mu_{{y_{p}{}}}^{\mathrm{ref}}\\}$ (with
${y_{p}{}}\in Y_{p}$) defining the reference equilibrium to which the open RD
system would relax if the nonconservative forces vanished.

Physically, the second law (40) expresses (and quantifies) the different
sources of free energy (on the right-hand side of Eq. (40)) that can balance
dissipation (on the left-hand side of Eq. (40)) and hence maintain open RD
systems out of equilibrium. Let us illustrate this point by considering three
simple cases.

First, when the chemical potentials of the chemostatted species are
independent of time and such that $f^{\mathrm{nc}}_{y}(\boldsymbol{r})=0$
(which also implies homogeneous chemostatting), both work contributions vanish
and the second law (40) reduces to
$\mathrm{d}_{t}\mathcal{F}[\boldsymbol{c}]=-T\dot{\Sigma}[\boldsymbol{c}]\leq
0$. This, together with $\mathcal{F}[\boldsymbol{c}]$ being lower bounded (as
shown in App. B), implies that the RD system relaxes to equilibrium while
minimizing $\mathcal{F}[\boldsymbol{c}]$ which serves as a Lyapunov function.
This RD system is therefore detailed balanced despite being open. As shown in
App. B, the specific equilibrium steady state
$\boldsymbol{c}_{\mathrm{eq}}(\boldsymbol{r})$ to which the RD system relaxes
is defined by the chemical potentials $\\{\mu_{{y_{p}{}}}^{\mathrm{ref}}\\}$
and the abundances $\\{L^{\lambda_{u}}[\boldsymbol{c}]\\}$.

Second, when the chemical potentials of the chemostatted species are dependent
on time but $f^{\mathrm{nc}}_{y}(\boldsymbol{r})=0$ still holds, only the
nonconservative work rate vanishes from the second law (40), namely,
$T\dot{\Sigma}[\boldsymbol{c}]=-\mathrm{d}_{t}\mathcal{F}[\boldsymbol{c}]+\dot{W}_{\mathrm{driv}}[\boldsymbol{c}]\geq
0$, and the RD system is still detailed balanced. Indeed, if the time
dependence is very slow, the RD system will quasi-statically follow the
changing equilibrium state. For faster time dependence, the RD system will be
prevented from reaching equilibrium.

Third, when the chemical potentials of the chemostatted species are
independent of time but $f^{\mathrm{nc}}_{y}(\boldsymbol{r})\neq 0$, detailed
balance is broken and the second law (40) reduces to
$T\dot{\Sigma}[\boldsymbol{c}]=-\mathrm{d}_{t}\mathcal{F}[\boldsymbol{c}]+\dot{W}_{\mathrm{nc}}[\boldsymbol{c}]\geq
0$. If the RD system eventually reaches a nonequilibrium steady state, the
second law is further reduced to
$T\dot{\Sigma}[\boldsymbol{c}_{\mathrm{ss}}]=\dot{W}_{\mathrm{nc}}[\boldsymbol{c}_{\mathrm{ss}}]\geq
0$. In this case, $\dot{W}_{\mathrm{nc}}[\boldsymbol{c}_{\mathrm{ss}}]$
quantifies the net free energy intake by the open RD system that is dissipated
to maintain the system in a nonequilibrium steady state.

###  IV.5 Example

We consider the RD system with the chemical species and reactions introduced
in Subs. II.2 and with the chemostats imposing $\boldsymbol{r}$-dependent
chemical potentials $\mu_{\ch{Y_{1}}}(\boldsymbol{r})$,
$\mu_{\ch{Y_{2}}}(\boldsymbol{r})$, and $\mu_{\ch{Y_{3}}}(\boldsymbol{r})$. We
choose $\ch{Y_{1}}$ as potential species (as discussed in Subs. III.3) and one
of the values of $\mu_{\ch{Y_{1}}}(\boldsymbol{r})$ as reference chemical
potential $\mu^{\mathrm{ref}}_{\ch{Y_{1}}}$. The corresponding nonconservative
work rate (whose general expression is give in Eq. (42)) specializes into the
sum of three contributions:

$\begin{split}\dot{W}_{\mathrm{nc}}[\boldsymbol{c}]=&\int_{V}\mathrm{d}\boldsymbol{r}\,(\mu_{\ch{Y_{1}}}(\boldsymbol{r})-\mu^{\mathrm{ref}}_{\ch{Y_{1}}})I_{\ch{Y_{1}}}(\boldsymbol{r})\\\
+&\int_{V}\mathrm{d}\boldsymbol{r}\,(\mu_{\ch{Y_{2}}}(\boldsymbol{r})-\mu^{\mathrm{ref}}_{\ch{Y_{1}}})I_{\ch{Y_{2}}}(\boldsymbol{r})\\\
+&\int_{V}\mathrm{d}\boldsymbol{r}\,(\mu_{\ch{Y_{3}}}(\boldsymbol{r})-\mu^{\mathrm{ref}}_{\ch{Y_{1}}})I_{\ch{Y_{3}}}(\boldsymbol{r})\,,\end{split}$
(45)

where we used the broken conservation law in Eq. (27). The first contribution
results from applying different chemical potentials to the species
$\ch{Y_{1}}$. The second (resp. third) contribution results from applying
different chemical potentials to the same moiety carried by both the species
$\ch{Y_{1}}$ and $\ch{Y_{2}}$ (resp. $\ch{Y_{3}}$) as illustrated in Fig. 3.

##  V Thermodynamic Properties of Active Self-Organization - Exact Results

The RD dynamics might reach a steady state where the chemical species are
inhomogeneous. We now use the thermodynamic theory introduced in Sec. IV to
prove that for two classes of CRNs, namely, pseudo detailed balanced (Subs.
V.1) and complex balanced (Subs. V.2) discussed in Fig. 4, stationary
concentration profiles are determined by the minimization of appropriately
constructed kinetic potentials acting as Lyapunov functions, when the chemical
potentials $\\{\mu_{y}\\}$ of the chemostatted species are homogeneous and
constant in time. Furthermore, for both classes of CRNs, diffusion processes
always equilibrate, which implies that the chemical potentials become
homogeneous at steady state. Additionally, for pseudo-detailed CRNs only, the
transient dynamics and steady-state concentration fields can be exactly
simulated also in reaction-diffusion systems relaxing towards equilibrium with
appropriately constructed chemical potentials.

The proofs of this result are derived in subs. V.1 and V.2 by assuming that
the chemostatted species are ideal, i.e., they do not interact with the
internal species, for the sake of simplicity. We can thus treat both
$\\{\mu_{y}\\}$ and $\\{c_{y}\\}$ (related via Eq. (30)) as control
parameters. We generalize our result to non-ideal chemostatted species in sub.
V.3.

Figure 4:  The gray box represents the set of all possible CRNs endowed with
all possible reaction fluxes. The blue box, a subset of the gray box,
represents the set of all possible complex balanced CRNs endowed with all
possible reaction fluxes. Pseudo detailed balanced CRNs are always a subset of
complex balanced CRNs. In Sec. V, we consider pseudo detailed balanced CRNs
with reaction fluxes of general form (green oval) and complex balanced CRNs
with Arrhenius-like reaction fluxes (red oval).

###  V.1 Pseudo Detailed Balanced CRNs

We define CRNs as pseudo detailed balanced when their stoichiometric matrix
for the internal species and lumped reactions does not admit any right-null
eigenvector. Namely, there is no $\\{\phi_{\varepsilon}\\}$ such that

$\sum_{\varepsilon>0}{S}_{x,\varepsilon}\,\phi_{\varepsilon}=0\,.$ (46)

We now show that the RD dynamics of pseudo detailed balanced CRNs is driven by
the minimization of a kinetic potential and relaxes to a steady state
$\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ characterized by a vanishing
diffusion entropy production rate, i.e.,
$\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}_{\mathrm{ss}}]=0$, independently
of whether the concentration fields are homogeneous or inhomogeneous.

Assumption. We consider reaction fluxes of the general form

$\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))=k_{\rho}\omega_{\mathrm{in}}(\boldsymbol{c}(\boldsymbol{r}),\\{\nu_{x,\rho}\\})\omega_{\mathrm{ch}}(\\{c_{y}\\},\\{\nu_{y,\rho}\\})\,,$
(47)

namely, the concentration-dependence of each flux
$\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))$ is univocally determined by
the stoichiometric coefficients via the unspecified functions
$\omega_{\mathrm{in}}(\boldsymbol{c}(\boldsymbol{r}),\\{\nu_{x,\rho}\\})$ and
$\omega_{\mathrm{ch}}(\\{c_{y}\\},\\{\nu_{y,\rho}\\})$. The latter depends
only on the concentrations of chemostatted species because they are assumed to
be ideal. Here, $k_{\rho}$ is an unspecified constant parameter. For example,
mass-action fluxes Laidler (1987); Pekař (2005) and Arrhenius-like fluxes
satisfy Eq. (47). Note that, while the latter read

$\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))=A_{\rho}\,e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}(\boldsymbol{r})){\nu}_{x,\rho}}{RT}}e^{\frac{\sum_{y}\mu_{y}{\nu}_{y,\rho}}{RT}}\,,$
(48)

with $A_{\rho}=A_{-\rho}$ being a constant parameter, the fluxes in Eq. (47)
can in general account for a concentration-dependent $A_{\rho}$.

Proof. We start by using ${S}_{x,\rho}={S}_{x,\varepsilon}$
$\forall\rho\in\mathcal{R}_{\varepsilon}$ to rewrite the contribution due to
the chemical reactions in Eq. (14) for the internal species as

$\sum_{\rho>0}{S}_{x,\rho}\,j_{\rho}(\boldsymbol{c}(\boldsymbol{r}))=\sum_{\varepsilon>0}{S}_{x,\varepsilon}\hat{j}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))\,,$
(49)

with

$\hat{j}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))=\sum_{\rho\in\mathcal{R}_{\varepsilon}}\big{\\{}\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))-\omega_{-\rho}(\boldsymbol{c}(\boldsymbol{r}))\big{\\}}$
(50)

being the lumped reaction $\varepsilon$ net current. By writing now each
reaction flux $\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))$ as the product
between its symmetric part
$s_{\rho}(\boldsymbol{c}(\boldsymbol{r}))\equiv\sqrt{\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))\omega_{-\rho}(\boldsymbol{c}(\boldsymbol{r}))}=s_{-\rho}(\boldsymbol{c}(\boldsymbol{r}))$
and its antisymmetric part
$\sqrt{\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))/\omega_{-\rho}(\boldsymbol{c}(\boldsymbol{r}))}$
specified by the local detailed balance condition (35),
$\hat{j}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))$ becomes

$\small\hat{j}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))=a_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))\Big{\\{}e^{-\frac{\sum_{x}\mu_{x}(\boldsymbol{c}(\boldsymbol{r})){S}_{x,\varepsilon}}{2RT}}-b_{\varepsilon}e^{-\frac{\sum_{x}\mu_{x}(\boldsymbol{c}(\boldsymbol{r})){S}_{x,-\varepsilon}}{2RT}}\Big{\\}}\,,$
(51)

where we factorized

$a_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))=\sum_{\rho\in\mathcal{R}_{\varepsilon}}s_{\rho}(\boldsymbol{c}(\boldsymbol{r}))e^{-\frac{\sum_{y}\mu_{y}{S}_{y,\rho}}{2RT}}\geq
0\,,$ (52)

and the constant quantity

$b_{\varepsilon}=\frac{a_{-\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))}{a_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))}=\frac{\sum_{\rho\in\mathcal{R}_{\varepsilon}}s_{\rho}(\boldsymbol{c}(\boldsymbol{r}))e^{-\frac{\sum_{y}\mu_{y}{S}_{y,-\rho}}{2RT}}}{\sum_{\rho\in\mathcal{R}_{\varepsilon}}s_{\rho}(\boldsymbol{c}(\boldsymbol{r}))e^{-\frac{\sum_{y}\mu_{y}{S}_{y,\rho}}{2RT}}}\geq
0.$ (53)

This results from the fact that i) the chemical potentials $\\{\mu_{y}\\}$,
and the corresponding concentrations $\\{c_{y}\\}$, are homogeneous and
constant, and ii) and the fluxes in $s_{\rho}(\boldsymbol{c}(\boldsymbol{r}))$
are univocally determined by the stoichiometric coefficients according to Eq.
(47) with ${\nu}_{x,\rho}={\nu}_{x,\varepsilon}$
$\forall\rho\in\mathcal{R}_{\varepsilon}$.

We now use the pseudo detailed balanced condition. Since there are no right-
null vectors of ${S}_{x,\varepsilon}$, we can always introduce the constant
coefficients $\\{\Delta_{x}\\}$ such that

$\ln b_{\varepsilon}=\sum_{x}\Delta_{x}\,{S}_{x,\varepsilon}/RT\,.$ (54)

Indeed, given the $n_{x}\times n_{e}$ matrix $\hat{\mathbb{S}}_{X}$ whose
entries are $\\{{S}_{x,\varepsilon}\\}$ with $x\in X$,
$\mathcal{E}\ni\varepsilon>0$, $\mathrm{rank}(\hat{\mathbb{S}}_{X})=n_{e}$.
Therefore, the set of internal species can always be split into two disjoint
sets $X_{I}$ and $X_{D}$ such that the matrix $\hat{\mathbb{S}}_{X_{I}}$ (with
entries $\\{{S}_{x,\varepsilon}\\}$ with $x\in X_{I}$ and
$\mathcal{E}\ni\varepsilon>0$) is square and invertible. This leads to

$\displaystyle\Delta_{x}$ $\displaystyle=RT\sum_{\varepsilon>0}\ln
b_{\varepsilon}(\hat{\mathbb{S}}_{X_{I}}^{-1})_{\varepsilon,x}\text{ }\text{
}\text{ }\text{ for }\text{ }\text{ }\text{ }x\in X_{I}\,,$ (55a)
$\displaystyle\Delta_{x}$ $\displaystyle=0\text{ }\text{ }\text{ }\text{ for
}\text{ }\text{ }\text{ }x\in X_{D}\,,$ (55b)

and allows us to write the lumped reaction $\varepsilon$ net current (51) as

$\small\hat{j}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))=\hat{\omega}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))-\hat{\omega}_{-\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))\,.$
(56)

Here, we introduced the lumped reaction $\pm\varepsilon$ fluxes that read

$\hat{\omega}_{\pm\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))\equiv
a_{|\varepsilon|}(\boldsymbol{c}(\boldsymbol{r}))e^{\frac{\sum_{x}\Delta_{x}{S}_{x,|\varepsilon|}}{2RT}}e^{-\frac{\sum_{x}\hat{\mu}_{x}(\boldsymbol{c}(\boldsymbol{r})){S}_{x,\pm\varepsilon}}{2RT}}\,,$
(57)

with $|-\varepsilon|=|\varepsilon|={\varepsilon}>0$, and satisfy a pseudo
local detailed balance condition

$RT\ln\frac{\hat{\omega}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))}{\hat{\omega}_{-\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))}=-\sum_{x}\hat{\mu}_{x}(\boldsymbol{c}(\boldsymbol{r})){S}_{x,\varepsilon}\,,$
(58)

with the chemical potentials

$\hat{\mu}_{x}(\boldsymbol{c}(\boldsymbol{r}))\equiv\mu_{x}(\boldsymbol{c}(\boldsymbol{r}))+\Delta_{x}\,,$
(59)

Note that Eq. (59) only implies a shift of the standard chemical potentials of
the internal species $\\{\mu_{x}^{\circ}\\}$.

Until now, we just re-wrote the lumped reaction net current
$\hat{j}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))$ entering the RD
equation (14) via Eq. (49) as in Eq. (56). However, this rewriting allows us
to show that the potential ${F}_{\mathrm{pdb}}[\boldsymbol{c}]$, obtained from
the free energy (28) of the RD system after applying the shift of the chemical
potentials in Eq. (59), decreases in time during RD dynamics (14). Indeed,
using $\mathrm{d}_{t}c_{y}=0$, its time derivative reads

$\mathrm{d}_{t}{F}_{\mathrm{pdb}}[\boldsymbol{c}]=-T\hat{\dot{\Pi}}_{\mathrm{rct}}[\boldsymbol{c}]-T\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}]\leq
0\,,$ (60)

where $\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}]\geq 0$ is given in Eq.
(37b) (since $\boldsymbol{\nabla}\mu_{y}=0$) and, using Eq. (58),

$\small
T\hat{\dot{\Pi}}_{\mathrm{rct}}[\boldsymbol{c}]=RT\sum_{\varepsilon>0}\int_{V}\mathrm{d}\boldsymbol{r}\,\hat{j}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))\ln\frac{\hat{\omega}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))}{\hat{\omega}_{-\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))}\geq
0\,.$ (61)

This, together with the fact that ${F}_{\mathrm{pdb}}[\boldsymbol{c}]$ is
lower bounded (since $F[\boldsymbol{c}]$ is lower bounded by assumption, see
Subs. IV.2), implies that the RD dynamics (14) of pseudo detailed balance CRNs
must relax towards a steady state
$\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ satisfying

$\hat{\dot{\Pi}}_{\mathrm{rct}}[\boldsymbol{c}_{\mathrm{ss}}]=0\text{ }\text{
}\text{ }\text{ and }\text{ }\text{ }\text{
}\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}_{\mathrm{ss}}]=0\,,$ (62)

or equivalently

$\sum_{x}(\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))+\Delta_{x}){S}_{x,\varepsilon}=0\text{
}\text{ }\text{ }\text{ and }\text{ }\text{ }\text{
}\boldsymbol{J}_{i}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))=0\,,$ (63)

for all $\varepsilon>0$ and for all $i$. $\square$

Let us now discuss the meaning and implications of Eqs. (62) and (63).

First, the dynamics of pseudo detailed balanced CRNs minimizes
${F}_{\mathrm{pdb}}[\boldsymbol{c}]$ which plays the role of a Lyapunov
function. However, ${F}_{\mathrm{pdb}}[\boldsymbol{c}]$ is not a thermodynamic
potential, but rather a kinetic potential: the chemical potentials (59) encode
kinetic properties, namely, the symmetric part of the reaction fluxes
$\\{s_{\rho}(\boldsymbol{c}(\boldsymbol{r}))\\}$ via $\\{\Delta_{x}\\}$ (see
Eqs. (54) and (53)).

Second, Eqs. (62) and (63) physically mean that diffusion processes are
equilibrated at steady state. This, together with Eq. (39b) and the fact that
$\hat{\mu}_{x}(\boldsymbol{c}(\boldsymbol{r}))$ and
$\mu_{x}(\boldsymbol{c}(\boldsymbol{r}))$ are related through the constant
coefficient $\Delta_{x}$, see Eq. (59), implies

$\boldsymbol{\nabla}\mu_{i}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))=0\quad\forall
i\,,$ (64)

namely, the chemical potentials of all species are homogeneous.

Third, spatial self-organization at steady state leads, in general, to
inhomogeneous reaction fluxes. Indeed, as long as reaction fluxes do not
exclusively depend on the chemical potentials, homogeneous chemical potentials
do not imply homogeneous reaction fluxes. However, this is not the case for,
for example, Arrhenius-like reaction fluxes (48) which become homogeneous when
the chemical potentials are homogeneous (see also Subs. VI.2).

Fourth, by plugging
$\boldsymbol{J}_{i}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))=0$ into Eq.
(14), we obtain that the concentration fields
$\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ are also a steady state of the
chemical dynamics for the internal species in every point $\boldsymbol{r}\in
V$, i.e.,

$\sum_{\rho>0}{S}_{x,\rho}\,j_{\rho}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))=0\,.$
(65)

Since both diffusion and the chemical reactions are at steady state, no
oscillations or moving structures (e.g., chemical waves) can emerge.

Fifth, $\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ is not, in general, an
equilibrium steady state: diffusion processes are equilibrated, i.e.,
$\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}_{\mathrm{ss}}]=0$, but chemical
reactions are not. Indeed,
$\hat{\dot{\Pi}}_{\mathrm{rct}}[\boldsymbol{c}_{\mathrm{ss}}]$ is not the
thermodynamic reaction entropy production rate (34a) and
$\hat{\dot{\Pi}}_{\mathrm{rct}}[\boldsymbol{c}_{\mathrm{ss}}]=0$ in Eq. (62)
does not imply $\dot{\Sigma}_{\mathrm{rct}}[\boldsymbol{c}_{\mathrm{ss}}]=0$.
The expressions of
$\hat{\dot{\Pi}}_{\mathrm{rct}}[\boldsymbol{c}_{\mathrm{ss}}]$ in Eq. (61) and
$\dot{\Sigma}_{\mathrm{rct}}[\boldsymbol{c}_{\mathrm{ss}}]$ in Eq. (37a) have
a similar structure, given in terms of the lumped reaction fluxes
$\hat{\omega}_{\varepsilon}$ and the actual reaction fluxes $\omega_{\rho}$,
respectively. Yet,
$\hat{\dot{\Pi}}_{\mathrm{rct}}[\boldsymbol{c}_{\mathrm{ss}}]$ does not
resolve the dissipation of the single reactions
$\rho\in\mathcal{R}_{\varepsilon}$ in the lumped reaction $\varepsilon$, and
it is always a lower bound of the
$\dot{\Sigma}_{\mathrm{rct}}[\boldsymbol{c}_{\mathrm{ss}}]$ (because of the
log sum inequality).

Sixth, by re-writing the lumped reaction net current
$\hat{j}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))$ according to Eq. (56)
with fluxes satisfying the pseudo local detailed balance condition (58), the
RD dynamics (14) for the internal species is mathematically equivalent to the
RD dynamics of a detailed balanced (and hence driven by purely passive
mechanisms) RD system. For this reason, we named this kind of CRNs pseudo
detailed balanced. Furthermore, Eq. (63) implies that
$(\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))-\Delta_{x})$ can be
written as a linear combination of unbroken conservation laws similar to what
happens at equilibrium where the equilibrium chemical potentials can be
written as a linear combination of conservation laws (see Eq. (39a) and App.
A).

We conclude this discussion by showing that the lumped reaction net currents
$\hat{j}_{\varepsilon}(\boldsymbol{c}(\boldsymbol{r}))$ cannot, in general, be
written as in Eq. (56) when CRNs are not pseudo detailed balanced. Namely,
there are no $\\{\Delta_{x}\\}$ satisfying Eq. (54). Indeed, if there were,
the existence of a right-null vector of ${S}_{x,\varepsilon}$ as defined in
Eq. (46) would imply

$\sum_{\varepsilon>0}\ln b_{\varepsilon}\phi_{\varepsilon}=0\,.$ (66)

However, Eq. (66) cannot hold in general since $\\{b_{\varepsilon}\\}$ is
defined using the symmetric part of the reaction fluxes so encoding kinetic
properties, while $\\{\phi_{\varepsilon}\\}$ is derived from the
stoichiometric matrix so encoding topological properties only.

In summary, the nonequilibrium dynamics of pseudo detailed balanced CRNs, as
well as the resulting self-organized structures, can be exactly obtained by
passive mechanisms too.

###  V.2 Complex Balanced CRNs

CRNs are said to be complex balanced Anderson  _et al._ (2010) when there
exists an homogeneous steady state $\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}$
of the chemical dynamics, i.e.,
$\sum_{\rho>0}{S}_{x,\rho}j_{\rho}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})=0$,
such that the sum of the net currents entering into a complex $k$ equals the
sum of the net currents exiting from it, namely,
$\sum_{\rho>0}\partial_{k,\rho}j_{\rho}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})=0$
$\forall k$, or equivalently

$\sum_{\rho}{\mathbb{1}}_{k,k(\rho)}j_{\rho}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})=0\text{
}\text{ }\text{ }\forall k\,.$ (67)

Deficiency-zero CRNs, where every cycle (i.e., every right-null eigenvector of
the substoichiometric matrix for the internal species) is also a right-null
eigenvector of the incidence matrix, are always complex balanced. If
deficiency is greater than zero, then CRNs are, in general, not complex
balanced (except for a set of kinetic parameters of null measure Cappelletti
and Wiuf (2016); Polettini  _et al._ (2015)).

We now show that the RD dynamics of complex balanced CRNs is driven by the
minimization of a kinetic potential and relaxes to a steady state
$\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ characterized by a vanishing
diffusion entropy production rate, i.e.,
$\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}_{\mathrm{ss}}]=0$, independently
of whether the concentration fields are homogeneous or inhomogeneous. We do so
by generalizing a previous proof for the existence of a Lyapunov function for
homogeneous complex balanced CRNs in dilute solutions undergoing mass-action
kinetics Anderson (2014).

Assumption. We consider reaction fluxes of Arrhenius-like form (48) and we
rewrite them as

$\omega_{\rho}(\boldsymbol{c}(\boldsymbol{r}))=a_{\rho}\,e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}(\boldsymbol{r})){\nu}_{x,\rho}}{RT}}\,,$
(68)

with $a_{\rho}\equiv
A_{\rho}\exp({({\sum_{y}\mu_{y}{\nu}_{y,\rho})}/{RT}})\neq a_{-\rho}$ being a
constant parameter accounting (also) for the chemostatted species. Note that,
since ${\nu}_{x,k(\rho)}={\nu}_{x,\rho}$ by definition of complexes, Eq. (67)
becomes

$\sum_{\rho}{\mathbb{1}}_{k,k(\rho)}\Big{\\{}a_{\rho}-a_{-\rho}e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}){S}_{x,\rho}}{RT}}\Big{\\}}=0\,,$
(69)

when endowed with Arrhenius-like reaction fluxes (68).

Proof. We start by introducing the potential

${F}_{\mathrm{cb}}[\boldsymbol{c}]\equiv
F[\boldsymbol{c}]-\sum_{x}\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})\int_{V}\mathrm{d}\boldsymbol{r}\,c_{x}(\boldsymbol{r})\,,$
(70)

and then we show that ${F}_{\mathrm{cb}}[\boldsymbol{c}]$ acts as a Lyapunov
function of the RD dynamics (14) with respect to a steady state characterized
by homogeneous chemical potentials. Note that the potential (70) is obtained
from the free energy (28) via a shift of the standard chemical potentials of
the internal species $\\{\mu_{x}^{\circ}\\}$ given by

$\Delta_{x}=-\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})\,,$ (71)

like what is done in Eq. (59) for pseudo detailed balanced CRNs.

By using Eqs. (14), (68), and $\mathrm{d}_{t}c_{y}=0$, the time derivative of
${F}_{\mathrm{cb}}[\boldsymbol{c}]$ reads

$\mathrm{d}_{t}{F}_{\mathrm{cb}}[\boldsymbol{c}]=-T{\dot{\Pi}_{\mathrm{cb}}[\boldsymbol{c}]}-T\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}]\,,$
(72)

where $\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}]\geq 0$ is given in Eq.
(37b) (since $\boldsymbol{\nabla}\mu_{y}=0$) , and

$\begin{split}T\dot{\Pi}_{\mathrm{cb}}[\boldsymbol{c}]\equiv-\int_{V}\mathrm{d}\boldsymbol{r}\sum_{\rho,\,x}\Big{\\{}(&\mu_{x}(\boldsymbol{c}(\boldsymbol{r}))-\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}))\\\
&{S}_{x,\rho}a_{\rho}e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}(\boldsymbol{r})){\nu}_{x,\rho}}{RT}}\Big{\\}}\,.\end{split}$
(73)

We have now to show that $-\dot{\Pi}_{\mathrm{cb}}[\boldsymbol{c}]\leq 0$. By
using i) Eq. (2), ii) multiplying and dividing by
$\exp{(\sum_{x}\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}){\nu}_{x,\rho}/{RT})}$,
and iii) $e^{a}(b-a)\leq e^{b}-e^{a}$ for any real number $a$ and $b$ (where
the equality sign is achieved if and only if $a=b$), we obtain

$\small\begin{split}-\frac{\dot{\Pi}_{\mathrm{cb}}[\boldsymbol{c}]}{R}&\leq\int_{V}\mathrm{d}\boldsymbol{r}\sum_{\rho}\Big{\\{}e^{\frac{\sum_{x}(\mu_{x}(\boldsymbol{c}(\boldsymbol{r}))-\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})){\nu}_{x,-\rho}}{RT}}\\\
-&e^{\frac{\sum_{x}(\mu_{x}(\boldsymbol{c}(\boldsymbol{r}))-\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})){\nu}_{x,\rho}}{RT}}\Big{\\}}a_{\rho}e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}){\nu}_{x,\rho}}{RT}}\,,\end{split}$
(74)

where the equality sign is achieved if and only if
$(\mu_{x}(\boldsymbol{c}(\boldsymbol{r}))-\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})){S}_{x,\rho}=0$.
Since $\sum_{\rho}=\sum_{k}\sum_{\rho}{\mathbb{1}}_{k,k(\rho)}$ and
${\nu}_{x,k(\rho)}={\nu}_{x,\rho}$ by definition of complexes, the two
integrands in Eq. (74) can be rewritten as

$\displaystyle\sum_{k}e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}(\boldsymbol{r})){\nu}_{x,k}}{RT}}\sum_{\rho}{\mathbb{1}}_{k,k(\rho)}a_{-\rho}e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}){S}_{x,\rho}}{RT}}\,,$
(75a)
$\displaystyle\sum_{k}e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}(\boldsymbol{r})){\nu}_{x,k}}{RT}}\sum_{\rho}{\mathbb{1}}_{k,k(\rho)}a_{\rho}\,,$
(75b)

respectively, and hence the r.h.s. of Eq. (74) reads

$\begin{split}-\int_{V}\mathrm{d}\boldsymbol{r}\sum_{k}&e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}(\boldsymbol{r})){\nu}_{x,k}}{RT}}\sum_{\rho}{\mathbb{1}}_{k,k(\rho)}\Big{\\{}a_{\rho}\\\
&-a_{-\rho}e^{\frac{\sum_{x}\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}){S}_{x,\rho}}{RT}}\Big{\\}}=0\,,\end{split}$
(76)

because of Eq. (69). Thus,

$\mathrm{d}_{t}{F}_{\mathrm{cb}}[\boldsymbol{c}]\leq 0\,.$ (77)

This, together with the fact that ${F}_{\mathrm{cb}}[\boldsymbol{c}]$ is lower
bounded (since $F[\boldsymbol{c}]$ is lower bounded by assumption, see Subs.
IV.2), implies that the RD dynamics (14) of complex balance CRNs must relax
towards a steady state $\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$
satisfying

$\dot{\Pi}_{\mathrm{cb}}[\boldsymbol{c}_{\mathrm{ss}}]=0\text{ }\text{ }\text{
}\text{ and }\text{ }\text{ }\text{
}\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}_{\mathrm{ss}}]=0\,,$ (78)

or equivalently

$\small\sum_{x}(\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))-\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})){S}_{x,\rho}=0\text{
}\text{ }\text{ }\text{ and }\text{ }\text{ }\text{
}\boldsymbol{J}_{i}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))=0\,,$ (79)

for all $\rho>0$ and for all $i$. $\square$

As we did in Subs. V.1, we now further discuss the conditions, the physical
meaning, and the implications of Eqs. (78) and (79).

First, as for pseudo detailed balanced CRNs, the dynamics of complex balanced
CRNs minimizes ${F}_{\mathrm{cb}}[\boldsymbol{c}]$ in Eq. (70), which plays
the role of a Lyapunov function. However, ${F}_{\mathrm{cb}}[\boldsymbol{c}]$
is not a thermodynamic potential, but rather a kinetic potential: it encodes
kinetic properties via the steady state concentrations in
$\\{\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})\\}$ (see Eq. (70)).

Second, as for pseudo detailed balanced CRNs, Eqs. (78) and (79) physically
mean that diffusion processes are equilibrated at steady state. This, together
with Eq. (37b), implies that the chemical potentials of all species are
homogeneous.

Third, spatial self-organization at steady state does not lead to
inhomogeneous reaction fluxes. Indeed, homogeneous chemical potentials at
steady state, together with reaction fluxes of Arrhenius-like form (68), imply
that the reaction fluxes are homogeneous too, unlike for pseudo detailed
balanced CRNs.

Fourth, as for pseudo detailed balanced CRNs, by plugging
$\boldsymbol{J}_{i}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))=0$ into Eq.
(14), we obtain that both diffusion processes and chemical reactions are at
steady state in every point $\boldsymbol{r}\in V$,. Consequently, no
oscillations or moving structures (e.g., chemical waves) can emerge.

Fifth, as for pseudo detailed balanced CRNs,
$\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ is not, in general, an
equilibrium steady state even if the diffusion processes are equilibrated,
i.e., $\dot{\Sigma}_{\mathrm{dff}}[\boldsymbol{c}_{\mathrm{ss}}]=0$.

Sixth, Eq. (79) implies that
$(\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))-\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}))$
can be written as a linear combination of unbroken conservation laws similarly
to what happens at equilibrium where the equilibrium chemical potentials can
be written as a linear combination of conservation laws (see Eq. (39a) and
App. A).

In summary, the nonequilibrium dynamics of complex balanced CRNs, as well as
the resulting self-organized structures, shares important features with
passive mechanisms, e.g., it is driven by the minimization of (appropriately
constructed) potentials and diffusion processes equilibrate.

###  V.3 Generalization to Non-ideal Chemostatted Species

We show here that the results of Subs. V.1 and V.2 can be generalized to the
case with non-ideal chemostatted species. Namely, when the chemical potentials
$\\{\mu_{y}(\boldsymbol{c}(\boldsymbol{r}))\\}$ (given in Eq. (29) and not in
Eq. (30)) are still homogeneous and constant in time, i.e.,

$\mu_{y}(\boldsymbol{c}(\boldsymbol{r}))=\mu_{y}\quad\forall\boldsymbol{r}\text{
and }\forall t\,,$ (80)

but the concentrations $\\{c_{y}(\boldsymbol{r})\\}$ are not. Note that this
implies that the exchange currents $\\{I_{y}(\boldsymbol{r})\\}$ in Eq. (14)
represent the mechanisms ensuring that Eq. (80) is satisfied, but
$\partial_{t}c_{y}(\boldsymbol{r})\neq 0$ in general.

To do so, we recognize that the exchange currents
$\\{I_{y}(\boldsymbol{r})\\}$ can always be represented as the (net) result of
exchange reactions $\rho_{y}\in\mathcal{R}_{\text{ex}}=\\{\pm 1,\pm
2,\dots,\pm n_{y}\\}$ between the species $y\in Y$ and ideal chemostatted
species $\tilde{y}\in\tilde{Y}$ according to

$Z_{y}\ch{<=>[\rho_{y}][-\rho_{y}]}Z_{\tilde{y}}\,.$ (81)

Indeed, as long as the exchange reactions $\rho_{y}\in\mathcal{R}_{\text{ex}}$
are fast enough (compared to the time scale of the other reactions
$\rho\in\mathcal{R}$ and of the diffusion processes), they are always
equilibrated and, consequently, Eq. (80) is satisfied (with $\mu_{y}$ playing
the role of the chemical potential of $\tilde{y}$). This implies that RD
systems with non-ideal chemostatted species are equivalent to RD systems with
$X\bigcup Y\bigcup\tilde{Y}$ species (where $X\bigcup Y$ are internal species
and $\tilde{Y}$ are ideal chemostatted species) and
$\mathcal{R}\bigcup\mathcal{R}_{\text{ex}}$ reactions. As proven in Subs. V.1
and V.2, the latter will relax to a nonequilibrium steady state with
homogeneous chemical potentials (i.e., where diffusion processes have
equilibrated) if their CRNs are either pseudo detailed balanced or complex
balanced and, therefore, also the former will relax to the same steady state.
Note that the CRN topology of RD systems with $X\bigcup Y\bigcup\tilde{Y}$
species might be different from the one of the original RD systems because of
the introduction of the exchange reactions (81).

###  V.4 Emergence of Spatial Structures by Diffusive Processes

In pseudo detailed balanced and complex balanced CRNs, the steady state
spatial structure $\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ is determined
by the minimization of appropriately constructed kinetic potentials obtained
by a simple shift of the standard chemical potentials of the internal species
$\\{\mu_{x}^{\circ}\\}$. This, together with the equilibration of the
diffusion processes at steady state, suggests that the chemical reactions
determine the total abundances of the chemical species, i.e.,
$\int_{V}\mathrm{d}\boldsymbol{r}\,\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$,
while diffusion determines the spatial structures, i.e., the
$\boldsymbol{r}$-dependence of $\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$.
Indeed, as proved in App. C, the steady-state concentrations
$\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ of pseudo detailed balanced and
complex balanced CRNs can also be obtained by minimizing the free energy of
pure diffusion, $F[\boldsymbol{c}]$ in Eq. (28), when the total abundances are
set to values equal to
$\int_{V}\mathrm{d}\boldsymbol{r}\,\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$.

###  V.5 Final Comment

We notice that the equilibration of diffusion processes in complex balanced
CRNs endowed with Arrhenius-like fluxes was also recently derived in Ref.
Miangolarra and Castellana (2023) with a different approach. Let us now stress
the differences between our work and Ref. Miangolarra and Castellana (2023).
First, our derivation is a closed derivation at the deterministic level of
description of the RD dynamics. Their derivation is based on the stochastic
dynamics of complex balanced CRNs. Second, our derivation is consistent with
nonlocal terms (e.g., of the form
$\boldsymbol{\nabla}c_{x}(\boldsymbol{r})\cdot
K_{x,{x^{\prime}}}\boldsymbol{\nabla}c_{x^{\prime}}(\boldsymbol{r})$, see Eq.
(82) in Sec. VI) in the non-ideal free energy. Their derivation determines the
Lyapunov function of the reaction dynamics (the equivalent of
${F}_{\mathrm{cb}}[\boldsymbol{c}]$ in Eq. (70)) by taking the macroscopic
limit of the steady-state rate function for homogeneous CRNs. Hence, their
Lyapunov function can consistently include only local terms (e.g., of the form
$c_{x}(\boldsymbol{r})M_{x,{x^{\prime}}}c_{x^{\prime}}(\boldsymbol{r})$, see
Eq. (82) in Sec. VI) in the non-ideal free energy (although they use nonlocal
terms in one of the examples). Third, our derivation can be applied to RD
systems with non-ideal chemostatted species (see Subs. V.3). Their derivation
requires ideal chemostatted species.

##  VI Thermodynamic Properties of Active Self-Organization - Numerical
Results

We now examine the thermodynamic properties (introduced in Sec. IV) of
specific RD systems. We first illustrate the exact results obtained in Sec V
for pseudo detailed balanced CRNs and complex balanced CRNs. We then examine a
CRN that is neither pseudo detailed balanced nor complex balanced. To do so,
we will numerically simulate their RD dynamics (14) to determine the final
steady state $\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r})$ starting from a
homogeneous steady state $\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}$ after the
concentration of each internal species $x$ is perturbed with a gaussian white
noise of zero average and $(0.1*c_{x}^{\mathrm{h}})^{2}$ variance. We consider
RD systems in a two-dimensional space with periodic boundary conditions.
Simulations are performed using the py-pde package developed in Ref. Zwicker
(2020).

###  VI.1 Model

We consider RD systems where the set of internal species includes a
nonreacting species, labeled nr in the following, while the chemostatted
species are assumed to be ideal, with homogeneous and time-independent
chemical potentials.

The non-ideal free energy, accounting for the interactions between the
internal species only, is assumed to have the form

$\begin{split}F^{\mathrm{ni}}[\boldsymbol{c}]=\frac{RT}{2}\sum_{x,{x^{\prime}}}\int_{V}\mathrm{d}\boldsymbol{r}\,\big{[}&c_{x}(\boldsymbol{r})M_{x,{x^{\prime}}}c_{x^{\prime}}(\boldsymbol{r})\\\
&+\boldsymbol{\nabla}c_{x}(\boldsymbol{r})\cdot
K_{x,{x^{\prime}}}\boldsymbol{\nabla}c_{x^{\prime}}(\boldsymbol{r})\big{]}\,,\end{split}$
(82)

where $\\{M_{x,{x^{\prime}}}\\}$ and $\\{K_{x,{x^{\prime}}}\\}$ are the
entries of the two (symmetric) matrixes $\mathbb{M}$ and $\mathbb{K}$
describing the mean-field molecular interactions and the cost of forming
interfaces, respectively. Hence, the chemical potentials of the internal
species are given by

$\begin{split}\mu_{x}(\boldsymbol{c}(\boldsymbol{r}))=\mu_{x}^{\circ}&+RT\ln
c_{x}(\boldsymbol{r})\\\
&+RT\sum_{{x^{\prime}}}\big{[}M_{x,{x^{\prime}}}c_{x^{\prime}}(\boldsymbol{r})-K_{x,{x^{\prime}}}\nabla^{2}c_{x^{\prime}}(\boldsymbol{r})\big{]}\,,\end{split}$
(83)

where we also used Eqs. (29) and (32). The chemical potentials of the
chemostatted species $\\{\mu_{y}\\}$ are treated as control parameter like in
Sec. V and are related to the corresponding concentrations $\\{c_{y}\\}$ via
Eq. (30).

Reaction fluxes are assumed to follow Arrhenius-like forms given in Eq. (68).
The entries of the Onsager matrix $\mathbb{O}(\boldsymbol{c}(\boldsymbol{r}))$
in the diffusion currents (36) are assumed to read

$O_{i,j}(\boldsymbol{c}(\boldsymbol{r}))=D_{i}c_{i}(\boldsymbol{r}){\mathbb{1}}_{i,j}\,,$
(84)

where $D_{i}$ is the diffusion coefficient of species $i$. Hence, the
diffusion currents for the internal species are given by

$\boldsymbol{J}_{x}(\boldsymbol{c}(\boldsymbol{r}))=-D_{x}c_{x}(\boldsymbol{r})\boldsymbol{\nabla}\mu_{x}(\boldsymbol{c}(\boldsymbol{r}))\text{
}\text{ }\forall x\in X\,.$ (85)

###  VI.2 Pseudo Detailed Balanced CRN

We consider a RD system with three internal species
$\\{\ch{X_{1}},\ch{X_{2}},\ch{X_{\text{nr}}}\\}$ and two chemostatted species
$\\{\ch{Y_{1}},\ch{Y_{2}}\\}$ which are interconverted via the reactions

$\begin{split}X_{1}+Y_{1}&\ch{<=>[$1$][$-1$]}X_{2}+Y_{2}\,,\\\
X_{1}&\ch{<=>[$2$][$-2$]}X_{2}\,.\end{split}$ (86)

The CRN (86) has three complexes $V_{1}=\ch{X_{1}}$, $V_{2}=\ch{X_{2}}$, and
$V_{3}=\ch{X_{\text{nr}}}$ and the corresponding graph of complexes reads

$\begin{split}V_{1}&\ch{<=>[$1$][$-1$]}V_{2}\,,\\\
V_{1}&\ch{<=>[$2$][$-2$]}V_{2}\,,\end{split}$ (87)

admitting only the lumped reaction from $V_{1}$ to $V_{2}$ and its backward
counterpart from $V_{2}$ to $V_{1}$. Hence, the stoichiometric matrix for the
internal species,

${\mathbb{S}_{X}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt\cr}}}}\right)$}}\,,$ (88)

becomes the matrix

$\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,$ (89)

when considering the only lumped reaction. Since the matrix in Eq. (89) does
not admit any right-null eigenvector, the CRN (86) is pseudo detailed balanced
(see Subs. V.1).

Figure 5: Steady state and corresponding thermodynamic properties of the
pseudo detailed balanced CRN (86). (a) Concentration field of $\ch{X_{2}}$.
(b) Density of reaction entropy production rate. (c) Difference between the
chemical potential of $i\in\\{\ch{X_{1}},\ch{X_{2}}\\}$ (labeled here
$\mu_{i}^{\text{ss}}(\boldsymbol{r})$) and its chemical potential at the
homogeneous steady state (labeled here $\mu_{i}^{\text{h}}$). (d) Net current
of the lumped reaction, i.e.,
$\hat{j}_{1}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))=j_{1}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))+j_{2}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))$.
Simulation parameters in arbitrary units: $RT=1$, $\chi=3$, $k_{1}=0.1$,
$k_{2}=0.05$, $\mu_{1}^{\circ}=0$, $\mu_{2}^{\circ}=-2$, $\mu_{\ch{Y}_{1}}=1$,
$\mu_{\ch{Y}_{2}}=0$, $A_{\rho}=1$ $\forall\rho$, $D_{x}=1$ $\forall x$,
$c_{\text{1}}^{\mathrm{h}}=0.087$, $c_{\text{2}}^{\mathrm{h}}=0.213$,
$c_{\text{nr}}^{\mathrm{h}}=0.7$.

We solve the RD dynamics (14) when the mean-field molecular interactions and
the cost of forming interfaces are given by the following matrixes

$\displaystyle\mathbb{M}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$\chi$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\chi$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\chi$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$\chi$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,\quad$
$\displaystyle\mathbb{K}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$k_{1}$\hfil\kern 5.0pt&5.0pt\hfil$k_{2}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$k_{2}$\hfil\kern 5.0pt&5.0pt\hfil$k_{1}$\hfil\kern
5.0pt&5.0pt\hfil$k_{2}$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$k_{2}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,.$ (98)

The typical steady-state solution and its thermodynamic properties are shown
in Fig. 5. We observe that the concentration fields spatially organize (Fig.
5a) into a pattern that is qualitatively analogous to a complete phase
separation obtained at equilibrium: the species X_2 accumulates in a single
droplet. The chemical reactions are however out of equilibrium, i.e.,
$\dot{\sigma}^{\mathrm{ss}}_{\mathrm{rct}}(\boldsymbol{r})>0$ (see Fig. 5b),
even if the net current along the lumped reaction, namely,
$\hat{j}_{1}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))=j_{1}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))+j_{2}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))$,
vanishes as predicted in Subs. V.1 (see Fig. 5d). Furthermore, in agreement
with our derivation in Subs. V.1, at steady state diffusion processes are at
equilibrium, i.e.,
$\dot{\sigma}^{\mathrm{ss}}_{\mathrm{dff}}(\boldsymbol{r})=0$, because the
chemical potentials are homogeneous (see Fig. 5c).

Note that the steady-state density of reaction entropy production rate in Fig.
5b is homogeneous because we used reaction fluxes of Arrhenius-like forms,
and, therefore the reactions fluxes become homogeneous at steady state (as
explained in Subs. V.1).

###  VI.3 Complex Balanced CRN

We consider a RD system with four internal species
$\\{\ch{X_{1}},\ch{X_{2}},\ch{X_{3}},\ch{X_{\text{nr}}}\\}$ and two
chemostatted species $\\{\ch{Y_{1}},\ch{Y_{2}}\\}$ which are interconverted
via the reactions

$\begin{split}X_{1}+Y_{1}&\ch{<=>[$1$][$-1$]}X_{2}+Y_{2}\,,\\\
X_{2}&\ch{<=>[$2$][$-2$]}X_{3}\,,\\\
X_{3}&\ch{<=>[$3$][$-3$]}X_{1}\,.\end{split}$ (99)

The CRN (99) has four complexes $V_{1}=\ch{X_{1}}$, $V_{2}=\ch{X_{2}}$,
$V_{3}=\ch{X_{3}}$, and $V_{4}=\ch{X_{\text{nr}}}$ and the corresponding graph
of complexes reads

$\begin{split}V_{1}&\ch{<=>[$1$][$-1$]}V_{2}\,,\\\
V_{2}&\ch{<=>[$2$][$-2$]}V_{3}\,,\\\
V_{3}&\ch{<=>[$3$][$-3$]}V_{1}\,.\end{split}$ (100)

Since the substoichiometric matrix for the internal species

${\mathbb{S}_{X}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,$ (101)

admits only one cycle that is also a right-null eigenvector of the incidence
matrix,

${\partial}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,$ (102)

the CRN (99) is deficiency-zero (see Eq. (11)) and, therefore, complex
balanced (see Subs. V.2).

Figure 6: Steady state and corresponding thermodynamic properties of the
complex balanced CRN (99). (a) Concentration field of $\ch{X_{2}}$. (b)
Density of reaction entropy production rate. (c) Difference between the
chemical potential of $i\in\\{\ch{X_{1}},\ch{X_{2}},\ch{X_{3}}\\}$ (labeled
here $\mu_{i}^{\text{ss}}(\boldsymbol{r})$) and its chemical potential at the
homogeneous steady state (labeled here $\mu_{i}^{\text{h}}$). (d) Difference
between the chemical potential of X_nr (labeled here
$\mu_{\text{nr}}^{\text{ss}}(\boldsymbol{r})$) and its chemical potential at
the homogeneous steady state (labeled here $\mu_{\text{nr}}^{\text{h}}$).
Simulation parameters in arbitrary units: $RT=1$, $\chi=3$, $k_{1}=0.5$,
$k_{2}=0.1$, $\mu_{1,2}^{\circ}=1$, $\mu_{3}^{\circ}=3$,
$\mu_{\text{nr}}^{\circ}=0$, $\mu_{Y_{1}}=4$, $\mu_{Y_{2}}=0$, $A_{\rho}=1$
$\forall\rho$, $D_{x}=1$ $\forall x$, $c_{\text{1}}^{\mathrm{h}}=0.007$,
$c_{\text{2}}^{\mathrm{h}}=0.274$, $c_{\text{3}}^{\mathrm{h}}=0.019$,
$c_{\text{nr}}^{\mathrm{h}}=0.7$.

We solve the RD dynamics (14) when the mean-field molecular interactions and
the cost of forming interfaces are given by the following matrixes

$\displaystyle\mathbb{M}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$\chi$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\chi$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$\chi$\hfil\kern 5.0pt&5.0pt\hfil$\chi$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$\chi$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$\chi$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt\cr}}}}\right)$}}\,,$ $\displaystyle\mathbb{K}=\hbox{}\vbox{\kern
0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$k_{1}$\hfil\kern 5.0pt&5.0pt\hfil$k_{2}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$k_{2}$\hfil\kern 5.0pt&5.0pt\hfil$k_{1}$\hfil\kern
5.0pt&5.0pt\hfil$k_{2}$\hfil\kern 5.0pt&5.0pt\hfil$k_{2}$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$k_{2}$\hfil\kern
5.0pt&5.0pt\hfil$k_{1}$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$k_{2}$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$k_{1}$\hfil\kern
5.0pt\cr}}}}\right)$}}\,.$ (113)

The typical steady-state solution and its thermodynamic properties are shown
in Fig. 6. As for the RD system in Subs. VI.2, we observe that the
concentration fields spatially organize (Fig. 6a) into a pattern that is
qualitatively analogous to a complete phase separation obtained at
equilibrium: the species X_2 accumulates in a single droplet. The chemical
reactions are however out of equilibrium, i.e.,
$\dot{\sigma}^{\mathrm{ss}}_{\mathrm{rct}}(\boldsymbol{r})>0$ (see Fig. 6b),
while diffusion processes are at equilibrium, i.e.,
$\dot{\sigma}^{\mathrm{ss}}_{\mathrm{dff}}(\boldsymbol{r})=0$, since the
chemical potentials are homogeneous (see Fig. 6c) as predicted by our
derivation in Subs. V.2.

Furthermore, according to Eq. (79), the difference
$(\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))-\mu_{x}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}))$
for every internal species can be written as a linear combination of unbroken
conservation laws (defined in Subs. III.2). This implies that in general

$\begin{split}(\mu_{\ch{X_{1}}}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))-\mu_{\ch{X_{1}}}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}))&=(\mu_{\ch{X_{2}}}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))-\mu_{\ch{X_{2}}}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}))\\\
&=(\mu_{\ch{X_{3}}}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))-\mu_{\ch{X_{3}}}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}))\\\
&\neq(\mu_{\ch{X_{\text{nr}}}}(\boldsymbol{c}_{\mathrm{ss}}(\boldsymbol{r}))-\mu_{\ch{X_{\text{nr}}}}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}))\,,\end{split}$
(114)

since the stoichiometric matrix for the internal species (101) admits the
following left null eigenvectors

$\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt\cr}}}}\right)$}}\,,\quad\quad\quad\hbox{}\vbox{\kern
0.86108pt\hbox{$\kern 0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr}}}}\right)$}}\,.$ (115)

The constraint in Eq. (114) is confirmed by the simulation (compare Fig. 6(c)
and Fig. 6(d)).

###  VI.4 Brusselator

We consider a RD system of three internal species
$\\{\ch{X_{1}},\ch{X_{2}},\ch{X_{\text{nr}}}\\}$ and four chemostatted species
$\\{\ch{Y_{1}},\ch{Y_{2}},\ch{Y_{3}},\ch{Y_{4}}\\}$ which are interconverted
via the reactions

$\begin{split}Y_{1}&\ch{<=>[$1$][$-1$]}X_{1}\,,\\\
2X_{1}+X_{2}&\ch{<=>[$2$][$-2$]}3X_{1}\,,\\\
Y_{2}+X_{1}&\ch{<=>[$3$][$-3$]}X_{2}+Y_{4}\,,\\\
X_{1}&\ch{<=>[$4$][$-4$]}Y_{3}\,.\end{split}$ (116)

The CRN (116) is known as Brusselator and can be represented in terms of the 6
complexes $V_{1}=\ch{X_{1}}$, $V_{2}=\ch{X_{2}}$, $V_{3}=\ch{2X_{1}+X_{2}}$,
$V_{4}=\ch{3X_{1}}$, $V_{5}=\ch{X_{\text{nr}}}$, and $\emptyset=Y_{1,2}$ and
the following graph of complexes:

$\begin{split}\emptyset&\ch{<=>[$1$][$-1$]}V_{1}\,,\\\
V_{3}&\ch{<=>[$2$][$-2$]}V_{4}\,,\\\ V_{1}&\ch{<=>[$3$][$-3$]}V_{2}\,,\\\
V_{1}&\ch{<=>[$4$][$-4$]}\emptyset\,.\end{split}$ (117)

The stoichiometric matrix for the internal species

${\mathbb{S}_{X}}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt\cr}}}}\right)$}}\,$ (118)

admits 2 cycles, i.e.,

$\boldsymbol{\phi}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr}}}}\right)$}}\quad\text{ and
}\quad\boldsymbol{\phi}^{\prime}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern
0.0pt\kern 2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle$\hfil\kern 5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr}}}}\right)$}}\,,$ (119)

but only the first one is a right-null vector of the incidence matrix

${\partial}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left(\kern
0.0pt\kern-2.5pt\kern-6.66669pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt\crcr\kern 0.86108pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$-1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt\cr
5.0pt\hfil\@arstrut$\scriptstyle\color[rgb]{0.50,0.50,0.50}\definecolor[named]{pgfstrokecolor}{rgb}{0.50,0.50,0.50}\pgfsys@color@gray@stroke{0.50}\pgfsys@color<EMAIL_ADDRESS>5.0pt&5.0pt\hfil$-1$\hfil\kern 5.0pt&5.0pt\hfil$0$\hfil\kern
5.0pt&5.0pt\hfil$0$\hfil\kern 5.0pt&5.0pt\hfil$1$\hfil\kern
5.0pt\cr}}}}\right)$}}\,.$ (120)

Hence, the deficiency reads $\delta=1$ (see Eq. (11)) and, in general, the CRN
(116) is neither pseudo detailed balanced nor complex balanced. Namely, any
homogeneous steady state $\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}}$ of the
chemical dynamics satisfying
$\sum_{\rho>0}{S}_{x,\rho}j_{\rho}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})=0$
does not satisfy
$\sum_{\rho>0}\partial_{k,\rho}j_{\rho}(\boldsymbol{c}_{\mathrm{ss}}^{\mathrm{h}})=0$
too (except for a set of parameters, e.g., $\\{A_{\rho}\\}$, $\\{\mu_{y}\\}$,
$\mathbb{M}$, and $\mathbb{K}$, of null measure Cappelletti and Wiuf (2016);
Polettini  _et al._ (2015)).

Figure 7:  Steady state concentration fields of (a) $\ch{X_{1}}$ and (b)
$\ch{X_{2}}$ of the CRN (116). Simulation parameters in arbitrary units:
$RT=1$, $\chi=0$, $k_{1}=0.1$, $k_{2}=0.05$, $\mu_{x}^{\circ}=0$, $\forall x$,
$\mu_{Y_{1}}=\ln 2$, $\mu_{Y_{2}}=\ln 5$, $\mu_{Y_{3}}=\ln 0.01$,
$\mu_{Y_{4}}=\ln 0.1$, $A_{\rho}=1$ $\forall\rho$, $D_{1}=1$, $D_{2}=10$,
$D_{\text{nr}}=1$, $c_{\text{1}}^{\mathrm{h}}=1.05$,
$c_{\text{2}}^{\mathrm{h}}=5.76$, $c_{\text{nr}}^{\mathrm{h}}=1.00$.  Figure
8:  Steady-state thermodynamic properties of the CRN (116). (a) Chemical
potential of $\ch{X_{1}}$ (b) Density of entropy production rate
$\dot{\sigma}^{\mathrm{ss}}(\boldsymbol{r})=\dot{\sigma}^{\mathrm{ss}}_{\mathrm{dff}}(\boldsymbol{r})+\dot{\sigma}^{\mathrm{ss}}_{\mathrm{rct}}(\boldsymbol{r})$.
(c) Density of diffusion entropy production rate
$\dot{\sigma}^{\mathrm{ss}}_{\mathrm{dff}}(\boldsymbol{r})$. (d) Density of
reaction entropy production rate
$\dot{\sigma}^{\mathrm{ss}}_{\mathrm{rct}}(\boldsymbol{r})=\sum_{\rho=1}^{4}\dot{\sigma}^{\mathrm{ss}}_{\rho}(\boldsymbol{r})$
and (e)-(f)-(g)-(h) its decomposition into the contribution of each reaction
$\rho=1,\,2,\,3,\,4$ as explained in Subs. IV.2. Simulation parameters are the
same as in Fig. 7.

We solve the RD dynamics (14) when the mean-field molecular interactions and
the cost of forming interfaces are given by the same matrixes as in Eq. (98).
The typical steady-state solution and its thermodynamic properties are shown
in Figs. 7 and 8, respectively. We observe that the concentration fields
spatially organize (Fig. 7) into a pattern that is striking different from the
previous ones examined in Subs. VI.2 and VI.3 (see Figs. 5 and 6,
respectively). In this case, X_1 accumulates in multiple droplets and X_2
remains almost homogeneous. Furthermore, both reactions and diffusion
processes are out of equilibrium since the chemical potentials are not
homogeneous (Fig. 8a):
$\dot{\sigma}^{\mathrm{ss}}_{\mathrm{rct}}(\boldsymbol{r})>0$ (Fig. 8d) and
$\dot{\sigma}^{\mathrm{ss}}_{\mathrm{dff}}(\boldsymbol{r})>0$ (Fig. 8c). This
has two main consequences. First, the dissipation due to the chemical
reactions, quantified by $\dot{\sigma}_{\mathrm{rct}}(\boldsymbol{r})$, is not
homogeneous either: reactions dissipate more inside the droplets (Fig. 8d).
However, not all reactions contribute in the same way to the total reaction
entropy production rate. Indeed, as show in Fig. 8e-h, reaction $\rho=1$
dissipates more outside the droplets, while $\rho=2$, $3$, and $4$ dissipate
more inside the droplets. Furthermore, $\rho=2$ (representing an autocatalytic
reaction) is the reaction dissipating the most. Second, the diffusion
processes are out of equilibrium and dissipate mostly at the interface between
the droplets and the bulk as shown in Fig. 8c.

Note that the steady-state density of nonconservative work (which satisfies
$\dot{w}^{\mathrm{ss}}_{\mathrm{nc}}(\boldsymbol{r})=\dot{\sigma}^{\mathrm{ss}}(\boldsymbol{r})$)
is peaked inside the droplets (see Fig. 8b) and, therefore, most of the free
energy maintaining nonequilibrium self-organization is provided from inside
the droplets (even if the chemostatted species are homogeneous).

##  VII Conclusions

We provided a general thermodynamic framework to study the energetics of self-
organization in open non-ideal RD systems. We showed that it can be used to
study in great detail, with spatial resolution, the contribution of the
reactive and diffusive processes to the dissipation of self-organized
structures. For special classes of RD systems, our framework allowed us to
construct potentials that are minimized by the dynamics (Lyapunov functions)
and whose minima correspond to the steady state. Thermodynamic potentials
(independent of the kinetics) are found for detailed balance RD systems that
relax to equilibrium, while kinetic potentials are found for pseudo detailed
balanced and complex balanced CRNs which relax to nonequilibrium steady
states. In doing so, we demonstrated that thermodynamics is useful not only to
study the energetics of RD systems but also to make inferences on their
dynamics.

We emphasize that our local detailed balance condition (35) is in agreement
with the general derivation in Ref. de Groot and Mazur (1984), where the
thermodynamic forces conjugated to the reaction currents are (stoichiometric-
weighted combinations of) chemical potentials
$\\{\mu_{i}(\boldsymbol{c}(\boldsymbol{r}))\\}$. Other works have used instead
the exchange chemical potentials
$\\{\mu_{i}^{\mathrm{ex}}(\boldsymbol{c}(\boldsymbol{r}))\\}$, namely, the
differences between the chemical potentials and the chemical potential of
another species chosen as reference. This is only consistent with the local
detailed balance condition when dealing with unimolecular (or pseudo-
unimolecular) reactions.

Future studies will be needed to extend our framework to fluctuating RD
systems Tiani and Täuber (2023) and to explore the phenomenology of self-
organized structures in non-ideal RD systems following the line of Refs Menou
_et al._ (2023); Luo and Zwicker (2023). Another interesting extension
concerns microswimmers which undergo a self-propelled motion Gaspard and
Kapral (2019); Markovich  _et al._ (2021); Fodor  _et al._ (2022); Bebon  _et
al._ (2024); Chatzittofi  _et al._ (2024); here, the challenge is to
understand how to properly describe the coupling between chemical reactions,
occurring either at the surface of (e.g., Janus colloids in a fuel bath) or
within (e.g., bacteria fulled by nutrient consumption) microswimmers, and the
resulting self-propelled dynamics.

##  VIII Acknowledgments

FA is supported by the project P-DiSC#BIRD2023-UNIPD funded by the Department
of Chemical Sciences of the University of Padova (Italy). TA, ÉF and ME are
supported by the Fond National de la Recherche—FNR, Luxembourg: T.A. by the
Project ThermoElectroChem (C23/MS/18060819), ÉF by the Project SMAC
(14389168), and ME by the Project ChemComplex (C21/MS/16356329).

##  Appendix A Minimization of the Free Energy

We prove here that the equilibrium steady state
$\boldsymbol{c}_{\mathrm{eq}}(\boldsymbol{r})$ of closed RD systems (defined
in Eq. (39)) corresponds to the minimum $\boldsymbol{c}^{*}(\boldsymbol{r})$
of the free energy (28) satisfying the conserved quantities of the RD dynamics
(14) (i.e., the abundances of the moieties (19) as discussed in Subs. III.2).

We start by recognizing that $\boldsymbol{c}^{*}(\boldsymbol{r})$ is defined
by the minimum $(\boldsymbol{c}^{*}(\boldsymbol{r}),\boldsymbol{f}^{*})$ of
the Lagrangian functional

$\small\begin{split}\mathcal{L}[\boldsymbol{c},\boldsymbol{f}]{}=F[\boldsymbol{c}]&-\sum_{\lambda}f_{\lambda}\Big{(}\sum_{i}{\ell}^{\lambda}_{i}\int_{V}\mathrm{d}\boldsymbol{r}\,c_{i}(\boldsymbol{r})-{L}^{\lambda}\Big{)}\,,\end{split}$
(121)

with $\boldsymbol{f}=(\dots,f_{\lambda},\dots)$ being the vector of the
Lagrange multipliers and $\\{{L}^{\lambda}\\}$ being the actual values of the
moiety abundances. Namely,
$(\boldsymbol{c}^{*}(\boldsymbol{r}),\boldsymbol{f}^{*})$ is the solution of
the following equations

$\displaystyle\frac{\delta\mathcal{L}[\boldsymbol{c},\boldsymbol{f}]}{\delta
c_{i}(\boldsymbol{r})}\bigg{|}_{*}$