Thermodynamic and Stoichiometric Laws Ruling the Fates of Growing Systems

We delve into growing open chemical reaction systems (CRSs) characterized by autocatalytic reactions within a variable volume, which changes in response to these reactions. Understanding the thermodynamics of such systems is crucial for comprehending biological cells and constructing protocells, as it sheds light on the physical conditions necessary for their self-replication. Building on our recent work, where we developed a thermodynamic theory for growing CRSs featuring basic autocatalytic motifs with regular stoichiometric matrices, we now expand this theory to include scenarios where the stoichiometric matrix has a nontrivial left kernel space. This extension introduces conservation laws, which limit the variations in chemical species due to reactions, thereby confining the system's possible states to those compatible with its initial conditions. By considering both thermodynamic and stoichiometric constraints, we clarify the environmental and initial conditions that dictate the CRSs' fate-whether they grow, shrink, or reach equilibrium. We also find that the conserved quantities significantly influence the equilibrium state achieved by a growing CRS. These results are derived independently of specific thermodynamic potentials or reaction kinetics, therefore underscoring the fundamental impact of conservation laws on the growth of the system.


I. INTRODUCTION
Chemical thermodynamics provides solid physical principles for explaining the energetics and for predicting the fate of chemical reaction processes [1,2].Applications of these principles to autocatalytic reactions are essential to elucidate physical constraints on the capability of self-replication [3,4].In the self-replication process, both the chemical components and the encapsulating volume of the system have to grow in a coherent manner [5][6][7][8].This consideration introduces a unique theoretical challenge to establish the thermodynamic consistency between autocatalytic reactions and volume expansion because the conventional chemical thermodynamics is based solely on the density (concentration) of chemicals, presuming a constant volume [9][10][11][12][13][14][15][16][17][18].
We have recently established a thermodynamic theory for growing systems, in which the number of chemicals and the volume are treated based on the rigorous thermodynamic basis [19].Accordingly, the theory generally formulates physical conditions to realize the growth of the system, identifies several thermodynamic constraints for the possible states of the growing system, and derives the form of entropy production and heat dissipation accompanying growth.However, this theory can only address systems with regular (full-rank) stoichiometric matrices, thus limiting its applicability to a set of minimum autocatalytic motifs [20].Given that chemical reaction systems (CRSs) are subject to stoichiometric constraints in general, and various biological functions are robustly realized by specific stoichiometric properties [21,22], it becomes imperative to delve into the influence of stoichiometry to growing systems for comprehensive understanding of the thermodynamics of self-replication.
This study elucidates the complex interplay between thermodynamics and stoichiometric conservation laws for the growing system and its consequence to the fate of the growing systems.By disentangling the geometric relationship among chemical numbers, densities, and potentials, we establish the conditions for the system to grow, shrink, or equilibrate, while simultaneously satisfying the second law of thermodynamics and the conservation laws.Furthermore, we show that the existence of conservation laws qualitatively alter the fate of the system and the geometric properties of the equilibrium state.
This paper is organized as follows.We devote Sec.II to outline our main results accompanying with an illustrative example of a growing CRS.Here, we clarify the environmental and initial conditions that determine the fate of the system: growth, shrinking or equilibration.In Sec.III, we explain the geometric structure of the growing system to characterize the equilibrium state as a maximum of the total entropy function, which the system may achieve.In Sec.IV, the form of the total entropy function is further investigated to determine if the system grows or shrinks by following the second law of thermodynamics.In Sec.V, we provide a mathematical proof of our main results outlined in Sec.II.Finally, we summarize our work with further discussions in Sec.VI.For better readability, we list the symbols and notations in Appendix L.

II. OUTLINE OF THE MAIN RESULTS
We outline our main results before presenting their derivations.In Sec.II A, we give a thermodynamic setup of a growing chemical reaction system (CRS).By introducing the chemical potential space, we analyze candidates of chemical equilibrium states in Sec.II B. By computing accessible regions of the system, we characterize the equilibrium state as an intersection of the accessible region and the candidates in Sec.II C. In Sec.II D, we present our main claims with the above preparation.We comment our claims for a special class of the growing CRS in Sec.II E. These results are summarized in Table I.Then, we demonstrate the claims by numerical simulations using a simple example of the growing CRS in Sec.II F. Finally, in Sec.II G, we outline the derivations of our results, which will be presented in the subsequent sections.

A. Thermodynamic setup
We consider the following thermodynamic setup in the present paper (Fig. 1).A growing open chemical reaction system is surrounded by an environment.We assume that the system is always in a well-mixed state (a local equilibrium state), and therefore, we can completely describe it by extensive variables (E, Ω, N, X).Here, E and Ω represent the internal energy and the volume of the system.N = {N m } ∈ R NN >0 denotes the number of chemicals that can move across the membrane between the system and the environment, which we call open chemicals hereafter.X = {X i } ∈ X = R NX >0 is the number of chemicals confined within the system; the indices m and i respectively run from m = 1 to N N and from i = 1 to N X , where N N and N X are the number of species of the open and the confined chemicals.The environment is characterized by intensive variables ( T , Π, μ), where T and Π are the temperature and the pressure; μ = {μ m } ∈ R NN is the chemical potential correspond-ing to the open chemicals.Also, ( Ẽ, Ω, Ñ ) denote the corresponding extensive variables.
In thermodynamics, the entropy function is defined on (E, Ω, N, X) as a concave and smooth function Σ[E, Ω, N, X].We also write the entropy function for the environment as Σ T , Π,μ [ Ẽ, Ω, Ñ ], and the total entropy can be expressed as where we use the additivity of the entropy.Further, the entropy function for the system has the homogeneity.Therefore, without loss of generality, we can write it as where σ[ǫ, n, x] is the entropy density and (ǫ, n, x) := (E/Ω, N/Ω, X/Ω).In this work, we consider only a situation without phase transition, and therefore, we assume that σ[ǫ, n, x] is strictly concave.The dynamics of the number of confined chemicals X(t) is described by where J(t) = {J r (t)} represents the chemical reaction flux; S = {S i r } denotes the stoichiometric matrix for the confined chemicals (see Fig. 1).The index r runs from r = 1 to N R , where N R is the number of reactions.Also, Einstein's summation convention has been employed in Eq. (3) for notational simplicity.By integrating Eq. (3), we obtain where X 0 = {X i 0 } denotes the initial condition and Ξ(t) = {Ξ r (t)} is the integration of J(t), i.e., the extent of reaction.
In this work, we assume that the stoichiometric matrix S does not have the right kernel space, dim Ker[S] = 0. ( By contrast, the left kernel space exists and generates a set of conserved quantities: i.e., the vector Here, U is a basis matrix: {U l i } := (U 1 , U 2 , ...) T whose row vectors U l form the basis of Ker[S T ], i.e., U S = 0, and l runs from 1 to dim Ker[S T ].Since we are mainly interested in the growth of the system, the time evolution of X(t) should perpetually increase the number of all confined chemicals while satisfying the conservation laws (Eq.( 6)).In order to realize this situation, the stoichiometric matrix S should satisfy which is known as the condition for the productive S [20].The productive S excludes the mass conservationtype laws, the existence of which trivially precludes the number of confined chemicals to increase continuously.
The number of open chemicals changes through the reactions and the diffusions.By defining the stoichiometric matrix for open chemicals as O = {O m r } and the diffusion fluxes as J D (t) = {J m D (t)}, we have Example 1: To give an illustrative example, we con-sider the following chemical reactions: Since N R = 2 in this case, we can represent the number of confined chemicals by introducing Ξ (11) For this example, dim Ker[S T ] is one, and U = (−1, 0, 1) satisfies U S = 0. Thus, the reactions have the conserved quantity L = U X = −X 1 + X 3 .We can intuitively see the conserved quantity from the chemical equations in Eq. ( 9) as follows.The numbers of A 1 and A 3 , i.e., X 1 and X 3 , change in the same manner: each chemical is consumed by one molecule in the forward reaction of R 1 and is produced by one in that of R 2 .Thus, the difference between X 1 and X 3 is conserved when each of the reactions occurs.
In this paper, we assume that the time scale of the chemical reactions is much slower than that of the others, i.e., J ≪ J E , J Ω , J D .This means that we consider the isothermal and isobaric process with a balance of chemical potentials for open chemicals.This assumption allows us to separately analyze the slow dynamics of chemical reactions from the fast dynamics of the energy, the volume, and the chemical diffusion.As shown in Appendix A, our dynamics is effectively governed by Eq. ( 3), and the other variables (E, Ω, N ) are rapidly relaxed to the values at the equilibrium state of the fast dynamics.As a result, (E, Ω, N ) = (E(X), Ω(X), N (X)) is obtained as a function of X.
With the above assumptions, we obtain the following expression for the total entropy (see Eq. (A11) in Appendix A), where X = X 0 + SΞ.Here, ϕ(•) is the partial grand potential density, which is given by a Legendre transformation of σ [ǫ, n, x]: Note that ϕ(x) is strictly convex because σ is strictly concave.The volume Ω(X) is variationally determined under isobaric conditions as Example 2: If we assume that the system and the environment are composed of ideal gas or solution, the partial grand potential density ϕ(x) is expressed as where R is the gas constant, and ν o ( T ) = {ν o i ( T )} and µ o ( T ) = {µ o m ( T )} are the standard chemical potentials of the confined and the open chemicals, respectively [78].
Since we assume that the environment also consists of the ideal gas, the chemical potential μ in Eq. ( 15) is represented as where ñ = {ñ m } is the density of the open chemicals in the environment.Furthermore, Eq. ( 14) gives the following expression for the volume Ω, which corresponds to the equation of state.
From Eq. ( 14), the density x ∈ X = R NX >0 is obtained by a nonlinear function ρ X (X) of X as Since Ω(X) is homogeneous, we have ρ X (αX) = ρ X (X) for α > 0. This fact implies that ρ X induces one-to-one correspondence between rays in the number space X and points in the density space X .Here, we write the corresponding ray to a density x, i.e., the fiber of x for ρ X , as

B. Candidates of chemical equilibrium states
We define the full grand potential density ϕ * (y) = ϕ * [ T , μ; y] by the Legendre transformation of ϕ(x) in Eq. ( 13): We mention that ϕ * (y) is strictly convex.Because of the one-to-one correspondence of the Legendre transformation by ϕ(x) and ϕ * (y), a state of the system is equivalently specified either by the density x ∈ X = R NX >0 of the confined chemicals or by its Legendre dual variable y = ∂ϕ(x) ∈ Y = R NX [79].The thermodynamic interpretation of y is the corresponding chemical potential to x.In addition, ϕ * (y) can be interpreted as the pressure of the system at the state y.Recalling the one-to-one correspondence by ρ X between rays in X and points in X , we notice that there exists the one-to-one correspondence between rays in X and points in Y. Here, we write the corresponding ray to a chemical potential y as The chemical equilibrium states are given by the balance of chemical potentials between reactants and products [78]: Thus, we can obtain the candidates of chemical equilibrium states by the solutions to Eq. ( 22): which we term the equilibrium manifold.Here, we note that M Y EQ (μ) = ∅, because we have assumed dim Ker[S] = 0 in Eq. ( 5), i.e., dim Im[S T ] = N R .

C. An equilibrium state as the intersection
The reaction dynamics in Eq. (3) with its conserved quantities L restricts the accessible region in X as This subspace is known as the stoichiometric compatibility class [80].By using the mappings ρ X in Eq. ( 18) and ∂ϕ, the accessible regions in the density space X and the chemical potential space Y are respectively obtained as Since the system evolves only on M Y STO (L) ∈ Y, the equilibrium state to which the system converges is represented by the intersection: if it is not empty.If the intersection is empty, the system never converges to the equilibrium state and the growth of the volume may happen.

D. Main claims of the present paper
In this subsection, we state our two main claims of this work before presenting their derivations in the following sections.
Our first claim provides the conditions to determine the fate of the system, i.e., growth, shrinking or equilibration.To formulate our claim, we define y min as Note that the productivity of S in Eq. ( 7) guarantees the existence and the uniqueness of y min (see the latter half of Appendix F).This y min gives the chemical equilibrium state whose pressure ϕ * (y min ) is the minimum in the candidates, M Y EQ (μ).
Claim 1 The fate of the system is classified by y min and the conserved quantities L in Eq. ( 6) as follows (see also Table I): 1.If ϕ * y min − Π > 0, the system grows and finally diverges for any L.
Here, L = 0 represents the case in which L l = 0 for all the components of L = {L l } in Eq. ( 6), and L = 0 indicates that at least one of the components is not zero.
We first notice that one of the conditions in the claim is given by the sign of ϕ * (y min ) − Π.Here, ϕ * y min represents the minimum pressure of the system within M Y EQ (μ) (see Eq. ( 28)), whereas Π is the pressure of the environment.The growth of the system occurs independently of the value of L when ϕ * (y min ) > Π, because the environmental pressure can not prevent the system from growing by counteracting any internal pressures ϕ * (y) for y ∈ M Y EQ (μ).The equilibration becomes possible only when ϕ * (y min ) − Π ≤ 0 is fulfilled.For the singular case of L = 0, the equilibration occurs only when the minimum internal pressure perfectly balances with the external one, namely, ϕ * (y min ) − Π = 0.When ϕ * (y min ) − Π < 0, the system shrinks and X(t) approaches 0 as t → ∞ because the external pressure dominates the internal one at the equilibrium state.For the generic L = 0, the fate of the system is qualitatively altered.The system equilibrates even if ϕ * (y min ) − Π < 0 owing to the nonsingular conservation laws, which prevent X(t) from approaching 0, i.e., shrinking.In fact, the equilibrium state y ∈ M Y EQ (μ) can have a higher pressure in this case than TABLE I.The fate of the system is classified by the sign of ϕ * (y min ) − Π and the conserved quantities L = {L l }.Equilibration is indicated by EQ(D) and EQ(I).For EQ(D), the equilibrium state depends on the initial conditions and the functional form of the reaction flux J(t).For EQ(I), it is independent for a given L. Here, y min = {y min i } is defined in Eq. ( 28) and it is identical with y EQ i = −μmO m r (S −1 ) r i for regular S.
ϕ * (y min ), and it balances with the external pressure Π.
From claim 1, the system converges to an equilibrium state when (1) ϕ * (y min ) − Π = 0 and L = 0, and (2) ϕ * (y min ) − Π < 0 and L = 0.For these cases, our second claim describes how the equilibrium states appear in the number space X.
Claim 2 In the number space X, the equilibrium states appear as follows (see also Table I): 1.When ϕ * (y min ) − Π = 0 and L = 0, the equilibrium state to which the system converges depends on the initial condition X 0 ∈ M X STO (L = 0) and the functional form of the reaction flux J(t).Such an equilibrium state lies on the ray r Y (y min ) in Eq. (21).
2. When ϕ * (y min ) − Π < 0 and L = 0, the equilibrium state is uniquely determined, irrespective of the initial condition X 0 ∈ M X STO (L = 0) and the form of the reaction flux J(t).

E. Regular stoichiometric matrix S
The fates of the system for the singular case of L = 0 are basically the same as the case that S is regular; no conservation laws exist so that dim Ker[S T ] = 0, and dim Ker[S] = 0 from Eq. ( 5) [81] (see Sec. V C for the reason why regular S can be interpreted as L = 0).Thus, the two claims clarify that the existence of conservation laws is a fundamental determinant of the fate of the growing system.We investigated the case of regular S in our previous work [19] (see also [82]).For comparison, we summarize the results here (see the cases for L = 0 and regular S in Table I).
For a regular S, the candidates of chemical equilibrium states, M Y EQ (μ) in Eq. ( 23), consist of precisely one point Here, we note that y min is identical with y EQ by Eq. ( 28).Then, we obtain Corollary 1 When the stoichiometric matrix S is regular, the fate of the system is classified by y EQ i = −μ m O m r (S −1 ) r i as follows (see also Table I): 1.If ϕ * y EQ − Π > 0, the system grows and finally diverges.
2. If ϕ * (y EQ ) − Π = 0, equilibrium states form the ray r Y (y EQ ) in Eq. (21).The system converges to one of them, depending on the initial condition X 0 ∈ X and the functional form of the reaction flux J(t).
This statement corresponds to claim 1 in [19].

F. Numerical demonstrations
In this subsection, we demonstrate our claims by using the specific reactions in Example 1 with the ideal gas case (see Example 2 ).
We first calculate the equilibrium manifold M Y EQ (μ) in Eq. ( 23) to obtain y min in Eq. ( 28).For the specific reactions in Eq. ( 9) and the stoichiometric matrices in Eq. ( 10), the simultaneous equations in Eq. ( 22) are written as The solutions y = {y i } are expressed as where h ∈ R is the coordinate of Ker[S T ].The set of solutions in Eq. ( 30) represents the equilibrium manifold M Y EQ (μ).Second, using Eqs ( 15) and ( 20), the full grand potential density ϕ * (y) is obtained for the ideal gas or solution as Using the solutions in Eq. ( 30), ϕ * (y) is represented on In Fig. 2(a), we numerically plot the pressure φ * (h) on M Y EQ (μ).Here, we have a unique h min = (−ν o 1 + ν o 3 − μ1 + 2μ 2 )/2 that attains the minimum φ * (h).Substituting it into Eq.( 30), we obtain y min and the minimum pressure ϕ * (y min ).With the parameters given in the caption, ϕ * (y min ) = 14.25.
To demonstrate that the fate of the system is classified by y min and the conserved quantities L in claim 1, we show numerical simulations of the dynamics, Eq. ( 3).From the solution X(t), the volume Ω(X) of the system is determined by the equation of state, Eq. ( 17).To numerically solve Eq. ( 3) with the stoichiometric matrix S in Eq. ( 10), we employ the mass-action kinetics with the local detailed balance condition to specify the flux function J(t) (see Appendix B).
Claim 1: In Fig. 2(b-d), we show numerical trajectories of the volume Ω(X) from three initial conditions for different pressures Π.The color of curves corresponds to each initial condition.The fate of the system changes with Π as follows.When Π = 13 < ϕ * (y min ) = 14.25 (Fig. 2(b)), the volume of the system increases with time from any initial conditions.As Π increases to Π = ϕ * (y min ) = 14.25 (Fig. 2(c)), the trajectory from an initial condition with L = 0 in light blue achieves an equilibrium state, whereas trajectories from the other two with L = 0 continue growing.As Π increases further, Π = 16 > ϕ * (y min ) = 14.25 (Fig. 2(d)), the trajectory with L = 0 in light blue decreases and the system shrinks.By contrast, those from the other two with L = 0 achieve equilibrium states.These numerical results demonstrate the claim 1.
Claim 2: In Fig. 2(e,f), we demonstrate the time evolution of the number X.When Π = ϕ * (y min ) = 14.25 and L = 0, the equilibrium state to which the system converges depends on the initial conditions in M X STO (L = 0).Furthermore, it indeed lies on the ray r Y (y min ), which is given by Eq. ( 21) (see Fig. 2(e)).In our numerical simulation, we compute it as follows.From y min , we obtain the corresponding density )/R T } by using Eq.(31).The ray r Y (y min ) is plotted by the one which includes x min .When Π = 16 > ϕ * (y min ) and L = 52 = 0, the system indeed converges to the unique equilibrium state, irrespective of the initial conditions in M X STO (L = 52) (see Fig. 2(f)).

G. Outline of the derivations of claims 1 and 2, and Corollary 1
The fates of the system stated in claims 1 and 2, and Corollary 1 are summarized in Table I.Before closing this section, we outline their derivations, which will be presented in the subsequent sections.They will be derived in two steps.In the first step, we investigate the existence of an equilibrium state, that is the intersection M Y STO (L) ∩ M Y EQ (μ) in Eq. ( 27).In Sec.III, we will show that the intersection is not empty in the case 2 for The time evolution of X is shown in the number space X for three initial conditions.(e) We set that all the initial conditions have L = 0 and are given by (X 1 0 , X 2 0 , X 3 0 ) = (48, 1, 48) (light blue), (25, 60, 25) (light green), (10, 50, 10) (light purple).Each initial condition is indicated by a circle in the corresponding color.From each initial condition, the curve is constrained in M X STO (L = 0) and finally converges to the square when Π = 14.25.All the squares are located on the red ray r Y (y min ).The black and red circles on this ray indicate the origin X = 0 and We set that the three initial conditions have L = 52 and are given by (X 1 0 , X 2 0 , X 3 0 ) = (148, 25, 200) (orange), (10, 175, 62) (blue), and (1, 10, 53) (purple), respectively.Each of them is indicated by a circle in the corresponding color.The curve is constrained in M X STO (L = 52) and finally converges to the unique point (square).The parameters are fixed as L = 0 and in the case 3 for L = 0 of claim 1, and is empty for the other cases.This will be summarized in Theorem 1.Following Theorem 1, we derive claim 2 from the correspondence in Eq. ( 21) between the number space X and the chemical potential space Y.In the second step, when the intersection is empty, we investigate the landscape of the total entropy function Σ tot to classify whether the system grows or shrinks in Sec.IV.We will show that the system has two possibilities for L = 0: it grows when Σ tot is not bounded above, or it shrinks when the supremum of Σ tot is at the origin X = 0.By contrast, the system always grows for L = 0 when the intersection is empty because Σ tot is not bounded above.This will be summarized in Theorem 2. Mathematical proofs of Theorems 1 and 2, and Corollary 1 will be presented in Sec.V.

III. GEOMETRIC REPRESENTATION OF EQUILIBRIUM STATES
For the derivation of claims 1 and 2, we analyze the accessible regions in the number space X, the density space X and the chemical potential space Y by using a geometric technique.In Sec.III A, we give the assumptions to the full grand potential density ϕ * (y) for the following mathematical treatment.In Sec.III B, we show that the accessible region is restricted by the isobaric condition, which we term the isobaric manifold.This manifold is further partitioned by the conservation laws in Sec.III C.After defining the equilibrium manifold as the candidates of chemical equilibrium states in Sec.III D, we characterize the equilibrium state to which the system converges as an intersecting point of the accessible region and the equilibrium manifold in Sec.III E. Furthermore, by using numerical simulations, we identify the conditions for the existence of the intersecting point, which we state as Theorem 1.
A. Assumptions to ϕ * (y) The thermodynamic property of the system is encoded in the full grand potential density ϕ * (y) (see Eq. ( 20)).
Here, y ∈ Y = R NX is the chemical potential of the confined chemicals, and ϕ * (y) represents the pressure of the system at the state y.In this work, we assume that ϕ * (y) is a smooth and strictly convex function on Y = R NX , and the image of the associated Legendre transformation, is equal to X = R NX >0 .In addition, we assume the following properties; (1) ϕ * (y) strictly increases with y i for an arbitrary fixed {y j } j =i .(2) ∂ i ϕ * (y) → 0 when y i → −∞.They are satisfied in most thermodynamic systems such as the ideal gas in Eq. (31).
Furthermore, we assume that the mapping ∂ϕ * : R NX → R NX >0 is bijective, and we consider that its inverse map is given by ∂ϕ : R NX >0 → R NX (see Eq. ( 13) for ϕ(x)).By using this, the above property (2) is rephrased as Finally, we assume that the pressure of the environment Π is greater than Π min , to guarantee that the volume Ω(X) is uniquely determined (see Appendix C).Here, Π min denotes the minimum pressure that the system can take as The second equality holds from the above property (1).

B. The isobaric manifold and the projection ρX
We remind that the density x is given by the mapping ρ X in Eq. ( 18): The range of the mapping ρ X describes possible states of the density x under the isobaric condition.To calculate the range, we compute the volume Ω by solving the minimization in Eq. ( 14).Its critical equation is given by Therefore, the possible region in X under the isobaric condition is constrained in which we call the isobaric manifold.Since the range of ρ X is represented in Eq. ( 37), we find that the mapping ρ X (X) is a projection of X ∈ X into I X ( Π, μ) along the ray which includes X.
Example 3: For the ideal gas case, by substituting ϕ(x) in Eq. ( 15) into Eq.( 37), we obtain where ñm is the density of the open chemicals in the environment (see Eq. ( 16)).Thus, the isobaric manifold in X represents a simplex for the ideal gas case (Fig. 3(a)).
The mapping ρ X is the projection of X ∈ X = R NX  (a) For the ideal gas case, the isobaric manifold I X ( Π, μ) is represented by a simplex in gray.(b) The mapping ρX (X) gives the projection of X ∈ X into I X ( Π, μ) along the ray which includes X.The thick arrow is pointing the corresponding density x = ρX (X).(c) For the reactions in Eq. ( 9), the stoichiometric manifold M X STO (L) in Eq. ( 40) represents a two dimensional plane with U = (−1, 0, 1).By applying the mapping ρX (X) to each point X ∈ M X STO (L), we get the corresponding manifold, M X STO (L) (see the upper region of the simplex in (d) colored by orange).(d) In the present case, the isobaric manifold I X ( Π, μ) is partitioned by M X STO (L > 0), M X STO (L < 0), and M X STO (L = 0).Here, M X STO (L > 0) and M X STO (L < 0) are given by the orange and the green regions, respectively, and their interface is given by M X STO (L = 0) in the light blue line.(e) The stoichiometric manifold M X STO (L) is characterized by a ray r in the space L of the conserved quantities (see the line below Eq. ( 6) for L).Since dim Ker[S T ] = 1 in the present example, we have three characterizations L = 0, L > 0 and L < 0. In Sec.I of the Supplementary material, we investigate another example with a higher dimension of Ker[S T ].(f) Given a ray r, any plane M X STO (L) for L ∈ r is mapped to the same region on the simplex by the mapping ρX .Any plane with L > 0 (e.g. in orange) is mapped to the upper region of the simplex (see (d)), whereas, any plane with L < 0 (e.g. in green) is mapped to the lower region.For L = 0, the plane in light blue is mapped to the interface between the regions for L > 0 and L < 0 (see the light blue line in (d)).
Finally, by using the mapping ∂ϕ, the isobaric manifold in Y is represented as (39) which is straightforward to interpret that the pressure of the system ϕ * (y) at the chemical potential y is balanced with Π in this manifold.
C. The stoichiometric manifold as partition of the isobaric manifold The accessible subspace of the state X(t) from an initial condition X 0 in Eq. ( 4) is which we call the stoichiometric manifold in X.
By applying the mapping ρ X to M X STO (L), we have the stoichiometric manifold in the density space X as where α > 0 is arbitrary.Using the mapping ∂ϕ, the stoichiometric manifold in the chemical potential space Y is given by (42) Example 4: For the specific stoichiometric matrix in Eq. ( 10), M X STO (L) is schematically depicted in Fig. 3(c,d) by applying the mapping ρ X to M X STO (L).
For further investigations, we define the following manifold in X as and that in Y as These manifolds correspond to the stoichiometric manifolds with conserved quantities L under the isochoric condition [78].Hereafter, we term the isochoric stoichiometric manifolds.Using these manifolds, the stoichiometric manifold under isobaric condition is represented as follows: and For arbitrary c > 0, we find M X STO (cL) = M X STO (L) and, correspondingly, M Y STO (cL) = M Y STO (L).Hence, each ray r in the space L of conserved quantities characterizes these manifolds (see also Example 5 below and Fig. 3(df)).To be more precise, for any L, L ′ ∈ r, M X STO (L) = M X STO (L ′ ) holds.We denote the stoichiometric manifold corresponding to a ray r as M X STO (L ∈ r).When L = 0, the stoichiometric manifold in X is written as Since the dimension of α in Eq. ( 45) vanishes for In addition, for any given L = 0, we find that MX STO (L/α) → MX STO ( L = 0) when α → ∞.Thus, the stoichiometric manifold M X STO (L = 0) forms an interface of all the stoichiometric manifolds M X STO (L = 0).These properties hold similarly in Y by replacing the superscript X with Y.
The above statements are summarized as follows.
Lemma 1 Let r denote a ray in the space L of conserved quantities.For every r, the corresponding stoichiometric manifold M X STO (L ∈ r) exists and the set of the stoichiometric manifolds partitions the isobaric manifold I X ( Π, μ).In addition, their interface is given by the stoichiometric manifold M X STO (L = 0).The isobaric manifold in the chemical potential space Y is similarly partitioned by stoichiometric manifolds (replace the superscript X with Y).
Example 5: In Fig. 3(d-f), we schematically depict the partition of the isobaric manifold I X ( Π, μ) by the stoichiometric manifolds M X STO (L ∈ r) for the ideal gas case.

D. The equilibrium manifold
In this subsection, we confirm that the candidates of the equilibrium states are given by Eq. (23).
The equilibrium states of the chemical reaction dynamics are given by the maxima of the total entropy with respect to the extent of reaction Ξ.From the differentiation of Eq. ( 12) and using Eq. ( 14), we have the critical equation as where X = X 0 + SΞ in Eq. ( 4) (see Appendix D).
Thus, at any equilibrium states, the density x satisfies By using the correspondence y i = ∂ i ϕ(x), it is represented as which confirms the balance of chemical potentials between reactants and products at the chemical equilibrium (see Eq. ( 22)).
The set of solutions to Eq. ( 49) represents the candidates of the equilibrium states in the density space X , which we call an equilibrium manifold: and in the chemical potential space Y as We note that the equilibrium manifold is an affine subspace in Y, whereas it is generally curved in X .In addition, using the basis matrix U of Ker[S T ] and solving Eq. ( 50), the manifold in Eq. ( 52) is rewritten as where y P is a particular solution to Eq. ( 50) and h = {h l } represents a coordinate of Ker[S T ].

E. The equilibrium state as an intersection
If the intersection M Y STO (L) ∩ M Y EQ (μ) is not empty, the system converges to the equilibrium state.We can compute the equilibrium state in Y by using Eqs.( 42) and (52).Similarly, we obtain the corresponding state in X as M X STO (L) ∩ M X EQ ( T , μ) by using Eq. ( 41) and (51).In this subsection, we numerically confirm the convergence in X by the dynamics in Eq. (3) using the specific setup in Sec.II F.
We first consider the equilibrium manifold.The points in M Y EQ (μ) are written in Eq. ( 30).Here, the second component y 2 is constant p := μ1 − μ2 .Thus, M Y EQ (μ) is a one-dimensional line located on the plane y 2 = p in Y (see the left panel of Fig. 4 (a) for the plane y 2 = p, and the red line in the left panel of Fig. 4 (b) for M Y EQ (μ)).By applying the mapping ∂ϕ * , we can compute M X EQ ( T , μ).For the ideal gas case in Eq. ( 31), the mapping is represented as ∂ i ϕ * (y) = exp{(y i − ν o i )/R T }.Thus, the points in M X EQ ( T , μ) are computed as where h ∈ R plays a role of the coordinate of M X EQ ( T , μ).We note that the second component x 2 is also constant q := ∂ϕ * (p) = exp{(μ 1 − μ2 − ν o 2 )/R T }.Thus, M X EQ ( T , μ) is a one-dimensional curve located on the plane x 2 = q.(see the red curve in the right panel of Fig. 4(b)).
We next consider whether the intersection I X ( Π, μ) ∩ M X EQ ( T , μ) is empty or not.As the pressure Π increases, the number of intersecting points change as follows.When Π = 13 < ϕ * (y min ) = 14.25, the intersecting point I X ( Π, μ) ∩ M X EQ ( T , μ) does not exist (see the right panel of Fig. 4(b)).When Π = ϕ * (y min ), the intersection appears as the unique contact point corresponding to y min .When Π = 16 > ϕ * (y min ), the intersection consists of two points.The above situation holds similarly in Y (see the left panel of Fig. 4(b)).
We further depict the intersection I X ( Π, μ) ∩ M X EQ ( T , μ) by the three dimensional plot in Fig. 5. From Fig. 3(d), I X ( Π, μ) is partitioned by the stoichiometric manifolds.When Π = ϕ * (y min ) = 14.25 (Fig. 5(b)), the intersection lies on M X STO (L = 0).When Π = 16 > ϕ * (y min ), the two intersecting points belong to each of M X STO (L) for L > 0 and L < 0, respectively (see Fig. 5(c)).Thus, for a given L, the intersecting point is uniquely determined and it depends on a ray in the space L of the conserved quantities (see Fig. 3(e) and Lemma 1).
In Fig. 5(b,c), we also show the time evolution of the system in the numerical simulation.Indeed, the system converges to the intersection M X STO (L) ∩ M X EQ ( T , μ).The above numerical result is an instance of the following theorem.

Theorem 1 The system converges to the intersecting point M Y
STO (L) ∩ M Y EQ (μ) if it exists.Its existence is classified as follows: where the intersecting point y EQ is uniquely determined for a given L, and it depends on a ray in the space L of conserved quantities.
A mathematical proof of this theorem will be shown in Sec.V.
From Theorem 1, we can derive claim 2 as follows.For the case 1 of claim 2 in which ϕ * (y min ) − Π = 0 (a) For the ideal gas case, the isobaric manifold I Y ( Π, μ) in the chemical potential space Y is a level set for ϕ * (y) (left), whereas I X ( Π, μ) represents a simplex in the density space X (right).We consider the planes y2 = p = μ1 − μ2 in Y and x 2 = q = exp{(μ1 − μ2 − ν o 2 )/R T } in X , respectively.The black curve in Y and the line in X indicate the cross section of the isobaric manifolds by the planes, respectively.(b) In the left panel, we plot M Y EQ (μ) (the red line) and I Y ( Π, μ) (gray curves) on the plane y2 = p.In the right panel, M X EQ ( T , μ) (the red curve) and I X ( Π, μ) (gray lines) are plotted on the plane x 2 = q.In Y (left), the three gray curves correspond to the isobaric manifolds I Y ( Π, μ) with the pressures Π = 13 (dotted), 14.25 (solid) and 16 (dashed), respectively.When Π = 13, the dotted curve does not have the intersection with M Y EQ (μ) (red line).When Π = 14.25, the solid curve has the unique intersecting point with the red line, denoted by the light blue square.It is located at y min = (−0.981,−0.288, 1.099), where y min has been obtained in Fig. 2. When Π = 16, the dashed curve has two intersecting points with the red line.The heatmaps represent the values of ϕ * (y) − ϕ * (y min ).In the right panel, we show the corresponding manifolds in X .Here, the light blue square is located at xmin = (3.00,5.25, 3.00), which is obtained by and L = 0, Theorem 1 states that the intersection M Y STO (L) ∩ M Y EQ (μ) consists of y min .From Eq. ( 21), this y min is mapped to a ray r Y (y min ) in the number space X.In addition, M X STO (L = 0) from Eq. ( 40) indicates that αX ∈ M X STO (L = 0) if X ∈ M X STO (L = 0) for arbitrary α > 0. This means that all the points on the ray r Y (y min ) is in M X STO (L = 0).Thus, all the points on the ray are equilibrium states and the system converges to one of them.The state is not uniquely determined only by the entropy function and depends on the initial conditions and the functional form of J(t).For the case 2 of claim 2 in which ϕ * (y min ) − Π < 0 and L = 0, Theorem 1 states that the intersection M Y STO (L) ∩ M Y EQ (μ) FIG. 5. (a) When Π = 13 < ϕ * (y min ) = 14.25, M X EQ ( T , μ) (red curve) and I X ( Π, μ) (gray simplex) do not have an intersecting point.The black line represents I X ( Π, μ) on the plane x 2 = q (see also Fig. 4).(b) When Π = 14.25 = ϕ * (y min ), the intersection consists of the unique point xmin.It lies on the stoichiometric manifold M X STO (L = 0) (the light blue line on the simplex).In numerical simulations, the time evolution of the system is also shown in a light blue line with arrows.From the initial condition (circle), it converges to the point xmin (square).(c) When Π = 16 > ϕ * (y min ), the intersection consists of two points (orange and green squares).They are located in the upper and the lower regions of the simplex, which correspond to M X STO (L) for L > 0 and L < 0, respectively (see Fig. 3(d)).We also show the time evolution of the system from the two initial conditions (circles).They respectively converge to the intersecting points (squares).In the numerical simulation, the initial conditions are given by (X 1 0 , X 2 0 , X exists as y EQ .This also corresponds to a ray r Y (y EQ ) in the number space X.Since M X STO (L = 0) is an affine subspace which does not include the origin X = 0, the intersection between M X STO (L = 0) and an arbitrary ray is represented by a point.Thus, the equilibrium state is uniquely determined in this case.
Theorem 1 states the cases when the system converges to an equilibrium state in our claim 1.In the next section, we will investigate the case when the system does not converge to an equilibrium state, that is

IV. FORM OF THE TOTAL ENTROPY FUNCTION
When the intersection M Y EQ (μ) ∩ M Y STO (L) is empty, the system cannot converge to an equilibrium state, and it is expected that the system grows or shrinks.In this section, we observe if the system grows or shrinks by numerical plots of the total entropy function.
If M X STO (L) ∩ M X EQ ( T , μ) = ∅, the total entropy function Σ tot (Ξ) can not have the maximum inside of M X STO (L) (see Sec. III D).Thus, the functional form of Σ tot (Ξ) has the following two possibilities: (1) Σ tot (Ξ) is not bounded above.(2) Σ tot (Ξ) is bounded above, but the supremum exists on the boundary of M X STO (L).In the first case, we find |Ξ(t)| → ∞ for t → ∞, because the second law requires that the system must climb up the unbounded landscape of the total entropy function.
Since we assume dim Ker[S] = 0, we get |X| → ∞ for |Ξ| → ∞.In addition, Ω(X) → ∞ for |X| → ∞ (see Appendix E).Thus, if Σ tot (Ξ) is not bounded above, the volume Ω grows.In the second case, we can further show that the supremum is possible only at the origin X = 0 (see Sec. V B).Thus, the volume Ω is shrinking and finally vanishes.
The above numerical result demonstrates the following theorem.
3. If ϕ * (y min ) − Π < 0 and L = 0, the supremum of Σ tot (Ξ) is at the origin of the number space X.
A mathematical proof of this theorem will be shown in Sec.V. Theorem 1 and 2 lead to the claim 1 and we find the fate of the system under isobaric condition.
Before closing this section, we confirm the claim 2 in the left panel of Fig. 6(b) and the right panel of Fig. 6(c).In the former panel, the red ray represents the set of the maxima of Σ tot (Ξ), which is r Y (y min ).We can see the system converges to a point on the ray.In the latter case, the maximum of Σ tot (Ξ) is the unique point and the system converges to it.

V. PROOF OF THE THEOREMS
In this section, we prove Theorems 1 and 2, and Corollary 1.For L = 0, Σ tot (Ξ) has maxima which form the red ray.The system converges to a point (the light blue square) on the ray.For L = 0, the heatmap shows unbounded landscape and the system goes to infinity.(c) For L = 0, Σ tot (Ξ) has the supremum at the origin X = 0 and the system is approaching this point.For L = 0, Σ tot (Ξ) has the unique maximum (the orange square) and the system converges to the point.See the caption in Fig. 2 for specific values of the parameters.

A. Proof of Theorem 1
We first provide a general proof valid for arbitrary dimensions.After concluding the proof, we schematically illustrate it for three dimensional space of Y in Fig. 7.
To investigate the existence of the intersection M Y STO (L) ∩ M Y EQ (μ), we introduce Birth's trajectory.Let y B ( L) represent the unique intersecting point between the isochoric stoichiometric manifold MY STO ( L) in Eq. ( 44) and the equilibrium manifold M Y EQ (μ) in Eq. ( 52): The uniqueness of y B ( L) is guaranteed by Birch's theorem and the point is known as Birth's point [12,17,78,83,84] (see Appendix F).We define Birch's trajectory by a collection of Birch's point along a ray in the space of conserved quantities.For given L = 0, it is defined as We find that y B (L/α) → y B (0) for α → ∞.Furthermore, from Eq. (F5) in Appendix F, we get Thus, the point y min is located at one of the end points of Birch's trajectory, but is not included in the trajectory.For this Birch's trajectory, the following lemma is satisfied.
Lemma 2 ϕ * (y) is a strictly increasing function along Birch's trajectory from the starting point y min .
A proof of this lemma is shown in Appendix G.
From Lemma 2, we can determine the existence of the intersection M Y STO (L) ∩ M Y EQ (μ) by the position of the starting point y min .When ϕ * (y min )− Π > 0, the starting point y min is located in the superlevel set: {y|ϕ * (y) > Π}.Taking into account the fact that I Y ( Π, μ) is the level set, the intersecting point M Y STO (L)∩M Y EQ (μ) does not exist for any L. When ϕ * (y min ) − Π = 0, the starting point y min is on I Y ( Π, μ).For L = 0, the intersecting point uniquely exists as y min because of Eq. (59).For L = 0, Lemma 2 indicates ϕ * (y) > ϕ * (y min ) = Π at any point y on Birch's trajectory b(L).Thus, from Eq. ( 58), the intersecting point does not exist.When ϕ * (y min ) − Π < 0, the starting point y min is located in the sublevel set: {y|ϕ * (y) < Π}.For L = 0, the intersecting point does not exist.For L = 0, the intersecting point uniquely exists.
This concludes the proof of Theorem 1.
In Fig. 7, we provide schematic illustrations in three dimensional space of Y for the case ϕ * (y min ) − Π < 0. To demonstrate the generality of the proof, we visualize two cases with dim Ker[S T ] = 2 and dim Ker[S T ] = 1 in the panels (a) and (b), respectively.From Eq. ( 6), the case in (a) has two conserved quantities, whereas the case in (b) has one conserved quantity.We also have dim Im[S] = 1 and dim Im[S] = 2 for the two cases, respectively, because dim Y = dim Ker[S T ]+dim Im[S] and dim Y = 3.The dimensions for the example in Eqs ( 9) and ( 10) correspond to the case (b), whereas the case (a) corresponds to the one investigated in Sec.I of the Supplementary material.
In the former case (Fig. 7(a)), the equilibrium manifold M Y EQ (μ) is represented by the two dimensional plane in red, whereas the isochoric stoichiometric manifold MY STO ( L) is given by a one-dimensional curve in light blue because dim M Y EQ (μ) = dim Ker[S T ] and dim MY STO (L/α) = dim Im[S].Accordingly, the black circles represent the intersecting points y B ( L) in Eq. ( 55) for L = L/α with varying α.The green curve indicates Birch's trajectory b(L) in Eq. ( 56), and one of the end points is located at y min in Eq. ( 57).The isobaric manifold I Y ( Π, μ) is represented by the gray surface, which is a level set for ϕ * (y).When ϕ * (y min ) − Π < 0, y min is in the sublevel set {y|ϕ * (y) < Π}.The stoichiometric manifold M Y STO (L) is given by the orange curve from Eq. ( 46).Thus, the intersection M Y STO (L) ∩ M Y EQ (μ) is represented by the red circle, which is also the intersecting point between b(L) and I Y ( Π, μ) in Eq. ( 58).
In the latter case (Fig. 7(b)), M Y EQ (μ) is represented by the red line, whereas MY STO ( L) is given by the two dimensional surface in light blue.The intersecting point y B ( L) is indicated by the black circle in Eq. ( 55) for L = L/α with a given α.In this case, the Birch's trajectory b(L) in Eq. ( 56) overlaps with the red line, and one of the end points is located at y min in Eq. (57).
stoichiometric manifold M Y STO (L) is given by the orange surface from Eq. ( 46) [85].Thus, the intersection M Y STO (L) ∩ M Y EQ (μ) is represented by the red circle.If M X STO (L) ∩ M X EQ ( T , μ) = ∅, the system never relaxes to the equilibrium state and have the following two possibilities: (1) Σ tot (Ξ) is not bounded above, and |Ξ(t)| → ∞ with time.As we noted in Sec.IV, the volume Ω grows in this setup.(2) Σ tot (Ξ) is bounded above and the supremum exists on the boundary of M X STO (L).First, we consider the case L = 0, where M X STO (L = 0) does not contact with the origin.In this case, Σ tot (Ξ) is not bounded above.To prove that, we will show that the supremum of Σ tot (Ξ) does not exist on the boundary of M X STO (L = 0) as follows.
The extent of reaction Ξ provides a coordinate on the stoichiometric manifold M X STO (L).Since X = R NX >0 , the domain of Ξ has a boundary.We denote the coordinate of a point on a boundary by Ξ B .Here, at least, one component of the vector X B := X 0 + SΞ B is zero.We define the index set by I B (X B ) = {i|X i B = 0}.Since the origin X = 0 is excluded, I B (X B ) never has all the indices.Also, the density x B = X B /Ω(X B ) has zero components as x i B = 0 for i ∈ I B (X B ).To investigate the gradient of Σ tot (Ξ) on the boundary, we calculate the thermodynamic force as If the supremum of Σ tot (Ξ) does not exist on the boundary, we can find a direction from the boundary to the interior of M X STO (L) such that Σ tot (Ξ) increases.To be more precise, by using f r (Ξ), this statement is mathematically phrased as follows: If the supremum of Σ tot (Ξ) does not exist on the boundary, a vector j = {j r } exists such that f r (Ξ B )j r > 0 and S i r j r > 0 for all the indices i ∈ I B .In the next paragraph, we prove the above statement.
From the assumptions to ϕ * (y) in Sec.III A, for every i, y i = ∂ i ϕ(x) → −∞ when x i → 0. Thus, ∂ i ϕ(x B ) → −∞ for all i ∈ I B .Choose a vector j satisfying S i r j r > 0 for i ∈ I B .Then, we can show that, when Ξ → Ξ B , f r (Ξ)j r → ∞ from Eq. ( 60).This indicates that the supremum of Σ tot (Ξ) does not exist on the boundary.
Second, we consider the case L = 0.By employing rays on M X STO (L = 0), we investigate whether Σ tot (X) is bounded above or not.
Since we have assumed dim Ker[S] = 0, we can rewrite Σ tot (Ξ) as a function of X.Using a particular solution y P to Eq. ( 50), we get + ΠΩ(X) where (see Appendix H for the details of the calculation).Here, D Y [y P ||y] is the Bregman divergence [86][87][88] induced by ϕ * (y), which is defined as and y(X) = ∂ϕ • ρ X (X).Although Eq. ( 61) is apparently a function of the particular solution y P , the value of Σ tot does not depend on the choice of y P as shown in Appendix I. Since only the first term Ω(X)K Y (y P ; y(X))/ T depends on X, we consider its landscape.
The stoichiometric manifold with L = 0, M X STO (L = 0), can be described by the collection of rays in X because it contacts with the origin X = 0. Since a ray in the number space X corresponds to a point in the chemical potential space Y as in Eq. ( 21), the value of K Y (y P ; y) is constant on each ray.In addition, Ω(X) is increasing along each ray from the origin.Therefore, by considering the sign of K Y (y P ; y), we can determine if Σ tot (X) increases along each ray in M X STO (L = 0).If K Y (y P ; y) is positive for a ray, Σ tot (X) increases along the ray.If K Y (y P ; y) = 0, Σ tot (X) is constant along the ray.If K Y (y P ; y) is negative for a ray, Σ tot (X) decreases along the ray.
When ϕ * (y min ) − Π > 0, we can find a point y ∈ M Y STO (L = 0) such that K Y (y P ; y) is positive (see Appendix J).If such a point y ∈ M Y STO (L = 0) exists, we can consider a corresponding ray in X to the point y.Along the ray, Σ tot (X) increases and, thus, is not bounded above.When ϕ * (y min ) − Π = 0, K Y (y P ; y) is negative for any y ∈ M Y STO (L = 0) except y = y min , and K Y (y P ; y min ) = 0. Thus, the maxima of Σ tot (X) forms the ray corresponding to y min , that is r Y (y min ) in Eq. (21).When ϕ * (y min ) − Π < 0, we can show that any y ∈ M Y STO (L = 0) gives negative K Y (y P ; y) (see Appendix K).Thus, Σ tot (X) decreases along any ray and the supremum is located at the origin X = 0.
Summarizing the cases for L = 0 and L = 0, we obtain Theorem 2.

C. Proof of Corollary 1
When S is regular, the accessible region covers the whole space of X, that is M X STO = X.Thus, M X STO can be described by the collection of rays in X.This situation coincides with the one for L = 0 in the previous subsection.Therefore, by employing a similar proof, we can classify the fate of the system as Corollary 1 (see also [19]).

VI. SUMMARY AND DISCUSSIONS
In this work, we have considered chemical thermodynamics of growing CRSs with stoichiometric conservation laws in which the stoichiometric matrix S has a nontrivial left kernel space.By clarifying the conditions for the fate of the system to grow, shrink, or equilibrate, we show that the existence of conservation laws qualitatively alter the fate of the system.It is emphasized again that our results are derived by general thermodynamic and stoichiometric structures without assuming any specific thermodynamic potentials or reaction kinetics; i.e., they are obtained based solely on the second law of thermodynamics and conservation laws.
Since we are mainly interested in the growth of the system, we have assumed the productivity of the stoichiometric matrix S: Im[S] ∩ R NX >0 = ∅ (see Eq. ( 7)).This condition allows the time evolution of X(t) so as to perpetually increase the number of all confined chemicals while satisfying the conservation laws.In the following two paragraphs, we comment on the fate of the system with the other classes of the stoichiometric matrix.Here, the complement of the productive S is defined as Im[S] ∩ R NX >0 = ∅.We divide it into two cases: . In the first case, (1) Im[S] ∩ R NX ≥0 = ∅, a vector v exists in Ker[S T ] such that all the components of v are positive.As a result, the conserved quantity v i X i = v i X i 0 implies that the total mass of the confined chemicals are conserved, and the system must converge to the equilibrium state, i.e., the growth of the system is not possible.We refer to this type of S as non-productive.A numerical example for this case is shown in Sec.III of the Supplementary material.
In the second case, (2) Im[S] ∩ R NX ≥0 = ∅ and Im[S] ∩ R NX >0 = ∅, a vector v exists in Ker[S T ] such that all the components of v are non-negative, but Ker[S T ] does not have a vector v ′ such that all the components of v ′ are positive.This means that the mass conservation law exists for a proper subset of the confined chemicals.Thus, the volume growth of the system can still be allowed with the increase of the confined chemicals that do not participate in the mass conservation law.Although such a situation may not be biologically relevant, we refer to this type of S as semi-productive.A numerical example for this case is shown in Sec.IV of the Supplementary material.According to the example, it is expected that the claims 1 and 2 still hold for the semi-productive cases.
We can possibly correlate our theoretical findings with empirical observations.Our findings suggest that the fate changes with the initial conditions (conserved quantities L) even in the same environment.We especially find that, as shown in Table I, a CRS does not vanish and is more likely to grow when it has non-zero L than when it has L = 0 or a regular S. If we suppose a population of the CRSs with varying initial conditions, and consider their fates in a slowly changing environment, surviving CRSs can be selective depending on L, where a CRS with non-zero L might be advantageous.Identifying a possible mechanism to survive in a varying environment is biologically essential, and our findings may help explain the outcome.
Nevertheless, several extensions of our theory remain to be essential as future work to understand actual experimental situations of growing protocells or biological cells.
The first extension is to consider the case when the system may relax to a state that continuously produces entropy with constant volume, namely, the conventional nonequilibrium steady state (NESS) [9,10,[12][13][14][15][16][17].It is possible when the matrix S has a nontrivial right kernel space, i.e., the assumption in Eq. ( 5) is not satisfied.It is a major challenge to understand how the growth of the system can be characterized and realized in such situations, and compatible with the NESS.
The second one is to consider a different hierarchy of the timescales in which the relaxation of the variables (E, Ω, N ) may not be rapid compared to the timescale of chemical reactions.Although the present work assumes that the chemical reactions are slow (see the assumptions made above Eq.( 12)), there may be cases in which the timescale of other processes is slow and/or comparable with that of the chemical reactions.A typical example is when we consider cell transport processes [90].In this case, the timescale of material exchanges between the system and the environment is relevant, and the processes may also couple with chemical reactions.In addition, the osmotic gradient (pressure difference) can play a driving force to change the volume, whereas isosmotic volume changes are also commonly known and widely investigated as a relevant process by alterations in intracellular solute content [91,92].It is of interest to see how the conditions for the growth of the system change with relevant processes of the thermodynamic setup.
The third one is to consider explicitly the effect of the membrane molecules.Although we neglect the tension of the membrane in this paper (see Fig. 1), it could be relevant and changes with time because the membrane molecules themselves are produced and supplied by the CRSs in biological cells.It is important to investigate if and how the effect changes the results by extending our theoretical framework.
Alternatively, it would also be promising to apply our framework to fit into an experimental technique using membrane-free compartments such as droplets based on liquid-liquid phase separation [93].Since its nature, a rigorous thermodynamic treatment would require an extended framework in which the entropy function is not strictly concave (see the assumption made below Eq. ( 2)).
In all the above extensions, the present work provides basic results to clarify the role of conservation laws to the growth of protocells and biological cells.
where N = (N 1 , N 2 ) denotes the number of B = (B 1 , B 2 ) in the system.The rate constants w r + and w r − satisfy the local detailed balance condition [11,14,16,78,84]: For the ideal gas case, the density N/Ω of the open chemicals in the system is constant and fixed to that of the environment ñ.Also, the volume Ω is given by the equation of state in Eq. ( 17).Thus, we can rearrange Eq. (B1) as where we absorb the constant densities of the open chemicals into the rate constants as ŵr + and ŵr − .For these effective rate constants, the local detailed balance condition in Eq. (B2) is written as log ŵr In the number space X, we can show that the volume function Ω(X) satisfies the homogeneity: Ω(αX) = αΩ(X) for α > 0. Thus, for a given X r , we have ∂Ω(αX r ) ∂α = ∂ ∂α {αΩ(X r )} = Ω(X r ). (E1) This implies that the volume function Ω(X) increases with the rate Ω(X r ) along the ray including X r .Therefore, the volume Ω(X) goes to infinity for |X| → ∞. that the stoichiometric matrix S is productive, therefore, M Y STO (L = 0) = ∅ (see Appendix J).Consider any point y ∈ M Y STO (L = 0).At the point y, we can consider the tangent plane of the level set {y|ϕ * (y) = Π}, which divides the space of Y into two regions.Since K Y (y P ; y) represents the inner product of the two vectors ∂ϕ * (y) and (y P − y) in Eq. (J1), the sign of K Y (y P ; y) is determined from the sides of the tangent plane to which the two vectors point.We also note that any point y P ∈ M Y EQ (μ) is on the one side of the tangent plane because the plane is parallel to the equilibrium manifold M Y EQ (μ) (see Fig. 9).When ϕ * (y min ) − Π < 0, the equilibrium manifold M Y EQ (μ) (the red line in Fig. 9) is located both in the sublevel and the superlevel sets.From the convexity of the function ϕ * (y), any point in the sublevel set is located on the opposite side of the tangent plane to the orientation of the normal vector ∂ϕ * (y) at the point y (see Fig. 9).Since any point y ∈ M Y EQ (μ) is on the one side of the tangent plane, we can show that the particular solution y P ∈ M Y EQ (μ) is located on the opposite side of the tangent plane to the orientation of ∂ϕ * (y).Thus, for y P ∈ M Y EQ (μ), the inner product between ∂ϕ * (y) and (y P i − y i ) is negative.Therefore, for y ∈ M Y STO (L = 0), K Y (y P ; y) is always negative.The energy exchange rate with the environment Sec.II A/Fig.

JΩ(t)
The volume expansion rate Sec.II A/ The partial grand potential density Sec.II A/Eq.( 13) y The Legendre dual variable of x Sec.II B ϕ * (y) The full grand potential density Sec.II B/Eq.( 20) Πmin The minimum pressure that the CRS can take Sec.III A/Eq.( 34) fr(Ξ) The thermodynamic force Sec.V B/Eq.(60)

FIG. 1 .
FIG. 1. Diagrammatic representation of open CRSs.The chemical reactions occur with the reaction flux J(t) = {J r (t)}, the rth reaction of which is represented as the chemical equation at the bottom.Here, A = {Ai} is the label of the confined chemicals, and B = {Bm} is the label of the open chemicals, which can move across the membrane with the diffusion flux JD(t) = {J m D (t)}.In addition, JE(t) denotes the energy exchange rate with the environment, and JΩ(t) is the volume expansion rate.Let X = {X i } and N = {N m } denote the numbers of the confined and the open chemicals in the system, respectively.Also, (S+) i r and (O+) m r denote stoichiometric coefficients of the reactants in the rth reaction, whereas (S−) i r and (O−) m r are the ones of the products.The stoichiometric matrices are given as S i r = (S−) i r − (S+) i r and O m r = (O−) m r − (O+) m r .For theoretical simplicity, we ignore the tension of the membrane and assume that the membrane never bursts.The stoichiometric matrix S = {S i r } has only a left kernel space.By representing the basis matrix of Ker[S T ] by U , L = U X is conserved under the chemical reaction dynamics.
) Here, three confined chemicals A = (A 1 , A 2 , A 3 ) and two open chemicals B = (B 1 , B 2 ) are involved in two reactions R 1 and R 2 [77].The stoichiometric matrices S and O for the confined and open chemicals, respectively, are

FIG. 3 .
FIG. 3.(a) For the ideal gas case, the isobaric manifold I X ( Π, μ) is represented by a simplex in gray.(b) The mapping ρX (X) gives the projection of X ∈ X into I X ( Π, μ) along the ray which includes X.The thick arrow is pointing the corresponding density x = ρX (X).(c) For the reactions in Eq. (9), the stoichiometric manifold M X STO (L) in Eq. (40) represents a two dimensional plane with U = (−1, 0, 1).By applying the mapping ρX (X) to each point X ∈ M X STO (L), we get the corresponding manifold, M X STO (L) (see the upper region of the simplex in (d) colored by orange).(d) In the present case, the isobaric manifold I X ( Π, μ) is partitioned by M X STO (L > 0), M X STO (L < 0), and M X STO (L = 0).Here, M X STO (L > 0) and M X STO (L < 0) are given by the orange and the green regions, respectively, and their interface is given by M X STO (L = 0) in the light blue line.(e) The stoichiometric manifold M X STO (L) is characterized by a ray r in the space L of the conserved quantities (see the line below Eq. (6) for L).Since dim Ker[S T ] = 1 in the present example, we have three characterizations L = 0, L > 0 and L < 0. In Sec.I of the Supplementary material, we investigate another example with a higher dimension of Ker[S T ].(f) Given a ray r, any plane M X STO (L) for L ∈ r is mapped to the same region on the simplex by the mapping ρX .Any plane with L > 0 (e.g. in orange) is mapped to the upper region of the simplex (see (d)), whereas, any plane with L < 0 (e.g. in green) is mapped to the lower region.For L = 0, the plane in light blue is mapped to the interface between the regions for L > 0 and L < 0 (see the light blue line in (d)).

FIG. 4 .
FIG. 4.(a) For the ideal gas case, the isobaric manifold I Y ( Π, μ) in the chemical potential space Y is a level set for ϕ * (y) (left), whereas I X ( Π, μ) represents a simplex in the density space X (right).We consider the planes y2 = p = μ1 − μ2 in Y and x 2 = q = exp{(μ1 − μ2 − ν o2 )/R T } in X , respectively.The black curve in Y and the line in X indicate the cross section of the isobaric manifolds by the planes, respectively.(b) In the left panel, we plot M Y EQ (μ) (the red line) and I Y ( Π, μ) (gray curves) on the plane y2 = p.In the right panel, M X EQ ( T , μ) (the red curve) and I X ( Π, μ) (gray lines) are plotted on the plane x 2 = q.In Y (left), the three gray curves correspond to the isobaric manifolds I Y ( Π, μ) with the pressures Π = 13 (dotted), 14.25 (solid) and 16 (dashed), respectively.When Π = 13, the dotted curve does not have the intersection with M Y EQ (μ) (red line).When Π = 14.25, the solid curve has the unique intersecting point with the red line, denoted by the light blue square.It is located at y min = (−0.981,−0.288, 1.099), where y min has been obtained in Fig.2.When Π = 16, the dashed curve has two intersecting points with the red line.The heatmaps represent the values of ϕ * (y) − ϕ * (y min ).In the right panel, we show the corresponding manifolds in X .Here, the light blue square is located at xmin = (3.00,5.25, 3.00), which is obtained by x i min = ∂ i ϕ * (y min ) = exp{(y min FIG. 5. (a) When Π = 13 < ϕ * (y min ) = 14.25, M XEQ ( T , μ) (red curve) and I X ( Π, μ) (gray simplex) do not have an intersecting point.The black line represents I X ( Π, μ) on the plane x 2 = q (see also Fig.4).(b) When Π = 14.25 = ϕ * (y min ), the intersection consists of the unique point xmin.It lies on the stoichiometric manifold M X STO (L = 0) (the light blue line on the simplex).In numerical simulations, the time evolution of the system is also shown in a light blue line with arrows.From the initial condition (circle), it converges to the point xmin (square).(c) When Π = 16 > ϕ * (y min ), the intersection consists of two points (orange and green squares).They are located in the upper and the lower regions of the simplex, which correspond to M X STO (L) for L > 0 and L < 0, respectively (see Fig.3(d)).We also show the time evolution of the system from the two initial conditions (circles).They respectively converge to the intersecting points (squares).In the numerical simulation, the initial conditions are given by (X 1 0 , X 2 0 , X30 ) = (48, 1, 48) (light blue circle), (148, 25, 200) (orange circle) and (90, 5, 20) (green circle).

FIG. 6 .
FIG. 6.The entropy landscapes Σ tot (Ξ) are shown for (a) Π = 13 < ϕ * (y min ) = 14.25,(b) Π = ϕ * (y min ) and (c) Π = 16 > ϕ * (y min ).The left and right panels show the cases L = 0 and L = 52 = 0, respectively.The heatmaps represent the values of Σ tot (Ξ).Here, we note that the domain of (Ξ 1 , Ξ 2 ) is restricted by X ∈ R N X >0 .The boundaries at which X 2 = 0 and X 1 = 0 are represented by dashed and dotted black lines, respectively.The black circles in the left panels denote the origin X = 0.The time evolutions of the system are shown by light blue and orange curves for the two initial conditions: (left) X0 = (48, 1, 48) with L = 0 and (right) X0 = (148, 25, 200) with L = 52.(a) The heatmaps in both panels indicate unbounded landscapes and the system goes to infinity.(b)For L = 0, Σ tot (Ξ) has maxima which form the red ray.The system converges to a point (the light blue square) on the ray.For L = 0, the heatmap shows unbounded landscape and the system goes to infinity.(c) For L = 0, Σ tot (Ξ) has the supremum at the origin X = 0 and the system is approaching this point.For L = 0, Σ tot (Ξ) has the unique maximum (the orange square) and the system converges to the point.See the caption in Fig.2for specific values of the parameters.

FIG. 7 .
FIG. 7. Schematic illustration of the manifolds in Y for the case ϕ * (y min ) − Π < 0. In this illustration, we set dim Y = dim Ker[S T ] + dim Im[S] = 3.The conserved quantities L are generally represented by L = {L l = U l i X i 0 } (l = 1, ..., dim Ker[S T ]).The panels (a) and (b) represent the case for dim Ker[S T ] = 2 and dim Im[S] = 1, and the case for dim Ker[S T ] = 1 and dim Im[S] = 2, respectively.We note that dim M Y EQ (μ) = dim Ker[S T ] and dim MY STO (L/α) = dim Im[S].The dimensions in the case (b) correspond to the illustrative example of chemical reactions in Eqs.(9) and (10), whereas those in the case (a) correspond to the one in Sec.I of the Supplementary material.(a) In this case, dim M Y EQ (μ) = 2 and dim MY STO (L/α) = 1.Thus, M Y EQ (μ) is denoted by the red plane, and the light blue curves represent MY STO (L/α) for α ∈ (0, ∞).The gray surface is I Y ( Π, μ).From Eq. (46), M Y STO (L) is given by the orange curve.Equation (56) enables us to depict Birch's trajectory b(L) as the green curve.From Eq. (58), the red circle denotes the intersecting point M Y STO (L) ∩ M Y EQ (μ) to which the system converges as the equilibrium state.(b) In this case, we consider dim M Y EQ (μ) = 1 and dim MY STO (L/α) = 2. Thus, M Y EQ (μ) is denoted by the red line, and MY STO (L/α) represents the two dimensional manifold in light blue.Thus, M Y STO (L) is given by the orange surface.We note that Birch's trajectory b(L) overlaps with the red line of M Y EQ (μ).The red circle represents the equilibrium state.

FIG. 8 .
FIG. 8. Illustration of the proof for the existence of y ∈ M Y STO (L = 0) such that K Y (y P ; y) > 0 when ϕ * (y min ) − Π > 0. The isobaric manifold I Y ( Π, μ) (solid black curve) represents the level set {y|ϕ * (y) = Π}, which divides the space Y into the sublevel set (lower left) and the superlevel set (upper right).The equilibrium manifold M Y EQ (μ) (red line) is located in the superlevel set.The black vector is the normal vector ∂ϕ * (y) of the level set.Since the range of ∂ϕ * (y) = x covers the positive orthant, y ∈ I Y ( Π, μ) exists such that it satisfies U ∂ϕ * (y) = 0, shown by the black circle.This point is on the stoichiometric manifold M Y STO (L = 0).The orange vector represents (y P − y), and the dashed line expresses the tangent plane of the level set at the point y.One can find y ∈ M Y STO (L = 0) such that any y P ∈ M Y EQ (μ) is located on the same side of the tangent plane with the black vector ∂ϕ * (y).For such y ∈ M Y STO (L = 0), the inner product ∂ϕ * (y)(y P − y) is positive.

FIG. 9 .
FIG. 9.Illustration of the proof that K Y (y P ; y) = ∂ϕ * (y)(y P − y) is negative for y ∈ M Y STO (L = 0) when ϕ * (y min )− Π < 0. The isobaric manifold I Y ( Π, μ) (solid black curve) represents the level set {y|ϕ * (y) = Π}, which divides the space Y into the sublevel set (lower left) and the superlevel set (upper right).The equilibrium manifold M Y EQ (μ) (red line) is located both in the sublevel and the superlevel sets.The black vector is the normal vector ∂ϕ * (y) of the level set, which satisfies U ∂ϕ * (y) = 0. Thus, the point y shown by the black circle is on the stoichiometric manifold M Y STO (L = 0).The orange vector represents (y P − y) and the dashed line expresses the tangent plane of the level set at the point y.Any y P ∈ M Y EQ (μ) is located on the opposite side of the tangent plane to the orientation of the black vector.Thus, the inner product ∂ϕ * (y)(y P − y) is negative.