Motivating exampleOur objective is to find a connection between the rate of production of an autocatalytic process, its distance from equilibrium, and the stoichiometry of the underlying reaction network. Consider, for illustration, the following reactive system:$$\begin{array}{rcl}&&{\mathsf{F}}+\alpha \,{\mathsf{A}}{\leftrightarrows}_{{k}_{-1}}^{{k}_{+1}}{\mathsf{B}}\\ &&{\mathsf{B}}{\leftrightarrows}_{{k}_{-2}}^{{k}_{+2}}\beta \,{\mathsf{A}}+{\mathsf{W}},\end{array}$$
(1)
with β > α > 0. The topology of (1) is encoded in its stoichiometric matrix,$$\begin{array}{rcl} \nabla \!\!\!\!&=&\!\!\!\! \begin{array}{c}{\mathsf{F}}\\ {\mathsf{W}}\\ {\mathsf{A}}\\ {\mathsf{B}}\end{array} \left( \begin{array}{cc}-1 & 0 \\ 0 & 1 \\ -\alpha & \beta \\ 1 & -1 \end{array} \right).\\ &&\,\,\begin{array}{rcl} 1 \!\quad& 2 \end{array} \end{array}$$
(2)
The species \({\mathsf{F}}\) and \({\mathsf{W}}\) act as fuel and waste for the overall production of the other species. In what follows, we will treat their concentrations, f and w, as constants. Only a submatrix of (2) is required to capture the autocatalytic behavior of (1):$$\begin{array}{rcl}{\mathbb{S}}\!\!\!\!&=&\!\!\!\!\begin{array}{c}{\mathsf{A}}\\ {\mathsf{B}}\end{array}\left(\begin{array}{cc}-\alpha &\beta \\ 1&-1\end{array}\right).\\ &&\,\,\begin{array}{rcl} \!\!1 \quad & \!\!2 \end{array} \end{array}$$
(3)
This submatrix establishes a connection between the time derivative of the concentrations of the autocatalytic species a and b with the reaction fluxes,$${{{{{{{{\rm{d}}}}}}}}}_{t}\,\left(\begin{array}{c}a\\ b\end{array}\right)={\mathbb{S}}\cdot {{{{{{{\bf{j}}}}}}}},$$
(4)
where j is a vector containing the fluxes of reactions 1 and 2, \({{{{{{{\bf{j}}}}}}}}={({j}_{1}, \, {j}_{2})}^{\top }\). For ideal isothermal solutions, each of these reaction fluxes can be decomposed as the difference of two one-way fluxes obeying the law of mass action:$${j}_{+1}={k}_{+1}\,f\,{a}^{\alpha },\qquad {j}_{-1}={k}_{-1}\,b,$$
(5)
$${j}_{+2}={k}_{+2}\,b,\qquad {j}_{-2}={k}_{-2}\,{a}^{\beta }\,w;$$
(6)
and the components of j are ji = j+i − j−i. Nothing prevents this reactive system from reaching the equilibrium state, since the concentrations of both \({\mathsf{A}}\) and \({\mathsf{B}}\) are unconstrained. To drive the network (1) away from its equilibrium state, we introduce a control mechanism which maintains the concentrations a or b constant thanks to an outgoing flux.When the concentration of species \({\mathsf{A}}\) is controlled (chemostatted), the dynamics of the system is entirely controlled by the evolution law for b, dt b = j1 − j2. The steady state is such that j1 = j2, which shows that there is a tight coupling26 between the two reactions. Here, j1 and j2 are also equal to the production rate \({{{{{{{\mathcal{J}}}}}}}}\) of the overall reaction:$${\mathsf{F}}+\alpha \,{\mathsf{A}} \longrightarrow \beta \,{\mathsf{A}}+{\mathsf{W}}.$$
(7)
This reaction has an overall affinity (which is equal to the opposite of the Gibbs free energy ΔG) \({{{{{{{\mathcal{A}}}}}}}}={\mu }_{{\mathsf{F}}}-{\mu }_{{\mathsf{W}}}-(\beta -\alpha )\,{\mu }_{{\mathsf{A}}}\), where μi is the chemical potential of species i. Since \({\mu }_{{\mathsf{F}}}\) and \({\mu }_{{\mathsf{W}}}\) are fixed, fixing \({\mu }_{{\mathsf{A}}}\) is essential to maintain the system in a non-equilibrium state where \({{{{{{{\mathcal{A}}}}}}}}\) is non-zero. Otherwise, the system will reach equilibrium, where the affinity vanishes. Solving for the steady-state concentration provides the expression for the macroscopic flux of production of species \({\mathsf{A}}\):$${{{{{{{\mathcal{J}}}}}}}}=\frac{{k}_{+1}\,{k}_{+2}}{{k}_{-1}+{k}_{+2}}\left[{a}^{\alpha }\,f-\left(\frac{{k}_{-1}{k}_{-2}}{{k}_{+1}{k}_{+2}}\right)\,{a}^{\beta }\,w\right].$$
(8)
Introducing the reaction quotient Q = a(β−α) w/f and the equilibrium constant K = k+1k+2/k−1k−2 of the global reaction (7), it is easy to show that the global flux \({{{{{{{\mathcal{J}}}}}}}}\) has the same sign as the affinity, since \({{{{{{{\mathcal{A}}}}}}}}=\ln \left(K/Q\right)\) (we work with units where R T = 1). Furthermore, the current goes through a maximum as a function of a whenever$$\frac{{a}^{(\beta -\alpha )}\,w}{f}=\frac{\alpha }{\beta }\,\frac{{k}_{+1}\,{k}_{+2}}{{k}_{-1}\,{k}_{-2}},$$
(9)
which corresponds to the condition Q = Q* = α/β × K. Thus, the maximum rate is reached when the chemical affinity becomes$${{{{{{{\mathcal{{A}}}}}}}^{* }}}=\ln \frac{\beta }{\alpha }.$$
(10)
Hence, the distance from equilibrium at which the autocatalytic network achieves its optimal production rate, namely \({{{{{{{\mathcal{{A}}}}}}}^{* }}}\), is not fixed by the values of kinetic constants or by the equilibrium constant of the global process. It only depends on the stoichiometry of the overall reaction. We illustrate these results in Fig. 1, for the case α = 1, β = 2. Note that this condition on the affinity can also be expressed in terms of the concentration a. Since \(\exp ({{{{{{{{\mathcal{A}}}}}}}}}^{* })=K/{Q}^{* }\), the point of maximum flux is given by$${\left(\frac{{a}^{* }}{{a}_{{{{{{{{\rm{eq}}}}}}}}}}\right)}^{\beta -\alpha }=\frac{\alpha }{\beta },$$
(11)
which means that optimality is reached when the concentration of the chemostat is at half its equilibrium value.Fig. 1: Selected properties of the autocatalytic system (1) for α = 1 and β = 2.a Representation of a reactor in which the concentrations of \({\mathsf{F}}\), \({\mathsf{W}}\), and \({\mathsf{A}}\) are maintained constant thanks to exchanges with external reservoirs (chemostats). b The corresponding global flux \({{{{{{{\mathcal{J}}}}}}}}\) as a function of Q vanishes when Q = K (on the black circle) and reaches a maximum at the same location Q* = K/2 (red triangles). c Species \({\mathsf{F}}\), \({\mathsf{W}}\) are still chemostatted, but now \({\mathsf{B}}\) is chemostatted instead of \({\mathsf{A}}\). d Now the maxima of the corresponding global reaction rate are reached at different locations such that 1/4 ≤ Q*/K ≤ 1/2 (shaded gray area). Note that all the fluxes still vanish when Q = K. e Cloud of points for the exponential of the global affinity \({{{{{{{{\mathcal{A}}}}}}}}}^{* }\), for randomly chosen sets of kinetic rate constants indexed by N. When species \({\mathsf{A}}\) is controlled (red triangles), the optimal affinity always equals \(\ln 2\), as explained in the main text. When species \({\mathsf{B}}\) is controlled (blue squares), the optimal affinity admits both a lower and upper bound, which corresponds to the gray area in the lower plot of b.We can carry a similar analysis if, instead of A, B is the chemostatted autocatalytic species. Now, the dynamics is ruled by dt a = β j2 − α j1, and the steady state is such that \({j}_{1}/\beta ={j}_{2}/\alpha ={{{{{{{\mathcal{J}}}}}}}}\), which is the overall production rate of species B,$$\beta \,{\mathsf{F}}+\alpha \,{\mathsf{B}} \longrightarrow \beta \,{\mathsf{B}}+\alpha \,{\mathsf{W}}.$$
(12)
The steady-state solution now involves polynomials of different order, making it impossible to find explicitly the conditions maximizing the rate of production with the previous method. Nonetheless, we can determine \({{{{{{{{\mathcal{A}}}}}}}}}^{* }\) by numerically finding the value of Q* that maximizes the global reaction rate for randomly generated values of the various kinetic constants. We find that the optimal affinity is now bounded from below (see Supplementary Note 12):$${{{{{{{{\mathcal{A}}}}}}}}}^{* }\ge \ln \frac{\beta }{\alpha }.$$
(13)
Here too, the constraint acting on the optimal distance from equilibrium solely depends on the stoichiometry of the overall reaction. This thermodynamic constraint translates into a threshold for the value of the concentration of the controlled species, b* as$${\left(\frac{{b}^{* }}{{b}_{{{{{{{{\rm{eq}}}}}}}}}}\right)}^{\beta -\alpha }\le \frac{\alpha }{\beta },$$
(14)
where beq is the equilibrium concentration of \({\mathsf{B}}\). Note the importance of the condition β > α > 0, which guarantees the invertibility of the stoichiometric matrix \({\mathbb{S}}\) and leads to equalities or inequalities Eqs. (10–14).Figure 1 summarizes the behavior of the network Eq. (1) in the particular case α = 1 and β = 2. It also reveals that, on top of the upper bound Eq. (14), the location of the global flux is also lower-bounded, Q*/K ≥ 1/4, yielding an upper bound for the affinity \({{{{{{{{\mathcal{A}}}}}}}}}^{* }\). This behavior can be recovered analytically for arbitrary values of α and β (see Supplementary Notes 11, 12).General approach for autocatalytic chemical reaction networks (CRNs)In what follows, we present a general approach that enables the computation of the overall affinity \({{{{{{{{\mathcal{A}}}}}}}}}^{* }\) corresponding to a condition of the local extremum of the global production rate. Our method relies on the stoichiometry of the reaction network, and circumvents the need to explicitly evaluate the steady states or reaction rates of the system under consideration.Following a recent stoichiometric characterization of autocatalysis21, a general autocatalytic network should contain one or several autocatalytic cores. An essential feature of these cores is that they are described by an invertible stoichiometric submatrix. The existence of this invertible submatrix is a sufficient condition for autocatalysis, which implies the absence of mass-like conservation laws13,27,28 thanks to Gordan’s theorem21. In the following, we assume that the full stoichiometric matrix ∇ contains a square submatrix \({\mathbb{S}}\) that is invertible:
(15)
Reactions of \({\mathbb{S}}\) will be called autocatalytic reactions, species of \({\mathbb{S}}\) will be called autocatalytic species, while other species will be called external. Note that we do not require that the matrix \({\mathbb{S}}\) corresponds necessarily to one of the autocatalytic cores listed in Table 1 of the Methods Section. Further, we also allow the network to contain catalytic reactions, (i.e., reactions where the same species can be both reactants and products) in contrast with the assumptions of ref. 21. The presence of a block matrix of zero in the lower right part of the matrix ∇ in Eq. (15) means that there is no boundary flux to or from external species. Such a splitting of species into internal and external ones is a standard assumption in metabolic network analysis29.Table 1 Bounds on the chemical affinities at the maximum of the macroscopic flux for the five autocatalytic cores21Similar to the example above, the concentrations of all external species are assumed to be constant, either because these species are in excess, which is typically the case for food or fuel species, or because they are in contact with a reservoir (chemostatted). Given the structure of the full stoichiometric matrix of Eq. (37), one can show that the network can satisfy detailed balance (see Eq. (38) and Supplementary Note 1). To break detailed balance, in addition, the concentration of at least one of the autocatalytic species should be controlled, by actively maintaining its concentration constant with an outgoing flux.We call this special species the X species, and we call the remaining non-controlled autocatalytic ones the Y species. The stoichiometric matrix \({\mathbb{S}}\) splits into a row vector SX and a matrix \({{\mathbb{S}}}^{Y}\)27:$${\mathbb{S}}=\left(\begin{array}{c}{{{{{{{{\bf{S}}}}}}}}}^{X}\\ {{\mathbb{S}}}^{Y}\end{array}\right).$$
(16)
The corresponding kinetic equations are given by$${{{{{{{{\rm{d}}}}}}}}}_{t}\,x={{{{{{{{\bf{S}}}}}}}}}^{X}\cdot {{{{{{{\bf{j}}}}}}}}+I=0$$
(17)
$${{{{{{{{\rm{d}}}}}}}}}_{t}\,{{{{{{{\bf{y}}}}}}}}={{\mathbb{S}}}^{Y}\cdot {{{{{{{\bf{j}}}}}}}},$$
(18)
where x denotes the concentration of species X and y is the concentration vector of all the Y species. The vector j contains the rates of the autocatalytic reactions, and I is a scalar function describing the exchange of matter with the chemostat.The inverse \({{\mathbb{S}}}^{-1}\) plays an important role:$${{\mathbb{S}}}^{-1}={\left\{{{{{{{{{\bf{g}}}}}}}}}_{\sigma }\right\}}_{{Z}_{\sigma }\in {{{{{{{\mathcal{Z}}}}}}}}},$$
(19)
where gσ is the column of \({{\mathbb{S}}}^{-1}\) associated to species Zσ, which denotes autocatalytic species of the X or of the Y type, and \({{{{{{{\mathcal{Z}}}}}}}}\) is the set of all the autocatalytic species. It follows from the property of the inverse that gσ represents a reaction pathway that produces a single unit of species Zσ without affecting the other species (see Supplementary Note 2). For this reason, gσ is the elementary mode of production of species Zσ21. From all the elementary modes of production, one can build a reaction vector g = ∑σgσ, which represents a combination of elementary modes that increases the amount of all the autocatalytic species by one unit:$${\mathbb{S}}\cdot {{{{{{{\bf{g}}}}}}}}={{{{{{{\bf{1}}}}}}}},$$
(20)
where 1 is a column vector full of ones. The existence of this vector is a sufficient condition for autocatalysis.The steady-state of Eq. (18) is a vector belonging to the right nullspace of \({{\mathbb{S}}}^{Y}\)27,28, which is spanned by gX, the elementary mode of production associated with species X. Thus, the stationary flux vector can be written as$${{{{{{{\bf{j}}}}}}}}={{{{{{{\mathcal{J}}}}}}}}\,{{{{{{{{\bf{g}}}}}}}}}_{X},$$
(21)
which we call a tight coupling condition, since all elementary fluxes are proportional to the global production rate \({{{{{{{\mathcal{J}}}}}}}}\) with constant coefficients of proportionality. Since in general, gX has rational components, it is convenient to rescale this vector to facilitate its interpretation as a mode of production. We show in the Methods section that, after the rescaling, the bounds on the global affinity are also rescaled by the same amount.We now introduce assumptions about the dynamics of the network at the level of elementary reactions. Every such reaction denoted ρ is assumed to be reversible, with a net flux given by$${j}_{\rho }={j}_{+\rho }-{j}_{-\rho }.$$
(22)
Here, one-way fluxes j±ρ obey mass-action kinetics30,31:$${j}_{\pm \rho }={k}_{\pm \rho }\,\prod\limits_{\sigma }\,{z}_{\sigma }^{{S}_{\pm \rho }^{\sigma }},$$
(23)
where zσ is the concentration of species Zσ, and \({{\mathbb{S}}}_{+}\) (resp. \({{\mathbb{S}}}_{-}\)) is the stoichiometric matrix associated to forward (resp. reverse) reactions, such that \({\mathbb{S}}={{\mathbb{S}}}_{-}-{{\mathbb{S}}}_{+}\). Note that the constant concentrations of the fuel/food species have been absorbed in the effective rate constants k±ρ. We choose to work with ideal solutions for clarity of presentation, but we show in Supplementary Note 4 that the results presented below remain valid for non-ideal solutions.Overall affinityThe affinity of an elementary reaction is connected to fluxes by means of the flux–force relationship: \({{{{{{{{\mathcal{A}}}}}}}}}_{\rho }=\ln ({j}_{+\rho }/{j}_{-\rho })\). Taking the linear combination of the elementary affinities, we obtain the global affinity$${{{{{{{\mathcal{A}}}}}}}}=\sum\limits_{\rho }{g}_{X}^{\rho }\,{{{{{{{{\mathcal{A}}}}}}}}}_{\rho }=\sum\limits_{\rho }{g}_{X}^{\rho }\,\ln \left(\frac{{j}_{+\rho }}{{j}_{-\rho }}\right)=\ln \left(\frac{K}{Q}\right),$$
(24)
where K and Q are, respectively, the equilibrium constant and the reaction quotient Q32. Here, the reaction quotient coincides with the concentration of the controlled autocatalytic species (Q = x) because the concentrations of food species have been incorporated in the rate constants, and the stationary global flux is a function of this quantity, \({{{{{{{\mathcal{J}}}}}}}}={{{{{{{\mathcal{J}}}}}}}}(Q)\).We expect that in the cases of interest, this flux will present at least one maximum when the CRN is brought out of equilibrium. Indeed, due to the tight coupling condition, the total entropy production rate (EPR) in the steady-state has a simple expression13,27,28:$$\Sigma =\sum\limits_{\rho }{j}_{\rho }\,{{{{{{{{\mathcal{A}}}}}}}}}_{\rho }={{{{{{{\mathcal{J}}}}}}}}{{{{{{{\mathcal{A}}}}}}}}.$$
(25)
The second law of thermodynamics dictates that Σ ≥ 0, and equality is achieved at equilibrium, where Q = K and both \({{{{{{{\mathcal{J}}}}}}}}\) and \({{{{{{{\mathcal{A}}}}}}}}\) vanish. It follows, from Eq. (24), that \({{{{{{{\mathcal{J}}}}}}}}(Q)\ge 0\) when Q ∈ [0, K]. Let us consider a class of autocatalytic networks, for which species X has a nucleating role24. In that case, the flux \({{{{{{{\mathcal{J}}}}}}}}(0)\) is zero when this species is absent (in other words, for Q = 0). This was the case, for example, with the simple autocatalytic system discussed earlier (see Fig. 1). For these networks, there must be at least one local extremum of the global flux in the interval Q ∈ [0, K]. This regime is optimal, in the sense that species X is produced at a maximal rate. Importantly, the extremum of the global flux and that of the EPR do not coincide because \({{{{{{{\mathcal{A}}}}}}}}\) in Eq. (24) is a non-linear function of Q which also enters in Eq. (25). This is in agreement with the observation that the optimal regime of operation of non-equilibrium systems generally does not correspond to a point where the EPR is extremum33.Response coefficients at maximum fluxWe call Q* the value of Q which makes the global flux maximal, with a zero derivative: \({{{{{{{{\rm{d}}}}}}}}}_{Q}\,{{{{{{{\mathcal{J}}}}}}}}=0\). Due to the tight coupling condition Eq. (21), all reaction fluxes are also at an extremum at Q*: dQ jρ = dQ j+ρ − dQ j−ρ = 0, for any reaction ρ. To characterize this configuration, we introduce the log-derivative of the stationary elementary fluxes, which we call F±ρ :$${F}_{\pm \rho }={{{{{{{{\rm{d}}}}}}}}}_{Q}\ln {j}_{\pm \rho }=\sum\limits_{\sigma }{S}_{\pm \rho }^{\sigma }\,{{{{{{{{\rm{d}}}}}}}}}_{Q}\ln {z}_{\sigma }.$$
(26)
These are the response coefficients of the steady unidirectional fluxes with respect to a change in Q. Because the F±ρs are log-derivatives, all the factors entering rate laws that do not depend on Q will not contribute to these coefficients. This includes the rate constants, which do not appear explicitly. Crucially, the coefficients F±ρs satisfy structural constraints related to the topology of the network. A graphical illustration of these structural relations is provided in Fig. 2 for the particular case of a linear and a branched reaction pathway. In the linear pathway, an arbitrary species Zi is transformed into a product species Zi+1 by a reversible and unimolecular reaction ρi, and then Zi+1 undergoes a similar reaction ρi+1. In such case, both \({j}_{-{\rho }_{i}}\) and \({j}_{+{\rho }_{i+1}}\) depend solely on the concentration zi+1 of species Zi+1. Consequently, \({F}_{-{\rho }_{i}}={F}_{+{\rho }_{i+1}}\), because both terms are equal to \({{{{{{{{\rm{d}}}}}}}}}_{Q}\ln {z}_{i+1}\). For the branched pathway, a species Z0 splits through the reaction + ρ0 into several products Zi with multiplicity \({S}_{-{\rho }_{0}}^{i}\), and Eq. (26) leads to \({F}_{-{\rho }_{0}}={\sum }_{i}{S}_{-{\rho }_{0}}^{i}\,{F}_{+{\rho }_{i}}\).Fig. 2: Structural constraints acting on the log-derivatives of the steady unidirectional fluxes with mass-action kinetics.a When a species Zi+1 is connected by unimolecular pathways upstream and downstream, the two outgoing one-way fluxes have equal log-derivative, \({F}_{-{\rho }_{i}}={F}_{+{\rho }_{i+1}}\). b In a branched pathway, a species Z0 is connected to various species Zi, and the sum of the log-derivative of the forward flux of all products balances the log-derivative of the reverse branched reaction, \({F}_{-{\rho }_{0}}={\sum }_{i}{S}_{-{\rho }_{0}}^{i}\,{F}_{+{\rho }_{i}}\).For a general network, these structural relations take a form analogous to Kirchoff’s laws (see Supplementary Note 3):$${{{{{{{{\bf{F}}}}}}}}}_{-}={{{{{{{{\bf{F}}}}}}}}}_{+}\cdot \left({{\mathbb{S}}}_{+}^{-1}\cdot {{\mathbb{S}}}_{-}\right),$$
(27)
provided \({{\mathbb{S}}}_{+}\) is invertible, which is the case in most networks of interest. Note that the definition of the coefficients F±ρ and the structural constraints Eq. (27) are valid even when Q ≠ Q*. The structural constraints acting on the F±ρs play a key role in our framework, because these coefficients are intimately related to the affinities at the optimal current. Indeed, when Q = Q* one has:$${j}_{+\rho }\,{F}_{+\rho }={j}_{-\rho }\,{F}_{-\rho },$$
(28)
because dQ jρ = 0. Using the definition of \({{{{{{{{\mathcal{A}}}}}}}}}_{\rho }\) and Eq. (28), we obtain:$${e}^{{{{{{{{{\mathcal{A}}}}}}}}}_{\rho }}=\frac{{F}_{-\rho }}{{F}_{+\rho }}.$$
(29)
At this point of optimal current, Eq. (27) defines a linear system of the form \({\mathbb{M}}\cdot {{{{{{{{\bf{F}}}}}}}}}_{{{{{{{{\boldsymbol{+}}}}}}}}}=0\), because the elements of F- can be expressed in terms of those of F+ by using the local affinities. The solutions of this system will be trivial if \({\mathbb{M}}\), which contains information on both the topology and the affinities, is not singular. Enforcing that the determinant of this matrix is zero results in a constraint involving solely the values of the optimal affinities and the elements of the stoichiometric matrices (see Supplementary Note 5). It does not require an explicit evaluation of the steady-state rates or concentrations, nor does it involve kinetic parameters. We illustrate this approach in the next section.The condition on the determinant is useful, but cannot easily be applied to large CRNs. However, using Eq. (29), the affinity of the overall reaction can be expressed only in terms of the F±ρ at the optimal flux:$${{{{{{{{\mathcal{A}}}}}}}}}^{* }=\sum\limits_{\rho }{g}_{X}^{\rho }\,\ln \left(\frac{{F}_{-\rho }}{{F}_{+\rho }}\right).$$
(30)
Finding a bound (a minimum or a maximum) for this global affinity corresponds to solving an optimization problem with Eq. (27) acting as linear constraints and additional constraints due to tight coupling and the second law as detailed in the Methods section below (see Eq. (39)). More details are provided in the Supplementary Note 11, where we establish conditions under which this optimization problem becomes concave. The solution to this optimization problem provides the thermodynamic bounds as well as information on the response coefficients using only the topology of the reaction network. As was the case with the method based on determinants, this link does not rely on an explicit evaluation of the steady-state fluxes or concentrations, nor does it require knowledge of kinetic parameters. In addition to this, the optimization approach can easily be used with large CRNs.Since the bounds on the optimal affinity do not depend explicitly on the expressions of the steady-state concentrations or reaction rates, our method is applicable to large and complex networks, in which these concentrations and rates are too complex to be computed. For the same reason, the bounds hold even if the system features multistability, which is often found in autocatalytic networks1. We show this explicitly for the bistable Schlögl model34 in Supplementary Note 15.As an illustration of our approach, we derive lower bounds satisfied by the global affinity of the Hinshelwood autocatalytic cycle35 with intermediate species, as represented in Fig. 3a. When the intermediate species, B or D, are controlled, the constraints on the F±ρs can be satisfied and the global flux has a zero-derivative maximum (see the blue curve in Fig. 3b). In that case, the affinity at the optimum of the current is lower-bounded,$${{{{{{{{\mathcal{A}}}}}}}}}^{* }\ge \ln 4,$$
(31)
and the bound can be approached as closely as desired with an appropriate choice of rate constants, as shown in Fig. 3c (blue triangles). When, instead, the concentration of A or C is maintained constant, the constraints acting on the F±ρs cannot be satisfied (see Supplementary Note 14). As a consequence, the global flux, \({{{{{{{\mathcal{J}}}}}}}}\), has no zero-derivative maximum and, then, is a decreasing function of Q whose largest value occurs at Q = 0 (see red curve in Fig. 3b), where the definition of the F±ρs and Eq. (28) do not apply. As a result, \(\exp ({{{{{{{{\mathcal{A}}}}}}}}}^{* })\) diverges. In addition to the Hinshelwood cycle, we also analyzed in Supplementary Note 13 all the autocatalytic cores21; the results are summarized in Table 1. Importantly, the bound found for the five autocatalytic cores remains valid if unimolecular segments of reactions are added to any of these networks (see Supplementary Note 10 and Supplementary Fig. 1).Fig. 3: Optimal properties of the Hinshelwood cycle.a Schematic representation of the Hinshelwood cycle with two types of species: red species (A or C) and blue species (B or D). b Global flux normalized by its maximum value \({{{{{{{{\mathcal{J}}}}}}}}}^{* }\) as a function of Q/K when a specific species of the cycle is chemostatted. This flux has a zero-derivative maximum (blue triangle) if a blue species is chemostatted or simply reaches its maximal value at Q = 0 (red square) if a red species is chemostatted. c Cloud of points for the exponential of the global affinity \({{{{{{{{\mathcal{A}}}}}}}}}^{* }\), for randomly chosen sets of kinetic rate constants indexed by N. If a blue species is chemostatted, points are lower-bounded by 4 (blue triangles) otherwise, the affinity diverges (red squares). d Global flux normalized by its maximum value \({{{{{{{{\mathcal{J}}}}}}}}}^{* }\) as a function of the degradation rate κ of a specific species. When a blue species is degraded, a zero-derivative maximum exists and is reached at a finite value of κ, while if a red species is degraded, the global flux is monotonously increasing and reaches its maximal value at infinity. Simulation parameters for b and d are k+1 = k+3 = 1 and k+2 = k+4 = k−1 = k−2 = k−3 = k−4 = 0.1.Bounds for type I and type II networksNow, two situations arise depending on whether species X appears or not as a reactant in the overall reaction. When species X is a reactant in the overall equation, the reaction vector gX defines a seed-dependent mode of production23,24. In that case, the overall reaction is$$\alpha \,{\mathsf{X}}+ (\diamond) \longrightarrow \beta \, {\mathsf{X}}+(\diamond),$$
(32)
with α and β being integers such that β > α > 0, and (♢) represents all the spectator species. Note that in our formalism, the concentrations of the external species are absorbed in the rate constants, so these species do not appear explicitly in the overall reaction. The simplest network obeying Eq. (32) is a generalized version of the Type I core presented in example (1) with an arbitrary number of intermediate species. As shown in Supplementary Note 8, the global production rate in a generalized Type I attains its maximum when$${{{{{{{{\mathcal{A}}}}}}}}}^{* }=\ln \left(\frac{\beta }{\alpha }\right).$$
(33)
Let us now consider a specific class of networks that we call non-intersecting, for which the stoichiometric matrix of the products, \({{\mathbb{S}}}_{-}\), is lower triangular up to permutations of species and reactions (see Supplementary Note 6). This means that in non-intersecting networks, each reaction only produces downstream species. We have shown that for non-intersecting networks, Eq. (33) represents the lowest achievable bound for \({{{{{{{{\mathcal{A}}}}}}}}}^{* }\) when the overall reaction is seed-dependent. We conjecture that this proposition is not restricted to non-intersecting networks but may hold in general.Conversely, when species X is not present as a reactant in the overall reaction, the reaction vector gX defines a seed-independent mode of production. In that case, the overall reaction reads$$(\diamond) \longrightarrow \beta \,{\mathsf{X}}+(\diamond),$$
(34)
where β is a strictly positive integer. It should be noted that Eq. (34) can describe overall reactions that are not necessarily derived from an autocatalytic network. However, the presence of an optimal flux and the bounds on the corresponding affinity can only be guaranteed if there is an underlying autocatalytic network present. The simplest topology compatible with Eq. (34) is a generalization of the Type II core:$$\alpha \,{\mathsf{X}}\rightleftarrows {{\mathsf{Y}}}_{1}\rightleftarrows \cdots \rightleftarrows {{\mathsf{Y}}}_{n}\rightleftarrows {{\mathsf{Y}}}_{1}+\beta \,{\mathsf{X}}.$$
(35)
This network is also non-intersecting (see Supplementary Note 9) and the global affinity associated with the maximum of \({{{{{{{\mathcal{J}}}}}}}}\) is$${{{{{{{{\mathcal{A}}}}}}}}}^{* }=\ln \left(\frac{\beta }{\alpha }+1\right).$$
(36)
Structural constraints limit the regime of optimal flux in autocatalytic reaction networks
