On non-ideal chemical-reaction networks and phase separation

Much of the theory on chemical-reaction networks (CRNs) has been developed in the ideal-solution limit, where interactions between the solutes are negligible. However, there is a large variety of phenomena in biological cells and soft-matter physics which appear to deviate from the ideal-solution behaviour. Particularly striking is the case of liquid-liquid phase separation, which is typically caused by inter-particle interactions. Here, we revisit a number of known results in the domain of ideal CRNs, and we generalise and adapt them to arbitrary interactions between the solutes which stem from a given free energy. Among these is the form of the steady-state probability distribution and Lyapunov functions for complex-balanced networks, where the creation and annihilation rates are equal for all chemical complexes which appear as reactants or products in the CRN. Finally, we draw a phase diagram for complex-balanced reaction-diffusion solutions based on the minimisation of such Lyapunov function with a rationale similar to that of equilibrium thermodynamics, but for systems that may sustain non-equilibrium chemical currents at steady state. Nevertheless, we find that complex-balanced networks are not sufficient to create diffusion currents at steady state.


I. INTRODUCTION
The cytoplasm of a cell does not behave like an ideal solution [1], since, in many cases, interactions among the solutes cannot be neglected. Indeed, in the cytoplasm there is a plethora of interactions among proteins, other macromolecules, and ions. Some of the most common interactions that are relevant in the cellular cytoplasm are steric and crowding effects [2,3], as well as electrostatic interactions [4,5]. Arguably, the most striking phenomenon caused by these interactions is the emergence of phase-separated condensates, also known as membraneless organelles in the cell-biology literature, which are now widely studied [6][7][8]. The composition of these membraneless organelles is different from the one of the cytoplasm, because they are typically enriched in a specific type of molecules while they exclude others [9]. Moreover, it has been hypothesised that such organelles spatially con-trol biochemical reactions, by modulating their rates and specificity within the condensate [10][11][12][13].
Given their important role in the internal spatial organisation of cells, the regulation of phase-separation phenomena is crucial for many cellular functions. One of the ways in which cells can dynamically control the onset, composition and function of membraneless organelles is through chemical reactions, notably post-translational modifications like phosphorylation [14,15] or methylation [16]. However, phase separation is also triggered by changes in the environment [17,18], thus establishing biological condensates as potential switch-like sensing and regulatory mechanisms.
While most of the insights outlined above are the result of extensive experimental efforts, the interplay between interactions within the solution and non-equilibrium chemical reactions has also been widely studied from the theoretical standpoint. Most of these efforts [19][20][21][22] have been based on effective reaction-diffusion models that can describe patterning and non-equilibrium phenomena in a simple way, but lack thermodynamic consistency. More precisely, in these approaches the reaction dynamics is modelled with mass-action kinetics (MAK), which implicitly assumes that the solution is ideal (while the interaction-influenced diffusion that drives phase-separation is not), which leads to the aforementioned lack of consistency. Early progress in reconciling the spatial patterns predicted by these models with a thermodynamically consistent description was limited to a linear-stability analysis of binary systems [23]. More recently, some works aimed at establishing a deterministic theory for non-ideal chemical-reaction networks (CRNs) [24], the relation between phase coexistence and chemical kinetics [25], and exploring minimal examples for pattern formation with non-ideal CRNs [26,27]. Nevertheless, the link between non-equilibrium CRNs and phase separation has not yet been elucidated in full generality.
Here, we aim at building a thermodynamically consistent framework for interacting reaction-diffusion systems which may exhibit phase-separation at steady state. Therefore, in this framework, in the same way diffusion is governed by a free energy (that takes into account the interactions), the dynamics of the chemical reactions must also reflect this freeenergetic dependency. Here, previous efforts are complemented by analysing the behaviour of non-ideal CRNs in the stochastic limit, in an effort to build a complete theory. This explicit description of non-ideality in the CRN allows us to naturally adapt and generalise the results from the well-established theory of ideal CRNs. We do so by first constructing a framework such that, in the absence of explicit non-equilibrium driving, the system relaxes to thermodynamical equilibrium. Then, we focus on complex-balanced networks, for which the steady-state creation and annihilation rate of each chemical complex are equal. For these type of CRNs, we derive the steady-state probability distributions and Lyapunov functionals, which allows us to obtain the steady-state concentration profiles.
The paper is organised as follows: In Section II we describe the dynamics of spatially homogeneous CRNs in the stochastic and deterministic limit, introduce the concept of complex balance, and recall the main features of MAK. In Section III, we impose thermodynamical constraints on the reaction rates for CRNs at equilibrium, by consistently relating these rates to the free energy, and discuss how they can be modified in non-equilibrium settings. In Section IV we generalise to non-ideal CRNs the known result for the steady-state distribution of complex-balanced networks. Building on this result, in Section V we propose a candidate Lyapunov function of complex-balanced systems. In the same Section, we generalise the previous Lyapunov function to systems with spatial inhomogeneities, and derive the resulting phase diagram for a non-equilibrium, complex-balanced, chemically reactive mixture. Finally, in Section VI we discuss the interpretations and implications of our results.

II. CHEMICAL-REACTION NETWORKS
A chemical-reaction network (CRN) is composed of N chemical species and M reaction pathways; which we assume reversible, for a better alignment with thermodynamic principles. A reaction within the CRN, denoted by the label ρ, is specified as follows: where X a , a = 1, · · · , N is one the N species in the network. In the rest of this paper, the indexes a and b will be used for chemical species only. The matrices r ρ a and s ρ a denote the number of particles of each species participating in the forward and backward reaction, respectively, i.e., r ρ a specifies the number of reactants of type a in the forward reaction ρ, and s ρ a , that of the products of type a in the backward reaction. Note that, given that the reactions are taken to be reversible, the distinction between reactants and products is arbitrary.
The amount of particles of species a created along the forward reaction is denoted by v ρ a = s ρ a − r ρ a . (2) We also introduce the vectors ρ = (v ρ 1 , · · · , v ρ N ), r ρ = (r ρ 1 · · · r ρ N ) and s ρ = (s ρ 1 , · · · , s ρ N ). , and the matrices Finally, we define a complex z as the number and type of particles that participate in a chemical reaction as either reactants (z = r ρ ) or products (z = s ρ ). A single complex z may appear in more than one reaction within the network. In order to clarify the definitions above, we illustrate them for the following network.

A. Example
For the CRN the vectors of reactants are r 1 = 1, 1, 0, 0 , r 2 = 0, 1, 0, 0 , and the product vectors are s 1 = 0, 0, 1, 0 , s 2 = 0, 0, 0, 1 , where each of the entries in the vector correspond to different species and the reactions are labelled by the superindices. Finally, we have the vectors that specify the net amount of particles of each species created by each of the reactions occurring once in the forward direction.
Building on the definitions above, in what follows we introduce the stochastic and deterministic description of a CRN, the starting point of the rest of this work.

B. Stochastic description
If the chemical species in the solution diffuse fast (with respect to the typical timescale of chemical reactions) and is stirred regularly, the system may be considered to be well mixed and it can be described in terms of a single homogeneous concentration of each of the species across space. Then, a state of the system-the number of particles of each type-is determined by the vector n = (n 1 , · · · , n N ).
Each state n of the system has a probability measure P (n, t) at any instant of time t. The dynamics for the probability of states of homogeneous CRNs is given by the Chemical Master Equation (CME), which reads [28]: where the summations over ρ run over all reactions in the CRN, v is given by Eq. (2), and the rate of the transitions in the network is given by the propensity functions f ±ρ , +ρ corresponding to the forward direction of the reaction and −ρ corresponding to the backward direction. Eq. (10) is different to other statements of the CME because we made explicit the fact that every reaction is reversible.

C. Deterministic description
For large particle numbers, by averaging both sides of Eq. (10) and assuming vanishing correlations according to the mean-field picture-assumptions that are supposed to be accurate in the large-particle number limit-one can work out a set of equations for the concentrations c a = n a /V in the macroscopic limit, where both n a and V are large. In this limit, the state of the system is specified by the concentrations c 1 , · · · , c N , and one obtains the following classical set of equations for the dynamics of the concentrations in a CRN [29]: where the currents J still need to be determined. While both the deterministic and stochastic descriptions refer to the same system, the former one is only accurate for large particle numbers, also known as the thermodynamic limit, where fluctuations are negligible.

D. Complex balance
In a CRN, a complex is a set of chemical species and their respective particle numbers, which take part in a reaction, as either reactants or products. Its most general expression is the vector z = (z 1 , · · · , z a , · · · ), where the index a runs over all chemical species and the integer z a is the number of molecules of the species X a in the complex z. Any CRN can be represented as a graph whose nodes denote the complexes that take part in the reactions, where there is an edge between two complexes if and only if there is a reaction z m z n in the CRN, and z m and z n are two different complexes-see for example   Fig. 1 as example. A) Detailed-balanced steady state. Here, the backward and forward rates of each reaction are equal. B) Complex-balanced steady state. There exist non-vanishing net currents (red arrows) but the steady state still satisfies the complex-balanced requirement, Eq. (15). It only allows for cycles at steady state that can be visualised from the network in terms of complexes (red cycle). C) General steady state. There exist non-vanishing net currents (red arrows) and the complex-balanced requirement does not hold. state which is not at equilibrium (i.e. is not entirely detailed-balanced) is by taking the net rate J +ρ − J −ρ in reaction 1 equal to the net rate of those of reactions 2 and 3, and thus reactions 4 and 5 must be detailed balanced (since the system has to be at steady state). This steady state is depicted in Fig. 2

B.
General steady state: The most general class of steady states is defined by the vanishing time derivatives of the dynamical equation (11). By splitting the contributions of each complex z, this condition can be rewritten as m z m a ρ|s ρ =zm where m is an index that labels each of the complexes in the network and the integer z m a represents its components. As a result, there are no constraints between the net rates of each reaction other than those imposed by the stationarity condition of Eq. (11).
In the example of Fig. 1, this implies that, at steady state, there can be current cycles where, for example, species A is created by reaction 1, but annihilated by reaction 4 through the complex A + D, thus breaking complex balance (see Fig. 2 C).
From this hierarchical classification, it can be clearly seen that detailed balance, Eq. (16), implies complex balance, Eq. (17), which, in turn, implies the steady-state condition, Eq. (18). However, the converse is not true: A general steady state is not necessarily complex balanced, and a complex-balanced steady state is not necessarily detail balanced. Therefore, complex balance is less restrictive of a constraint than detailed balance, but it is more restrictive than a generic steady state.

E. Complex balance in networks with mass-action kinetics
As a particular instance of special importance, in what follows we will discuss complex balance in ideal CRNs with MAK.
In short, MAK consists of the hypothesis that the rate of the chemical reaction is proportional to the product of the concentrations of the reactants: As a result, in the deterministic description, the MAK expressions for the currents read where k ±ρ are the rate constants. In what follows, we will denote by c * a the steady-state concentration of species a in the deterministic description. Then, for a deterministic system with MAK, the complex-balance condition is given by Conversely, in the stochastic description [with dynamics is given by Eq. (10)] the MAK expressions for f ±ρ become Previous studies [33] have shown that the steady state of complex-balanced CRNs with MAK is known to have a product-form expression in terms of independent Poisson distributions, and it reads linking the deterministic steady state (c * a ) to the stochastic steady state π. Furthermore, in Ref. [33] it is shown that, if the propensity functions take the more general form then the steady-state distribution reads where M is a normalisation constant and θ a function which maps the vector of integer numbers n into a real-valued positive number.
In what follows, we will demonstrate that the result (23) can be generalised to the nonideal case, i.e., to a class of propensity functions f which take into account the physical interactions between molecules.

NETWORKS
In the previous Section we introduced the general description of CRNs, both on a stochastic and deterministic level: in either cases, a choice for the propensity functions, or currents, must be made to set the network dynamics. For ideal solutions, the most common choice is MAK, as outlined above. However, in what follows we consider solutes which mutually interact and which are, therefore, not ideal, and specify the propensity functions.

A. Equilibrium systems
Here, we consider CRNs at thermodynamic equilibrium, i.e., systems which are not subject to external, non-equilibrium driving. Given that the system is an equilibrium one, at steady state the principle of detailed balance must hold for every reaction ρ: The probability flux across a reaction ρ in the forward direction must equal the one in the backward direction. In this Section, we will impose the detailed-balance condition on the propensity functions at thermal equilibrium and suggest a generalisation for systems out of equilibrium.

Stochastic description
In the stochastic description, the detailed-balance condition at steady state reads where the equilibrium probability distribution P eq (n) for closed stochastic systems-total number of particles fixed-is given by the canonical Boltzmann distribution: with β = 1/(k B T ), k B is the Boltzmann constant, T the temperature, F (n) the Helmholtz free energy of the system in state n, and Z a normalisation factor-the partition function in statistical physics. For systems that exchange mass with a single particle reservoir, the equilibrium distribution (28) is replaced by the distribution for the grand-canonical ensemble [35]. Combined with Eq. (28), the detailed-balance condition in Eq. (27) yields the following constraint for the propensity functions: Then, we choose the following functional form for the propensity functions: where k ρ is the reaction constant, which needs to be equal in both the forward and the backward reaction for Eq. (29) to be satisfied. Given that the free energy F may, in general, depend on the inter-particle interactions-such as steric, electrostatic, or other interactions-Eq. (29) implies that the chemical-reaction rates may depend on these interparticle interactions. The choices (30) and (31) for the propensity functions are not unique, but it is particularly appealing because it reduces to MAK for ideal systems. In fact, consider an ideal latticemodel solution with N particles including both solvent and solute-see Appendix A for details. The free energy is where N = a n a (including solvent particles in the sum) and µ 0 a is the standard-state chemical potential of species a (taken with respect a given reference state noted as '0'), which may depend on parameters like temperature or nature of the solvent and the solute a. Then, the rates take the following form: =k ρ e β a r ρ where (N − a r ρ a )!/N ! can be approximated by N − a r ρ a . Setting we obtain that f +ρ coincides with the MAK propensity function (21), and similarly for f −ρ and Eq. (22).
We conclude this Section with a remark on the reaction constant, k ρ : In Eqs. (30) and (31) we have assumed that k ρ is a constant of the reaction, independent on the state n of the system. However, in general k ρ may depend on n, because the system itself is part of the environment where the chemical reactions take place. These effects can be disregarded for most cases in ideal solutions (since they are usually dilute), but they may not be negligible in non-ideal systems. For instance, in the case of phase separation, the multiple phases of the system may constitute very different environments for the chemical reactions, accelerating them or slowing them down.
Independently of whether k ρ in Eqs. (30) and (31) depends on the system state or not, detailed balance, Eq. (29), must still hold. This means that the forward reaction constant for a state n must be equal to the backward reaction constant for a state n + v ρ . One way to ensure this equality while keeping the state-dependency of the reaction constants, is to make F with the entropic term a log n a ! removed [37]. The quantities F 0 (n) and F 0 (n + v ρ ) denote the free energy of the system before and after the reaction, subscripts I and II specify the system phase, and ∆F the height of the free-energy barrier. For each phase, there are two minima in the free-energy landscape, corresponding to whether the reaction has occurred or not, see the left-and right-hand minimum, respectively. In this example, ∆F depends on the phase the reaction takes place in, and, thus, the reaction constant k ρ would also depend on the environment in which the reaction occurs.
k ρ a function of the state n deprived of the reactant complex, i.e., n − r ρ , for the forward case, and of n + v ρ − s ρ for the backward one: where Eqs. (35) and (36) satisfy Eq. (29) because see Eq.
(2). The dependency above of k ρ on the system state can be pictured as follows. In analogy with the classical transition-state theory, we can think of the microscopic mechanism of a reaction as a random walk in a free-energy landscape [36,37], see Fig. 3. Then, the value of the rate constant k ρ depends on the height of the free-energy barrier ∆F of the reaction. While the free energies of reactants and products (the stable local minima in the reaction landscape) have free energies F 0 I , F 0 II defined by F , this is not the case for the barrier height ∆F . The dependency of the height of the barrier-and thus of k ρ -on the system state is precisely the one discussed in Eqs. (35) and (36), and it may strongly affect the CRN dynamics. In summary, we are connecting the chemical reaction rates to the free energy of the system F , but also to ∆F which sets the value of the reaction constants k ρ .

Deterministic description
When the particle numbers n are large compared to the number of reactants and products, r and s, respectively, the free-energy differences which appear in the rates (30) and (31) can be rewritten as where in the first line we expanded F to first order in r, in the second line we used the definition of the chemical potential of species a: and f(c) is the free energy per unit volume in the deterministic notation. Therefore, the currents in a deterministic, non-ideal CRN can be written as which is an expression conceptually similar to that given by other approaches to construct thermodynamically consistent dynamics for deterministic CRNs [24,26]. Once again, the currents (40) match their ideal MAK counterpart (19) if the chemical potentials used in the rates are those of an ideal solution, i.e., µ a = 1/β log c a + µ 0 a . Here and in the rest of the text, dimensional arguments of the logarithms remain due to the fact that we are absorbing the effect of the total concentration in µ 0 , i.e., the original chemical potential was µ a = 1/β log(c a /c tot ) + µ 0 a , where c tot = a c a (the sum includes the solvent), but since variations in c tot can be neglected − log(c tot ) is just a constant and is absorbed into µ 0 a (and, thus, into k ρ ).
As in the previous section, if we assume the rate constant is state-dependent then the currents are given by wherek ρ is still a constant and any dependency of the rate constant on the state is given by the function g ρ (c).

B. Non-equilibrium systems
So far we considered the propensity functions of equilibrium CRNs. Given the large number of physically interesting systems which are out of equilibrium, such as living beings, in what follows we will generalise the analysis of Section III A to a specific type of nonequilibrium systems: those in which the work is done by the chemostats they are connected to.
Let us assume that N out of the N species in the system are connected to multiple particle reservoirs-chemostats: In the stochastic and deterministic description, each chemostat keeps constant the chemical potential of the species to which it is connected. Then, in general, the system will not relax to equilibrium, because of the work done on it by the chemostats. In the stochastic and deterministic description, the dimensions of the space of states or concentrations, respectively, is reduced to N − N ≤ N , since the chemostatted species are no longer dynamical variables.

Stochastic description
We assume that connecting the system to several chemostats does not alter any of the mechanisms of the chemical reactions, since it only tunes the concentration of the species to which they are connected, in order to match a given value of chemical potential. Then, reactions that involve both chemostatted and non-chemostatted species are driven in one direction by the work done by the chemostats inserting and removing particles from the system (in order to keep their chemical potentials constant). Given that the mechanism of reaction remains the same, in line with the previous section the rates of these driven chemical transitions are taken to be where now F is the free energy of the N − N non-chemostatted species, n contains the particle numbers of the non-chemostatted species only, and the summation over b runs over the chemostatted species. For the sake of clarity, in what follows we will reserve the index b for the chemostatted species, and the index a for the non-chemostatted ones. The rationale behind these relations is that the chemical reaction is still driven by free energy differences except that now the the free energy differences due to the consumption of chemostatted species is just given by the chemical potential of the chemostats µ b . The terms b r ρ b µ b and b s ρ b µ b in the exponential represents the chemical work done by the chemostats (with chemical potentials fixed at µ b ) when a reaction ρ occurs, which pushes the system out of equilibrium. The effect of the non-chemostatted species is still given by the free energy differences F (n) − F (n − r ρ ) and F (n) − F (n − s ρ ).
This implicitly assumes that the chemostatted species are abundant (so that the chemical potential does not fluctuate) and that they are ideal (negligible interactions with the non-chemostatted species). If the chemostatted species were not ideal, then the concentration of species might dynamically vary to match the chemostatted chemical potential as the particle numbers in the system change. Here, we only consider the simpler case of ideal chemostatted species and refer the interested reader to Ref. [24], where the case of non-ideal chemostatted species was analysed.
As in Section III A [see Eqs. (35) and (36)], the rate constants in Eqs. (42) and (43) may be generalised in such a way that k ρ depends on the system state: wherek ρ is independent of n. Propensity functions of this form have been suggested before in other contexts, such as in the modelling of molecular motors [38].

Deterministic description
Proceeding along the lines of Section III A 2, in the deterministic limit the above propensity functions result in the currents IV. STEADY-STATE DISTRIBUTION FOR COMPLEX-BALANCED,

NON-IDEAL CRNS
In what follows, we will prove one of the central results of this work, i.e., that the complex-balance condition allows us to generalise to non-ideal CRNs the result (23) [33] for the steady-state distribution of the network, which is generally unique (for details see Refs. [39,40]).
Namely, we claim that CRNs for which the complex-balance condition (15) holds, the steady state of the stochastic dynamics (10) with propensity functions (44) reads where the parametersμ a depend on the chemostats to which the system is connected and on the reaction constants of the network, but not on F . These parameters can be obtained from the CRN in the ideal and deterministic limit, thus significantly simplifying the task of obtaining analytically the steady-state of the system. Note that we reserve µ b for the chemical potentials of the chemostats whileμ a is a parameter that describes how the particle numbers at steady state of the non-chemostatted species depend on the non-equilibrium driving of the system. An additional necessary condition to prove this result is that the function g ρ must be the same for all reactions, i.e., g ρ = g; the relaxation of this hypothesis will be discussed in Section VI. The proof follows closely that of Anderson, Craciun and Kurtz [33], and here we only present its main steps-see Appendix B for a full proof. We will substitute the steadystate (46) into the dynamical equations, look for solutions where the probability flux across complexes vanishes, and obtain the complex-balance condition for a network with MAK, Eq. (20). We can thus conclue that, if the network modelled deterministically with MAK is complex-balanced at steady state, i.e. Eq. (20) is satisfied, then Eq. (46) is the steady-state probability distribution of its stochastic non-ideal counterpart. Furthermore, the parameters µ a in Eq. (46) can be obtained by solving Eq. (20).
By inserting the ansatz (46) in Eq. (10) with propensity functions of the form (44) and g ρ = g for all reactions, at steady state we obtain The previous equation is satisfied if, for each complex z, we have Given that in the previous equation the complex z is fixed, it can be simplified and yields Eq. (49) can be shown to be equivalent to the complex-balance condition for deterministic CRNs with MAK, Eq. (20), with rate constants given by Eqs. (51) and (52). These rate constants include the contribution of the standard-state chemical potentials µ 0 a and the chemostats, as is usually the case in MAK [41] (although, without loss of generality, for the purposes of this result, all µ 0 a can taken to be 0). Hence, a CRN for which the deterministic steady-state is complex balanced allows for a steady state of the form (46) for its stochastic and non-ideal version. Solving Eq. (20) for the steady-state concentrations with MAK and rate constants (51) and (52) yields c * a and, thus, the parametersμ a [via Eq. (50)] which appear in the steady-state distribution (46). The exponential relationship between the concentrations c * a andμ a reflects the logarithmic contribution of concentrations in the ideal chemical potential of solutes: µ a = µ 0 a + 1/β log c a . Equation (46) shows that the steady-state distribution of a non-ideal complex-balanced CRN has the form of an effective Boltzmann distribution, with the standard-state chemical potentials µ 0 a shifted byμ a (typically µ 0 a would be included within F ). From the physical standpoint, it is interesting to note that in Eq. (46) the free-energetic contribution F and the non-equilibrium term aμ a n a factor out.
This result is similar to Theorem 6.6 of Ref.
[33] -here Eq. (26)-but we have generalised it slightly to include rates of the form (44), which includes the function g that could be of interest in phase-separated systems as it modulates the rates depending on the environment. Moreover, our approach relates both the rates (44) and the steady-state distribution (46) to thermodynamic quantities, like free energies and chemical potentials.
In what follows, we will illustrate the result (46) with a minimal working example of a complex-balanced CRN, and compare its predictions with numerical simulations.

A. Example
Let us consider the following CRN-see Fig. 4 for a graphical representation: B C, with a free energy taken from a regular-solution theory (where each particle, including the solvent, occupies a finite volume and thus total volume is linked to the total number of

particles), see Appendix A for details.
For the sake of simplicity, we assume that the solvent particle number, n sol , is conserved, and allow the total volume to vary: where the total number of particles of species A, B and C, is kept constant in the CRN defined in (53).
Since the CRN (53) is complex balanced (which can be checked a posteriori), its steady state in the stochastic description and with propensity functions (44) can be obtained from its deterministic dynamical equations (11). To achieve this, we write the stoichiometry matrices which, together with the reaction constants given by Eqs. (51) and (52) and the free energy (32), completely define an ideal CRN. For simplicity, we assume that the standard-state chemical potentials µ 0 take the value 0 and thatk ρ = 1 for every reaction ρ. Finally, as an example, we take the non-equilibrium contribution of the chemostats to be present only in the reaction C A+D, with b r b µ b = 0 and b s b µ b = log(5/2)β −1 . These considerations, together with MAK [Eq. (19)] and the dynamics (11), yield the following set of deterministic and ideal equations for the CRN: Note that in the system derived from the matrices (56) d t n A = d t n D . For the sake of concreteness, we take as initial conditions and the solution of the above system at steady state is n * . According to Eq. (50), we have the following identity: n * a = e −βμa , which enables us to obtain the values ofμ a and the steady-state probability (46).
Note that there are two conservation laws, Eqs. (58) and (59), and four chemical species: hence, π neq (n) is a distribution with only two independent variables.
In order to evaluate Eq. (46) explicitly, let us assume that the system has the following regular-solution free energy where the first addend is an entropic term, the second corresponds to the internal energies of the chemical species taken with respect to that of species B, and the third to an interaction between species A and C. Setting χ = 10 and µ 0 A = µ 0 C = µ 0 D = log 2, we obtain the bimodal steady-state probability depicted in Fig. 5, which closely matches the one obtained from a simulation of the same CRN using the Gillespie algorithm [42]. Simulations were started in parallel from random Poissonian initial conditions satisfying the constraints above, and the samples were obtained after the simulation relaxed to steady state. Note that, in order to arrive to the set of Eqs. (57), we assumed all µ 0 a = 0, while in the free energy we are giving them a different value. It would have been equivalent to insert these values of µ 0 a into the system (57) and omit them in the free energy (60).

V. LYAPUNOV FUNCTION FOR COMPLEX-BALANCED STEADY STATES
A Lyapunov function is a function that is minimised by the dynamics of the system and takes the value 0 at steady state. Under fairly general conditions, the logarithm of the steady-state probability distribution in the stochastic CRN is a Lyapunov function of the deterministic one [43,44].
While the exact form of the Lyapunov function has been obtained for ideal and complexbalanced CRNs [30,45], here we demonstrate that for non-ideal, complex-balanced CRNs the following function decreases with the dynamics where the factor 1/V has been inserted to maintain the magnitude intensive while V → ∞. Our approach generalises the results of Anderson and Nguyen [46] for product-form stationary states of CRNs. We will call the function (61) a Lyapunov function: This is a slight abuse of terminology, because we will only prove that L decreases with the dynamics, not that its value is zero at steady state. L only takes the zero value if where the asterisk denotes values at steady state. Given that Z is a normalisation factor for the stochastic complex-balanced CRN at steady state, see Eq. (46), it reads which, for a large (many particles) deterministic CRNs, can be evaluated using the saddlepoint approximation, where the sum is evaluated at the minimum value of the argument of the exponential. If the deterministic system is monostable, then the argument of the exponential has a single local minimum. Therefore, for monostable CRNs, this approximation will yield the correct value and the Lyapunov function Eq. (61) will take the value 0 at steady state. However, care must be taken when handling multistable CRNs in this way, which is why, in order to avoid this complexities, we will not prove that Eq. (61) takes the value 0 at steady state in general. Nevertheless, the fact that this function decreases with the dynamics is sufficient for our purposes.
In what follows we sketch the proof that, for complex-balanced non-ideal CRNs, the Lyapunov function (61) is a decreasing function of time-for a full step-by-step proof, see Appendix C.
Given that the normalisation factor Z does not depend on time but only on the nonequilibrium steady state, the time derivative of L is where in the second line we used Eqs. (11) and (45), together with the assumption g ρ = g. After adding and subtracting terms of the form a r ρ aμ a in the exponentials (of the form a s ρ aμ a for the second exponential), we repeatedly apply the inequality e x (y − x) ≤ e y − e x to the sums of chemical potentials, and obtain The expression in the right-hand side (RHS) above can be split in terms of the different complexes in the system: For a complex-balanced system, it can be shown that the expression in curly brackets in Eq. (66) vanishes for each complex z independently, as a consequence of the complex-balance condition for MAK systems, Eq. (20). Then and L decreases, or remains unchanged, along a trajectory. We conclude that, unlike in classical equilibrium systems, here it is not the F that is minimised by the dynamics, but a free energy (61) where the standard chemical potentials µ 0 a are shifted byμ a . This shift, which is entirely due to the non-equilibrium contribution of the chemostats, enables the system to present non-vanishing chemical-reaction net flows between species at steady-state, which, in the following, we will call chemical currents.

A. Spatially heterogeneous systems
In order to describe phase-separating systems, in what follows we will incorporate in our framework spatial inhomogeneities. In the deterministic description, concentrations are now a function of space, c a (x) within a volume Ω, and the free energy is a functional of these concentrations, F [c].
The time derivative of the concentrations is given by the following reaction-diffusion (RD) equation where the time dependence of c is omitted, the first term in the RHS of the equation represents diffusion, and the second one the chemical reactions. As in the linear irreversible thermodynamics framework [47], the driving force of the diffusion current J a is the gradient of chemical potentials, ∇µ a : the diffusive currents then read where M ak is the mobility matrix. We assume no-flux boundary conditions for the non-chemostatted species, where ∂Ω denotes the boundaries of the volume Ω.
We now consider a generalisation of the Lyapunov function (61) to inhomogeneous systems. In the following, we will show that the dynamics (68) for complex-balanced networks minimise the Lyapunov functional where F [c] is the free energy of the system, which depends on the concentration profile through The time derivative of the Lyapunov functional (71) yields where µ a (x) = δF [c(x)]/δc a (x) is the local chemical potential. By applying the results of Section V, Eq. (67), at every spatial point x, we obtain that the second term in the square brackets of the RHS of Eq. (73) is negative or zero. Therefore, to prove that dL/dt ≤ 0 it is sufficient to show that the first term in the square brackets of the RHS is also negative or zero. In this regard, we note that The first term in the RHS of the last equality vanishes due to the divergence theorem and Neumann boundary conditions (70). By observing thatμ does not depend on space, the addend containingμ in the second term vanishes (∇μ a = 0). Finally, if the Onsager reciprocal relations for the mobility matrix M ak hold [48], then the addend containing µ a (x) in the second term is necessarily positive, because it represents the entropy production of a diffusion process [37,47]. The Onsager reciprocity relations ensure that a system relaxes to equilibrium in the absence of external work. Thus, the condition that the Onsager relations hold is not a limitation of the result but a consequence of thermodynamical consistency. Combining the results above, we obtain that It follows that, for a non-ideal, complex-balanced system, L decreases, which we can now use to obtain useful information about the steady state, along the lines of the free energy minimisation for systems at thermodynamic equilibrium. Therefore, for a complex-balanced system, we can minimise L (subject to constraints in particle numbers) in order to obtain the concentration profile at steady state. This minimisation results in constraints for the steady-state profile of the form where λ is a Lagrange multiplier that enforces the particle-conservation constraint-for further details see Section V B. Equation (76) implies that at steady-state in a complexbalanced solution there cannot be any diffusive currents, since the chemical potential is constant throughout space and the force driving diffusion currents is ∇µ a (x). Nevertheless, chemical currents can exist at steady state, as noted in the previous section, and the concentration profile may not be homogeneous. This is a major consequence of the present work.

Example
Let us consider the following CRN (see Fig. 4 for a graphical representation): with a free energy taken from a regular solution theory, as before. We will first assume that the system is homogeneous and later we will analyse the full reaction-diffusion system. For simplicity, we assume that the system is driven out of equilibrium solely by imposing a non-equilibrium chemical potential difference in the transition from C We take all the reaction constants k equal to each other, and note that the network is necessarily complex-balanced, as all chemical reactions are unimolecular (in unimolecular networks each of the species is a complex, hence the steady-state condition is equivalent to the complex-balance condition, if g ρ = g for every reaction ρ).
Proceeding along the lines of Section IV A, in the stochastic description the steady-state of the CRN (77) with propensity functions (44) can be obtained from the following ideal and deterministic rate equations: The solution of the above system at steady-state is We obtain the values ofμ a by identifying n a with e −βμa [as in Eq. (50)]. Noting that we can express such potentials with respect to that of species C, we obtain the Lyapunov function of the system where a detailed expression of the normalisation constant Z is not essential here, because Z is constant along the dynamics, and it does not alter the location of the minima of L in the space of concentrations c.
Unlike above, we will now describe the amount of species with reference to the fraction of volume they occupy at each point of space φ(x). Then, where φ solv (x) is the volume fraction of the solvent and the sum runs over solutes only. Equation (79) states that the solution is incompressible and, thus, φ solv (x) = 1 − a φ a (x). The reason for using volume fractions instead of concentrations is threefold: It is the convention normally used in phase separation studies and regular solution models, it enforces incompressibility (which is the case in most liquids) and is dimensionless. For simplicity we will assume that the molecular volumes of every species is the same, so that φ a (x) is proportional to c a (x). Hence, the following regular solution free-energy density can describe spatial inhomogeneities in an incompressible solution: where the first two terms in the RHS are entropic terms and the following two represent the interactions between the solutes. The last addend represents the free-energetic cost of spatial inhomogeneities in the concentration profiles, and is known as Cahn-Hilliard term [49].
Assuming the system is one-dimensional, the resulting Lyapunov functional for the RD system is whereμ a are the ones obtained for the homogeneous system and do not depend on the coordinate x. We set χ AA = −2, χ AB = −7, ∆ = −2, κ A,A = κ B,B = 5, κ A,B = 1 and any other Cahn-Hilliard coefficient equal to 0. With this parameter set, the reaction-diffusion system exhibits phase separation at steady state (see Fig. 6). By entering this free energy in the RD equations (68) and assuming no state dependency of the reaction constants k ρ , we obtain a set of equations which describes the dynamics of the system. Figure 6 shows that the Lyapunov functional (81) is minimised by the dynamics, and that the non-equilibrium steady state is characterised by phase coexistence.
Finally, in Fig 6 C, the net reaction flux J +ρ − J −ρ at steady state as a function of the spatial coordinate is depicted. This net reaction flux is constant in space and, given the topology of the CRN (77), is equal for all reactions ρ. The fact that the net reaction flux is independent of the spatial coordinate despite the varying concentrations (see Fig. 6 B) is a result of chemical reaction fluxes being driven by the chemical potential, which, as argued above, is constant-see Eq. (76). Note that this is also a consequence of having dropped the dependency of the reaction constants on the environment via a function g(c). If all reaction constants were subject to this modulation (which has to be the same for every reaction for our results to hold), then the reaction rates at steady state could be space-dependent but the chemical potential would still be constant.

B. Phase Diagram of a chemically reactive mixture
Since the Lyapunov functional for complex-balanced systems discussed in Section V A is minimised by the dynamics, it carries plenty of information on the steady state.
Along the lines of phase separation for equilibrium systems, the steady state can be obtained by minimising L[c] with respect to c, subject to certain constraints, e.g., particle conservation. The concentration profiles which realise the absolute minimum of L[c] may be either spatially uniform, or depend on space, according to the system parameters. On a qualitative level, the phenomenology of a complex-balanced system does not change much with respect to that of a non-ideal solution at equilibrium, but the non-equilibrium terms may alter the phase diagram, thus tweaking the onset of phase separation.
To illustrate this point, in this Section we consider a non-ideal solution with the CRN (77) in the deterministic description, and obtain its phase diagram. Therefore, we minimise the Lyapunov functional (81) of Section V A 1, with the particle-conservation constraint where φ a (x) is the volume fraction of species a, and the constant φ N fixes the total volume fraction of the solutes. Then, the function that needs to be minimised is the Lagrangian where λ is the Lagrange multiplier associated with the conservation of solutes, and |Ω| the volume of Ω. For the sake of simplicity, we take the typical lengthscale of Ω to be large with respect to inter-species interfaces: as a result, the volume fractions φ(x) can be approximated by piecewise constant functions. If the system phase separates, we assume that only two homogeneous, distinct phases, which we denote by '1' and '2', will appear. Within this assumption, the Lagrangian (83) reads where Ω 1 and Ω 2 stand for the volumes phases 1 and 2, respectively, with |Ω 1 | + |Ω 2 | = |Ω|, and we consider the following free-energy density where all species except A are considered non-interacting and A interacts with itself (χ > 0 implies an effective attraction between A particles). The minimisation of Λ yields the phase diagram in Fig. 7, see Appendix D for details. Phase separation occurs in regions II and III of the phase diagram, as shown in the concentration profiles displayed in the insets.
From the form of the free energy (85), we can see that is the species A that drives phase separation, since for χ > 0 the free energy will favour segregating A from the rest of the solution. Thus, whether the steady-state displays one phase or a coexistence of phases also depends on the value of the non-equilibrium chemical potential difference ∆, which can alter the concentration of A at steady state and, hence, modulate phase separation, as can be seen in Fig 7 B.

VI. DISCUSSION
In this work, we have shown that for a chemically reactive non-ideal solution we can obtain results for complex-balanced networks analogous to those for ideal solutions, provided that the system is modelled in a thermodynamically consistent way. This implies that the rates of the chemical reactions incorporate the interactions between the species in the system and, therefore, mass-action kinetics (MAK) no longer holds. By generalising MAK to a non-ideal solution, we obtained the steady-state probability distribution for a stochastic complex-balanced CRN and the Lyapunov function of its deterministic counterpart, which determines the phase diagram of the system.
Our results are of particular importance for non-equilibrium phase-separating systems. By combining previous results from the mathematical theory of CRNs [33,46] and concepts of non-equilibrium thermodynamics [37,48], we found that the resulting complex-balanced RD system cannot sustain diffusion currents at steady state, see Eq. (76). Since, in many cases, diffusion currents are required for pattern formation in reaction-diffusion systems, breaking complex balance is a necessary condition to obtain such patterned steady states, at least when interactions are modelled in a thermodynamically consistent way unlike, e.g.,  those in Refs. [21,22]. In this regard, complex balance can be broken in two ways: First, by choosing a suitable network topology that allows for a steady state which is not complexbalanced, as in Ref. [23]. Second, in a system where different phases coexist, by allowing the reaction rates to depend differently on local environment: For example, in Ref. [27] a patterned steady-state is produced by allowing one (and only one) of the reaction constants to depend on the concentration of an enzyme which localises in one of the phases. Mathematically, this violates one of the necessary conditions for our results to hold, namely g ρ = g (see Section IV), thus allowing for more general steady states.
In biological cells, phase separation has been hypothesised to perform many functions, such as, accelerating biochemical reactions within the condensate irrespective of the rate of the reaction in the dilute phase [18,50]. The present work implies that, in order to control chemical reactions in each of the phases independently (at steady state) breaking complex balance is necessary, by virtue of Eq. (76). Indeed, in a complex-balanced system, the chemical potential of every species is constant throughout space. Then, given that the force driving the chemical reactions are the chemical potentials, the reaction rates in both phases are related, making it impossible to regulate the rates of chemical reactions in each phase in a fully independent way, and suggesting that breaking complex balance in one of the two ways outlined above is crucial for such control.
Overall, complex balance is known to be a key feature of CRNs which determines not only their behavior [32,33] but also their thermodynamic properties [34]. In this analysis, we further stress the connection between the characteristics of the reaction network and the thermodynamically consistent structure of the physical system, in an effort to generalise results from ideal CRNs, and explore non-equilibrium dynamics of complex-balanced networks. However, little is known about non-complex balanced systems and, given our results, further research regarding the behaviour of this type of networks out of thermodynamic equilibrium would be of the utmost importance, both from the physical [51] and biological [50] point of view.
• Authors' contributions: A.M.M. conceived the study and wrote the paper. M.C.
contributed to the discussions and revisions.

A. FORM OF THE PROPENSITY FUNCTIONS FOR A REGULAR-SOLUTION THEORY
We now consider a model of a solution based on a lattice where each chemical species (including the solvent) occupies one lattice site, thus neglecting differences in molecular volumes.
In a lattice with N sites (note that the number of sites is proportional to the volume) occupied by N different species, with N a=1 n a = N , the configurational entropy is given by where the argument of the logarithm is the number of microstates. The internal free energy of each species a is given by the standard-state chemical potential µ 0 a . We incorporate in this regular-solution model interactions among neighbouring sites, whose energy (in the mean-field approximation) reads a,k [n a n k (1 − δ a,k ) + n a (n a − 1)δ a,k ] = a,k χ ak 2N (n a n k − n a δ a,k ) , where χ ak represents the interaction energy between species a and k, and it can also be interpreted as the matrix of virial coefficients. Taken into account the previous considerations, the free energy for a homogeneous mixture of chemical species in the regular-solution model reads where the first two terms in the last line represent the ideal free energy (see Eq. (32) in the main text), while the last term is exclusively due to interactions between solutes. With this expression of the free energy we can now derive an expression for the propensity functions (30) and (31). The forward (or backward) rates are a function of the free energy difference of the complex: where ∆F id is the ideal part of the free energy difference of the complex, given by Eq. (33) in the main text. In the RHS of Eq. (89), only the first two terms are non-vanishing as we approach the thermodynamic limit (N → ∞, while keeping n a /N fixed): hence, for large systems, the rest of the interacting terms are negligible. However, for a unimolecular reaction, since the free energy difference takes a particularly simple form, we have that where a is the reactant of the reaction +ρ.
In the thermodynamic limit, Eqs. (30), (31) and (89) imply that the deterministic rates can be written as where the particle numbers have been replaced with concentrations (an additional logarithmic factor has been absorbed into the rate constant k +ρ , as explained in the main text, Section III A 2) and the part of the chemical potential representing the internal energy has also been absorbed in the rate constant k +ρ . Setting the rates (91) match the general expression given in the main text, Eq. (40).

B. PROOF OF THE COMPLEX-BALANCED DISTRIBUTION
In this Section we present the full proof of the result (46). At steady state, the CME with rates of the form (44) and g ρ = g for all reactions ρ, reads By dividing the previous expression by P (n) and substituting the ansatz (46) into it, we obtain ρk ρ g(n − s ρ )e β[F (n)−F (n−s ρ )+ a v ρ We now rewrite the relation above in terms of a summation over each of the complexes z ∈ C separately z ρ|s ρ =zk where the subscript 'ρ|s ρ = z' denotes that the sum runs only over reactions ρ whose product complex s ρ is equal to z. This previous equation will be satisfied if for every complex z. For any given complex z, Eq. (96) can be rewritten in the following form: We now divide both sides by g(n − z) exp{β[F (n) − F (n − z)]}, and obtain where we have substituted v ρ = s ρ − r ρ and, depending on the reactions over which the sum runs, one of this terms can be replaced by the complex z. Finally, given that in Eq. (98) z is fixed, we can divide both sides by exp(β a z aμa ), yielding , where in the third line we have used Eq. (11) with currents given by Eq. (45). By adding and subtracting terms of the form a r ρ aμ a in the exponentials (of the form a s ρ aμ a for the second exponential), we rewrite the previous equality as dL dt = ρ k k ρ g(c)(µ k +μ k )(s ρ k − r ρ k )e β[ a r ρ a (µa+μa)− a r ρ aμa We now consider the inequality e s (t−s) ≤ e t −e s -which results from 1+x ≤ e x , ∀ x ∈ R with x = t − s-and apply it to the sums of chemical potentials. In the first term in the RHS of Eq. (101), we set s = a (µ a +μ a )r ρ a and t = a (µ a +μ a )s ρ a , and conversely in the second term. We then obtain dL dt ≤ This expression can now be separated in terms of the different complexes in the system: In order to find the steady state of the system, we need to minimise the Lyapunov functional or, in the presence of particle-conservation constraints, the Lagrangian, (83). A substantial simplification can be made by neglecting the contribution of the interfaces, i.e., considering the system as composed of two homogeneous phases. In this approximation, the function which needs to be minimised is the Lagrangian (84), which depends on eight independent variables: φ p a for a = A, B, C, p = 1, 2, λ and |Ω 1 |.
First, we reduce the dimensionality of the problem by equating the derivatives of the Lagrangian with respect to the concentrations of the species: where a and k denote two chemical species, and p = 1, 2 refers to the phases. Equation (106) for a system at equilibrium yields the equality of chemical potentials (with their appropriate stoichiometry). Here, however, Equation (106) includes the shifted chemical potential term µ due to the out-of-equilibrium complex-balancing. For a simple free energy like Eq. (85), Eq. (106) implies which reduces the problem to just four variables: φ 1 A , φ 2 A , λ and |Ω 1 |. Finally, given that we are interested in the phase diagram of the mixture and not in the actual steady state of the solution (i.e., we do not need to know how much volume each of the phases occupies), we can avoid solving for |Ω 1 |. This can be achieved by enforcing the stationarity condition of the Lagrangian with respect to the volume: which, together with yields a fully determined system for the unknowns φ 1 A , φ 2 A and λ (the dependency on |Ω 1 | drops out). The resulting equations for such unknowns are transcendental equations which, in general, have no explicit analytical solution. Therefore, they need to be solved numerically. Even the numerical solution is involved as the parameters near criticality, which is why in Fig. 7 the density of data around the critical point decreases.