Origin of exponential growth in nonlinear reaction networks

Significance Natural systems (e.g., cells and ecosystems) generally consist of reaction networks (e.g., metabolic networks or food webs) with nonlinear flux functions (e.g., Michaelis–Menten kinetics and density-dependent selection). Despite their complex nonlinearities, these systems often exhibit simple exponential growth in the long term. How exponential growth emerges from nonlinear networks remains elusive. Our work demonstrates mathematically how two principles, multivariate scalability of flux functions and ergodicity of the rescaled system, guarantee a well-defined growth rate. By connecting ergodic theory, a powerful branch of mathematics, to the study of growth in biology, our theoretical framework can recapitulate various growth modalities (from balanced growth to periodic, quasi-periodic, or even chaotic behaviors), greatly expanding the types of growing systems that can be studied.

Supplementary Information for 7 8 Origin of exponential growth in nonlinear reaction networks 9 Wei-Hsiang Lin * , Edo Kussell, Lai-Sang Young and Christine Jacobs-Wagner * 10    Tables   54  76  Table S1: Common flux functions 54 77 Table S2: Reactions and flux functions in the autocatalytic biosynthesis model 55 78 Our goal is to characterize a mathematical class of flux networks that give rise to exponential 92 dynamics in the long term. While a system evolves and increases its size, the relative proportion 93 of each node also evolves. Since a growing process is not stationary, acquiring its long-term 94 statistical averaging is generally not possible. Here, we explored the class of scalable functions 95 (described in section 2) for constructing dynamical systems. The proportional scaling property 96 allows us to decouple the n-dimensional dynamics into 1-dimensional "bulk dynamics" and (n-1)-97 dimensional "relative proportion dynamics". The latter dynamics is confined in the (n-1)-98 dimensional unit simplex, which enables us to apply ergodic theory for long-term statistical 99

Methods: Simulation and analysis of scalable reaction networks
averaging. 100

101
In section 1, we introduce the basic mathematical structure for flux networks and the concept of 102 long-term growth rate. In section 2, we introduce scalable reaction networks and present the 103 main result (Theorem 2.3). This result guarantees, via the use of ergodic theory, converged long-104 term growth rates for many (though not all) initial conditions of scalable networks. Examples 105 are discussed. In section 3, we consider the addition of scalable noise and show that the 106 presence of such noise, together with a regenerative condition on the networks, guarantees 107 well-defined long-term growth rates for all initial conditions. Along with the theory, we 108 investigate concrete examples to illustrate non-steady-state growth modes and rich behaviors 109 (including chaotic dynamics) for growing reaction networks. In section 5, we generalize the class 110 to asymptotic scalable networks, which replaces the stringent "proportional scaling" condition to 111 a more relaxed "asymptotic scaling". In section 6, we study the necessary structural motifs for 112 the network to have a positive long-term growth rate. 113 114 1. Flux network ( , , ) and long-term growth rate 115 As in thermodynamics formulations, we consider a system that is connected to an environment. 116 The environment (denoted as ) acts as a reservoir and provides unlimited supply and removal 117 of materials. The system is composed by multiple nodes 0 , . . . , 2 and interconversion reactions 118 0 , . . . , 4 . A reaction can happen between the environment and system nodes, as well as 119 among multiple system nodes. In general, a reaction 5 can be represented by 120 121 5 : 9 + 0 0 + ⋯ + 2 2 → 9 + 0 0 + ⋯ + 2 2 (1) 122 where ? and ? ( = 0, . . . , ) represent the stoichiometry coefficients in this conversion. 123 Specifically, whenever the reaction happens, a proportional amount of , 0 , . . . , 2 is converted 124 into another proportional amount of , 0 , . . . , 2 , where both proportions follow the relative 125 ratio of D and D . It is clear that if ? > ? , the reaction 5 consumes ? and decreases the 126 amount of material on node ? . 127 128 The above description is simplified by classifying nodes that are increased, decreased, or 129 unchanged during the conversion. If we denote ? ≡ ? − ? and environment ≡ 9 , then the 130 reaction 5 can be rewritten as (2) 132 for = 0,1, . . , . Nodes ? with ? < 0 (or ? > 0) are referred to be at the upstream (or 133 downstream) of the reaction 5 . Conceptually, each conversion reaction can be regarded as an 134 "action" that decreases the amount of its upstream nodes and increases the amount of its 135 downstream nodes. This can be visualized as the network diagram: 136 137 where { S } and { T } are the upstream and downstream nodes of the reaction , respectively. 138 Note that the network diagram is more complex than the directed graph -it allows multiple 139 upstream/downstream nodes. 140

3) 149
The stoichiometry matrix and the network diagram are below: 150 151 152 Above, we used the stoichiometry matrix ?5 to describe the topology of the reaction network. 153 Each reaction 5 is further associated with a flux function 5 ( (4) 169 The flux magnitude 5 is usually a nonlinear function that depends on the amount or relative 170 amount of materials in system nodes. In general, we write 5 ( ) where = ( 0 , . . . , 2 ) p is a 171 column vector. For the general nonlinear function 5 ( ), there is no easy way to explicitly solve 172 the equation. Our approach is to constrain the function form of 5 ( ) and study the long-term 173 dynamics of the system. 174

175
Let ≡ 0 +. . . + 2 denote the total amount of material in the system, or the system size. We 176 are interested in flux networks that have asymptotic exponential growth of ( ). We define 177 long-term growth rate by the following limit (if it exists): 178 ≡ lim l→s 0 l log ( ) (5) 179 Note that for an arbitrary flux network, this limit may not exist. The asymptotic behavior of 180 0 l log ( ) can be different, such as (1) ( ) becomes negative in the long term, and hence is 181 not well-defined, (2) ( ) diverges at a finite time b and hence is not well-defined, and (3) 182 0 l log ( ) diverges to infinity as → ∞, and (4) 0 l log ( ) oscillates indefinitely and does not 183 converge to a real number as → ∞. Therefore, it is clear that constructing flux networks with 184 arbitrary flux functions may not give a well-defined long-term growth rate. 185 The long-term growth rate can be positive, negative or zero, which correspond to exponential 186 growth, exponential decay, or other sub-exponential dynamics, respectively. The value of may 187 also depend on the initial condition (0) of the system. Each trajectory can have a different . 188 We explicitly denote this by c(9), if necessary. 189 Exponentially growing systems are often studied by starting with a system of equations

203
In this manuscript, we denote 2 : { ∈ ℝ 2 | ? ≥ 0} as the non-negative quadrant and 204 2ƒ : { ∈ ℝ 2 | ? > 0} as the positive quadrant. For all discussions, we will assume that the 205 system ( ) has the initial condition (0) ∈ 2 , since the amount of starting material (e.g., 206 metabolites) should be nonnegative in living systems. Recall that a node ? is upstream of flux 207 5 , if the stoichiometric coefficient ?5 < 0.  SRNs have the nice property that the vector field satisfies the scaling property ? ( ) = ? ( ). 221 Therefore, once we specify the vector field on the − 1-dimensional unit simplex 2'0 ≡ { ∈ 222 2 | 0 +. . . + 2 = 1}, the entire vector field can be obtained from the scaling relation. 223 Furthermore, as we will show below, this scaling property is closely related to the asymptotic 224 exponential growth of the system. 225 To analyze the scalable network, we change variable from the -coordinate to ( , )-coordinate 226 where = 0 +. . . + 2 and ? = ? / for = 1, . . . , (see Fig. S1). The vector = 227 ( 0 , . . . , 2 ) p is the relative fraction of each node in the system. In general, the differential 228 equation Since is a function of , the dynamics of ( ) is determined by the dynamics of ( ). For 238 scalable fluxes, the trajectory ( ) is confined in the − 1-dimensional unit simplex space 239 2'0 ≡ { ∈ 2 | 0 +. . . + 2 = 1}. To see this, note that the upstream-limited condition (C2) 240 guarantees that whenever the upstream node of a flux is zero, the flux itself is also zero. Hence, 241 the non-negative quadrant 2 is forward invariant for ( ) (9) 245 Our next step is to replace the time average by the space average on the space 2'0 . To 246 proceed, we need to introduce some concepts about ergodic theory. Below, we adopt the 247 notation ( ) = ( 0 ( ), . . . 2 ( )) p as a vector function and denote = ( 0 , . . . , 2 ) p as a 248 vector on 2'0 . We denote ℬ( ) as the Borel sets of , and write a probability measure 249 ( ) as . Note 1: Given an initial condition (0) that is -regular, the ergodic probability measure can 273 be constructed as follows. Let « be the characteristic function such that: 274 where ⊆ 2'0 is a Borel set. We define a family of probability measures Conversely, for any initial condition (0), if l converges to an ergodic measure , then (0) is 283 -regular. We note, however, that l may not converge for some (0), and that the limit 284 measure (even when it exists) need not be ergodic and may depend on (0). 285 Note 2: If ¯( 0) is another initial condition with ( ( ),¯( )) → 0 as → ∞, then the growth 286 rate at ¯( 0) is also well-defined and equal to . 287 Note 3: It is possible that for some initial conditions (within a set of -measure zero) the growth 288 rate does not converge to the value given by Eq. 12. Here, we analyze a deterministic example to 289 illustrate this possibility. 290 We consider a deterministic case on the simplex space 2'0 with 291 n = 3 (see figure at right). The simplex space contains two fixed 292 points: 5 is unstable and ± is stable. For this system, we find 293 two ergodic measures on \ : the delta measures ( − 5 ) ≡ 5 294 and ( − ± ) ≡ ± . Since ± is a stable attractor, all trajectories 295 with the initial condition (0) ≠ 5 converge asymptotically to 296 the constant trajectory ( ) = ± and therefore have = ( ± ). 297 The trajectory with (0) = 5 has = ( 5 ). Thus, if ( 5 ) ≠ 298 ( ± ) then with respect to the ergodic measure ± we have that { 5 } is a set of measure zero for 299 which does not converge to the value given by Eq. 12; and the same is true, respectively, for the 300 ergodic measure 5 and the set \ \{ 5 }. 301 This discussion, together with Theorem 2.3, shows that for a scalable network, many initial 302 conditions (0) are guaranteed to have well-defined long-term growth rates. This includes all 303 (0) for which the rescaled initial conditions (0) are regular. 304 We now consider a few illustrative examples: 305 (1) For a flux network with ( ) converging to a fixed point * , the ergodic probability 306 measure is a Dirac -measure supported on * . The ergodic averaging for the long-term 307 growth rate is reduced to = ( * ). The long-term dynamics converges to balanced 308 exponential growth on the half line { = * , > 0}, with ( ) → * µ(™ * )l . 309 (2) For a flux network with ( ) converging to a limit cycle , the ergodic probability 310 measure is supported by the limit cycle. The long-term growth rate can be calculated by the 311 contour integral = ∮¸( ) , and the dynamics converges to exponential growth with 312 periodic functions ? ( ) as prefactors, namely, ? ( ) → ? ( )~l. 313 (3) There are attractors of ( ) that do not have an ergodic probability measure.
In the repressilator, the synthetic flux of nodes \ , ] ,^ are repressed by node ] ,^, \ , 336 respectively. We implemented this by defining the synthetic flux function \ , ] ,^ as decreasing 337 functions of \ , ] ,^. With ? denoting the fraction of ? , the flux functions are: 0 ( ) = 338 The autocatalytic repressilator network can exhibit different growth modes, depending on the 340 cooperativity coefficient (see Fig. 3A This system can be mapped to a network diagram by defining 12 scalable flux functions (see Fig.  350 S2a), three for the birth process and nine for the death process. The death process is driven by a 351 nonlinear flux function in the form of ( D ? )/ . 352 In our simulation, we find that the trajectory ( ) of the flux network exhibits similar behavior 353 as in the original May-Leonard system. Consistent with the analytical results of the May-Leonard 354 system, varying \ and ] changes the system behavior drastically (see Fig. S2): 355 (1) For \ + ] < 2 (Fig. S2B), the system has balanced growth and three species co-exist. 356 357 (2) For \ > 1 and ] > 1 (Fig. S2C), species 0 is dominant and grows exponentially, while 358 the other species \ and ] go extinct. The system still undergoes balanced growth, with 359 fixed point * ( ) = (1,0,0) on the boundary of \ . 360 361 (3) For \ + ] > 2 and \ < 1 < ] (Fig. S2D), ( ) asymptotically converges to a 362 heteroclinic cycle on the boundary of the system. The heteroclinic attractor ≡ { ∈ 363 \ | 0 \ ] = 0} is the omega-limit set of ( ), but ( ) is not regular with respect to any 364 ergodic probability measure. In this case, each component ? ( ) oscillates between 0 and 1 365 while each period becomes exponentially longer. On the other hand, the overall system ( ) 366 grows unboundedly without the long-term growth rate converging to a real number (see 367 Remark (3)  In this section, we study the relation between the SRN framework and chemostat-type 388 equations. The general SRN framework assumes unlimited external resources, which conflicts 389 with chemostat conditions in which an essential nutrient is in limited amount. However, for a 390 chemostat equation, there exists an equivalent SRN model that gives rise to the same dynamics. where is the dilution rate of the culture, È2 is the constant level of the external nutrient 402 supply, is the stoichiometry constant for the biomass yield, and ( ) = 45É L OEƒL is the 403 Monod's function, which describes the dependence of the nutrient uptake rate (influx) on the 404 extracellular nutrient level ( ). 405 406 The above chemostat equations are not scalable (e.g. ( ) does not satisfy ( ) = ( )). 407 However, since the equations are formulated in terms of concentrations ( ) and ( ), the 408 system dynamics should be independent of the chemostat volume. Intuitively, scaling the 409 chemostat volume to different sizes should yield identical dynamics of concentrations. 410 Therefore, if we reformulate the equation in terms of the total number of molecules, the system 411 should be an SRN. 412 413 To rigorously show this, we include the solvent molecule species in the equations. Let , , 414 and be the total numbers of nutrient molecules, cells, and solvent molecules in the system, 415 and let ¦ , c , Í be their unit volumes. The total culture volume > 0 is a time-independent 416 positive constant and is given by This imposes a constraint among , , and , and we have the relations = / , = / , 421 and = / . We refer to the quantities ( , , ) as intensive variables and ( , , ) as 422 extensive variables, similar to thermodynamics definitions. 423 424 Since is a constant, we have 425 Note that the variable

435
For the solvent variable ( ), by the constraint of constant total volumes , we have 436 Since both Ï ( ), Ï ( ) are scalable, ′( ) is also scalable. Therefore, the ODE of the new 438 variables ( , , ) is an SRN and is equivalent to the chemostat equations. The same mapping 439 can be applied for studying more complex chemostat dynamics, including multi-species 440 interactions (as shown in Fig. 4E and described in SI in the Method section). 441 442 Remark: For the above equivalent SRN model mentioned above, the long-term growth rate of 443 the system ( , , ) is zero. This is due to the constraint that the total volume = ¦ + 444 c + Í is a positive constant, and hence it is impossible to have > 0 ( diverges to 445 infinity) or < 0 ( converges to 0). 446 447 The SRN framework provides additional insight into how biomass growth rate, ( ), is related 448 to the dilution rate . We define 449 If the dilution rate is higher than the maximal biomass growth rate (i.e., > 45É ), then c < 451 0 and all cells in the chemostat will be washed away. In general, the constant dilution rate is equal to the average biomass growth rate ( / ) 462 with respect to an ergodic measure . For balanced growth (as a special case), we have = 463 ( * / * ) where ( * , * , * ) is the fixed point of the system. The same relation holds for 464 more complex chemostat dynamics. For example, when a multi-species chemostat exhibits 465 oscillatory behavior, the limit cycle has the property that the average biomass growth rate on 466 the limit cycle is equal to the dilution rate. This also applies for dynamics which include a limit 467 torus or chaos. For scalable reaction networks, ergodic theory offers a conceptual understanding of long-term 475 growth rates. Because the rescaled system takes values in a compact region, long-term growth 476 rates are guaranteed to exist for many initial conditions, namely those that are regular with 477 respect to some ergodic probability measure. Given a particular initial condition, however, 478 growth properties are not always easy to determine, because not every initial condition is 479 regular and more complicated systems often have multiple ergodic measures each with their 480 own set of regular points and growth rates. This less-than-ideal state of affairs is improved 481 greatly in the presence of random noise. Real-world systems are inherently noisy and 482 mathematically stochastic terms are often added to model uncontrolled fluctuations in physical 483 and biological systems. As we will show in this section, the averaging effect of random noise 484 leads to a simpler and more tractable dynamical picture. 485 To investigate how flux noise affects the system's behavior, we generalize the ordinary 486 differential equation Here, ? represents the standard Brownian motion, with integral in Stratonovich-sense. ∎ 501 In the above SDE, the deterministic part of the equation, i.e., ( ), has identical requirements 502 as in scalable flux functions. We consider a noise function ? which also scales with the system 503 size. This is motivated by the following reasons: 504 (i) In nature, the environment is never truly static. For example, the environmental 505 temperature has nonzero fluctuations and this globally affects all enzyme activities. Consider 506 a flux = ? where the rate constant = ( ) is a function of temperature . Fluctuations 507 of the environmental temperature create noise ≈ ′( ) on the rate constant . 508 Hence, this creates noise on flux ≈ ( ) ? = ( / ) . This noise scales with the flux 509 magnitude and therefore scales with system size. 510 (ii) Consider a system with finite space, for example, a bacterial culture in a turbidostat, the 511 population is kept at a constant size while the excess of growing bacteria is discarded by 512 medium dilution. The dilution process usually creates a sampling noise ∼ √ . In this case, 513 the growth dynamics is equivalent to an exponentially growing population with / copies 514 of subsystem with size . If the exponential growing population has size and a noise level 515 • √ , then each subsystem has a fluctuation √ . In this case, the equivalent exponentially 516 growing population has a noise level that scales with system size. 517 In addition to scalable noise, our formulation allows a system with > 0 to have other types of 518 noise that grow slower than the system size. These "sub-scalable noises" will decay to zero as 519 the system grows to infinity. Therefore, only the scalable noise term remains. 520 To establish the existence and uniqueness of solutions for all > 0 for the SDE in Eq. Theorem 5.6.1). Details are given in the Appendix. 526 As in the deterministic case, we now consider the rescaled SDE ( ) defined by ? ≡ ? /| |, 527 | | = 0 +. . . + 2 . We begin with a result asserting that starting from an initial condition where 528 ? > 0 for all , i.e., all the substances in the network are present in positive amounts, no 529 substance will be depleted in finite time. Proof: One way to prove this is to produce a Lyapunov function on the interior of 2'0 with 533 ( ) → ∞ as tends to 2'0 , the boundary of 2'0 . The function is to have the property 534 * ≤ for some fixed > 0 where * is the diffusion operator of the SDE; it controls the 535 speed with which solutions can approach the boundary. A complete proof is given in the 536 Appendix. ∎ 537 The SDE (19) generates a time-homogeneous Markov process on 2'0 . For 9 ∈ 2'0 and > 0, 538 we use L ( 9 , • ) to denote the transition probability starting from 9 , that is, if ⊆ 2'0 is a 539 Borel set, then L ( 9 , ) is the probability that ( ) ∈ given that (0) = 9 . A probability 540 measure is said to be stationary if ( ) = ∫ L ( , ) ( ) for all Borel sets and all > 0. 541 A stationary probability measure is called ergodic if there is no subset ⊂ 2'0 with 0 < 542 ( ) < 1 such that L ( , ) = 1 for ∈ and > 0. The Ergodic Theorem in section 2 applies 543 in this setting as it does in the deterministic case. 544 (19), and let be an ergodic probability measure of the rescaled system. Then we have either 546 ( 2'0 ) = 0 or 1. 547 Proof: This follows immediately from Proposition 3.2, which says that the interior of 2'0 is an 548 invariant set, so it must have measure 0 or 1 by the ergodicity of . ∎ 549 Below, we give a sufficient condition for the existence of long-term growth rates for noisy 550 scalable networks. Following the notation in Eq. (20), we call a system regenerative if for 551 every = 1, . . . , , the th component of the drift vector ? ( ) > 0 whenever ? = 0. 552 Mathematically, this means the drift term in the SDE points into the interior of the simplex 2'0 553 at all points on its boundary. 554 Theorem 3.4: Consider an SRN with noise as defined by Eq. (19) and assume that the system is 555 regenerative. Then, there is a number with the property that starting from any initial condition 556 (0) ∈ 2'0 , the long-term growth rate converges almost surely to . sets. Finally, the Birkhoff Ergodic Theorem applied to the function ( ) gives the result that the 567 long-term growth rate of every (0) ∈ 2'0 exists almost surely and is equal to = 568 ∫ ( ) ( ). ∎ 569 We remark that even though the invariant density ( ) is positive on all of 2'0 , it is in general 570 much more concentrated in some regions of 2'0 than in other regions, that is to say, certain 571 network configurations (in terms of proportions of the substances present) are much more 572 prevalent even though there is a small theoretical probability of exploring other parts of the 573 phase space. What makes the situation in Theorem 3.4 more tractable than the noise-free case 574 is the uniqueness of ergodic measure and the fact that every initial condition is regular with 575 respect to the unique ergodic measure. 576 The general case (of not-necessarily-regenerative networks): Not all realistic biological 577 networks are regenerative. Not being regenerative means that there is no explicit mechanism in 578 place for each substance to restore itself as it heads towards relative depletion. To be precise, 579 relative depletion here means the rescaled amounts ? ( ) (but not necessarily the actual 580 amounts ? ( )) become arbitrarily small. Ergodic measures always exist for the same reasons as 581 before, but unlike the regenerative case, there may exist ergodic measures supported on 2'0 . 582 For example, if ? ≡ 0 on { ? = 0}, then starting from (0) with ? (0) = 0, almost surely 583 ? ( ) = 0 for all > 0 (see condition (N2)), so there is a subsystem defined on the (n-2)-584 dimensional simplex 2'0 ∩ { ? = 0} with its own ergodic measures. 585 In this general case, let us consider only initial conditions starting from the interior of 2'0 . 586 There is the following dichotomy: 587 Case 1. Existence of an ergodic measure with a density ( ) > 0 on all of 2'0 . In this case, 588 long-term growth rates of all (0) in the interior of 2'0 converge almost surely to = 589 ∫ ( ) ( ) , as before. 590 Case 2. All ergodic measures are supported on 2'0 . In this case, starting from (0), in the 591 interior of 2'0 , one or more of the substances will come close to relative depletion (though 592 none can deplete in finite time according to Proposition 3.2). More precisely, given any 593 (0), for any > 0 we have Prob{ ( ( ), 2'0 ) > } → 0 as → ∞. We summarize our 594 findings below, leaving mathematical proofs to be published elsewhere. 595 Since trajectories spend nearly 100% of their time near 2'0 , their large-time behaviors are 596 dictated by the ergodic measures supported on 2'0 . Now 2'0 is a finite union of (n-2)-597 dimensional simplices and an analysis similar to that for 2'0 can be made for each of these 598 simplices. After successively reducing dimension, we arrive at the conclusion that i) the system 599 has at most a finite number of ergodic invariant measures, each supported on a lower 600 dimensional simplex and ii) in this lower dimension, the system has a density with respect to the 601 Lebesgue measure. 602 It may happen that ( ) converges almost surely to one of these ergodic measures , in which 603 case its long-term growth rate is given by = ∫ ( ) ( ). Another possibility is for ( ) to 604 explore multiple ergodic measures one after another (in a way similar to the heteroclinic loop 605 example in section 2). In this case, there may not be a single number that correctly describes the 606 trajectory's long-term growth rate. 607 608 4. An exponentially growing system with chaotic dynamics 609 In the following example, we study a flux network converging to a nonperiodic trajectory and 610 investigate the effect of noise on its long-term growth rate. 611 The double repressilator network. We construct a flux network with two repressilators 612 repressing each other (see Fig. 2F and Fig. S3). The double repressilator network consists of 613 seven nodes: 0 represents the common precursor, \ , ] , and ^ form the first repressilator, 614 and ë , ì , and í form the second one. In addition, the two repressilators are coupled by 615 mutual repression, where nodes \ and ë repress the influxes of each other. The full 616 differential equation is: respectively. We define ≡ 5 / ± as the relative repression strength between the two 623 repressilators. By varying , the system exhibits complex behaviors, including periodic, quasi-624 periodic, and non-periodic trajectories of ( ). For different values of parameter , the omega-625 limit set of the ( ) can have complex geometry such as a limit torus ( = 20), a complicated 626 limit cycle ( = 100) and a strange attractor ( = 500) (Fig. 2G). 627 To show how the attractor of ( ) transitions from a limit cycle to a limit torus to a more 628 complex attractor, we project the attractor on the \ − ] plane and get a Poincare section at 629 { \ = 0.2}. For various , we see that the system exhibits rich behavior, including period 630 doubling and transition to chaos (Fig. 2G, insets). For the full bifurcation diagram, see Fig. S3A. 631 To quantitatively analyze the attractor, we numerically calculated the largest Lyapunov 632 exponent (LLE) of the trajectories Fig. S3B). LLE quantifies how small perturbations propagate 633 along the solution trajectory, and is a useful indicator for chaotic dynamics: positive LLE indicates 634 strange attractor while LLE=0 indicates a limit torus, a limit cycle or a fixed point. We found that 635 for with values between 400-500, LLEs are significantly positive, which is consistent with what 636 is qualitatively observed in the bifurcation diagram. There are also sporadic "chaotic islands" for 637 = 120. 638 For comparison, we introduce a linear noise ? ( ) = ℎ ? ? on the SDE. We find that the LLEs 639 with and without noise (black versus light blue lines in Fig. S3B) are comparable, indicating that 640 the chaotic behavior is robust under perturbation. Although it is difficult to show the 641 convergence of in a nonperiodic, deterministic system, there is noise associated with fluxes in 642 most realistic systems. By introducing small noise, we can infer whether the long-term behavior 643 observed in the deterministic system is robust. By Theorem 3.4, converges under arbitrarily 644 small ℎ ? . We simulate long-term growth rate with ℎ ? = 0 and ℎ ? = 0.1 for various (Fig. S3C). 645 For non-chaotic regimes, noise has little or no effect on , while at the border between periodic 646 and chaotic regimes (~420), substantially differs whether noise is present or not. The noise 647 seems to affect the transition between different types of attractors. 648 The chaotic growing system indicates that an exponentially growing system can have a complex 649 trajectory that is unpredictable in the long-term. Still, the long-term growth rate of the system 650 converges. We conclude that while the trajectory is sensitive to small perturbations, the long-651 term growth rate, which is a statistical averaging on the projected system, can be still be robust. 652 653 5. Asymptotically scalable reaction networks 654 The condition (C3) of scalable flux requires that the flux function satisfies ( ) = ( ) for a 655 scalar > 0. This is a relatively strong condition that reduces the degrees of freedom of the 656 function by one dimension. Here, we show that a more general class of flux functions, referred 657 to as asymptotically scalable functions, gives a similar convergence property on long-term 658 growth rate. 659 . We obtain the following proposition, which is analogous to Theorem 2.3: 688

Proposition 5.3: Given an asymptotically scalable reaction network with the initial condition 689
(0) ∈ 2 , consider the long-term dynamics of ( ): For case (iii), we obtain = 0 directly since ( ) is bounded. 702 Finally, there are cases not listed above, which are 0 = liminf l→s ( ) < limsup l→s ( ) and 703 liminf l→s ( ) < limsup l→s ( ) = ∞; in these cases, convergence of is inconclusive. 704 6. Autocatalytic flux networks 705 One of the essential features of living systems is the existence of autocatalysis in suitable 706 environments. Scalable reaction networks with positive growth rates can be used to address a 707 range of questions related to autocatalytic systems. For example, what are the fundamental 708 motifs in an autocatalytic flux network, and can we draw inferences about the long-term growth 709 dynamics of a network by examining its topological structure? 710 In an SRN, one can consider the growth rates of individual nodes or fluxes as [ † ] ≡ 711 lim l→s (1/ )log † or [ 5 ] ≡ lim l→s (1/ )log 5 , when the limit exists. Given an SRN with > 712 0, not every node D in the network would be expected to grow as fast as the system size N. If a 713 node D grows slower than the system, we have D → 0 as → ∞. In this case, we say that the 714 node D does not contribute to the growth of the system. Our goal is to characterize the nodes 715 that contribute to system growth and to study the mutual dependence of growth rate between 716 fluxes and nodes. 717 To describe the structure of a flux network ( , ), we adopt the following terminology: A node 718 ? is upstream of a reaction 5 , if ?5 < 0. A node ? is downstream of a reaction 5 , if ?5 > 719 0. Given a reaction 5 , let ( 5 ), ( 5 ) denote the collection of its upstream and 720 downstream nodes, respectively. A reaction 5 is an influx of node ? , if ? ∈ ( 5 ). A 721 reaction 5 is an efflux of node ? , if ? ∈ ( 5 ). Given a node ? , let ( ? ), ( ? ) denote 722 the collection of its influxes and effluxes, respectively. 723 724

725
The above definition only involves the network connection; it is not related to the flux function 726 set . If we associate a flux function set to a flux network ( , , ), then we can introduce the 727 following concept of "maintenance set": 728 Definition 6.1: Let = { 0 , . . . , 2 } denote nodes in the system. We say a node ? maintains the 729 reaction 5 , if ? = 0 implies 5 ( ) = 0. The collection of nodes that maintains 5 is called the 730 maintenance set of reaction 5 , denoted by ( 5 ). ∎ 731 Remark: One reaction can have zero, one, or multiple maintenance nodes. For a scalable 732 reaction flux, all of its upstream nodes are in the maintenance set (by condition C2). Still, there 733 may be additional nodes that are not upstream of a reaction but do maintain it, such as, for 734 example, an enzyme that is necessary for a reaction flux but is not consumed in the reaction. 735 For SRNs, the only possible scenario for a reaction 5 to have no maintenance set is the case 736 where the only upstream node is the environment { }. For example, consider the reaction 737 5 : { } → ] , and the associated flux function 5 ( ) = 0 + \ . Biologically, 0 and \ are 738 redundant transporters that import ] from the environment into the system. For this case, 739 ( 5 ) is empty since 0 and \ are redundant in the flux function. In this case, if we split the 740 reaction 5 into two reactions 50 and 5\ and associate two flux functions 50 ≡ 0 , 5\ ≡ \ , 741 then we have ( 50 ) = 0 and ( 5\ ) = \ . 742 Practically, we can always modify the reaction set and split the flux function such that the 743 maintenance set is nonempty under this new representation. This can be done by noticing that 744 5 ( ) = 5 ( ) = 5 ( ) 0 + ⋯ + 5 ( ) 2 and by setting 5D ( ) = 5 ( ) D . It can be checked 745 that 5D ( ) are scalable functions and ( 5D ) is nonempty. In general, the method to split a 746 flux function is not unique, and we would like to choose the representation that best fits the 747 biological intuition. In the following discussion, we assume that all reactions have a nonempty 748 maintenance set. 749 and test if ( 5 ) ∈ £ D ¤. Remove every reaction 5 such that ( 5 ) ∉ £ D ¤, and 764 obtain a smaller reaction set Dƒ0 . 765 (2) Repeat the above procedure until either D is empty or an autocatalytic circuit is found. If 766 the algorithm reaches an empty set, then there is no autocatalytic circuit in 0 . 767 Proof: Suppose has no autocatalytic subset, the above algorithm will yield an empty set. 768 Suppose has an autocatalytic circuit ′, we need to show every reaction 5 ∈ ′ will not be 769 removed by this algorithm. Since ( 5 ) ⊆ ( ′) ⊆ ( ′), 5 will not be removed at step 770 = 1. Suppose that we have ′ ⊆ D up to step . We have ( 5 ) ⊆ ( ′) ⊆ ( ′) ⊆ 771 ( D ) and hence 5 will not be removed in step . Therefore, ′ will be a subset of Dƒ0 . By 772 induction on , the reactions in the subset ′ will not be removed by the algorithm. Hence, after 773 at most steps ( is the number of reactions in the system), the algorithm reaches a stable 774 reaction set * with ( * ) ⊆ ( * ) and ′ ⊆ * . ∎ 775 The following Lemma establishes a simple relation satisfied by scalable fluxes, with several 776 corollaries, which will be useful in proving our main result on autocatalytic circuits below (Thm. 777 6.5). 778 and since node growth rates are bounded by the system's long-term growth rate , which exists 797 and is finite by Thm. 2.3, the above limit is finite. ∎ 798 Corollary 6.4b: For an SRN with long-term growth rate > 0, if a flux 5 has growth rate , then 799 all of its maintenance nodes † grow at rate . 800 Proof: Let † be a maintenance node of 5 . By Lemma 6.4, 5 ( ) = 5, † ( ) † . Since 5, † ( ) is 801 a bounded non-negative function, positive growth of 5 implies that † must grow at least as 802 fast as 5 . Since 5 is growing at the same rate as the whole system and † is a component of the 803 system, † must be growing exactly at the rate . ∎ 804 Theorem 6.5: Consider an SRN ( , , ) with an initial condition (0), which is regular for an 805 ergodic measure . If the long-term growth rate is positive, then there exists at least one 806 autocatalytic circuit in the network such that every node D ∈ ( ) has a long-term growth 807 rate . 808 Proof: By Cor. 6.4a, all node growth rates [ † ] exist. We provide a recursive procedure to find 809 an autocatalytic circuit satisfying Eq. (29). Given a network with growth rate > 0, there is at 810 least one node † with [ † ] = , and it is clear that there exists at least one reaction (90) ∈ 811 ( † ) with [ 5 ] = . We denote 9 ≡ { (90) }. Now, we define 820 noting that all fluxes in Dƒ0 have growth rate . By this procedure, the size { D } monotonically 822 increases. Since the total number of reactions in the network is finite, there exists a maximal 823 collection such that is invariant under this procedure. Next, we must show that ( ) ⊆ 824 ( ), i.e., is an autocatalytic circuit. 825 Given an arbitrary reaction 5 ∈ , suppose that this reaction is recruited into at step for 826 the first time. For each node ? ∈ ( 5 ), there exists at least one reaction that will be 827 recruited into at step + 1, denoted as [ ? ] (if this reaction does not yet belong to D ). 828 Hence, ? ∈ ( [ ? ] ) ⊂ ( ). This shows ( 5 ) ∈ ( ) for an arbitrary reaction 5 ∈ 829 , and therefore, is an autocatalytic circuit. 830 By construction, every flux 5 associated with 5 ∈ has a growth rate > 0. By Cor. 6.4b, all 831 maintenance nodes of has a growth rate . ∎ 832 Remark 2: The condition > 0 is essential in Theorem 6.5. If the system has zero or negative 836 long-term growth rate, the statements in Propositions 6.4a and 6.4b no longer hold. That is, the 837 growth rate of each node in an autocatalytic circuit does not necessarily need to be the same. 838 This method is remarkable in the sense that we do not need to solve the differential equation 839 (either analytically or numerically). The procedure only requires that we inspect each flux 840 function and obtain its maintenance set. The remaining computation is purely algebraic. It is the 841 flux network connection and the property of positive growth rate that constrain the growth rate 842 between nodes and fluxes. 843 844 Example 6.6: Random networks 845 To characterize random networks, we define 5SlØ as the probability that a random reaction 846 network is autocatalytic, which depends on its statistical properties. Intuitively, for Eq. 6 to be 847 satisfied, one should decrease the average number of ( 5 ) and increase the average number 848 of ( 5 ). Furthermore, given a random reaction 5 , the chance of a given node to belong to 849 ( 5 ) or ( 5 ) could be different. Intuitively, if some nodes acted as "hubs" to maintain 850 many other nodes, the probability 5SlØ should be higher (Fig. S6A). In contrast, if a node 851 required many other upstream nodes (Fig. S6B), 5SlØ should be lower. 852 853 To test how network topology affects 5SlØ , we constructed an ensemble of random reaction 854 networks with a fixed number of nodes (n = 100), but with different connections of reactions. 855 Consistent with our expectation, we found that 5SlØ transitions from 0 to 1 as the number of 856 reactions increases in the network (Fig. S6). To see how hub-like connections affect 5SlØ , we 857 randomly sampled the maintenance and downstream nodes of every reaction from power-law 858 distributions with exponents 4l and T+ . We found that 5SlØ is strongly affected by both 859 exponents. For example, network connections with uniform probability ( 4l = T+ = 0) usually 860 required 150 or more reactions to be autocatalytic (Fig. S6A). However, for a hub-like 861 maintenance set with skewed probability 4l = 2, the number of reactions required could be 862 reduced to 100. This contrasts with a hub-like downstream set with skewed probability T+ = 1, 863 for which more than 400 reactions were required for the network to be autocatalytic. 864 865
The next theorem states that for an SDE satisfying the above criteria, the function ( ) in the 899 above theorem can be found explicitly. 900 .c 4 2 Èn0 : ≥ 0.
The quantity within the parenthesis in (A4) is called the Fichera drift.
To proceed, we show that the conditions (A1-A4) hold for the scalable SDE 1. Numerical solution of ordinary differential equations 922 923 All ordinary differential equations (ODE) were integrated by MATLAB function ode45 with 924 relative tolerance = 10 '] and absolute tolerance = 10 'ì . In Fig. 2b, c, and g, 925 the trajectory ( ) was also calculated by Mathematica function NDSolve with default accuracy, 926 yielding consistent results. In Fig. S2, the trajectory ( ) was integrated by a log-transform ODE 927 for improving the accuracy when ? is close to zero. For calculating the long-term growth rate 928 , we simulated the trajectory ( ) for a sufficiently long time such that ( ) converges to its 929 attractor. After convergence, we averaged the instantaneous growth rate ( ( D )) with time 930 step = 0.01 (time unit) to obtain growth rate . To simulate a system with noise, the differential 931 equation was integrated with an additional Gaussian noise term for every integration time step. 932 933 The autocatalytic single-repressilator model has four nodes 0 to ^ (see Fig 3A and  For each species, the influxes from the environment ( 0 to ] ) are modeled by a Lotka-Volterra-1031 like quadratic function. We modified the function to be scalable. The production of secreted 1032 metabolites are linear fluxes (^ to ì ) that are proportional to the species abundance. There is 1033 a total of six cross-feeding fluxes, with each metabolite being able to feed one or both species 1034 that do not produce this metabolite. The system dynamics is governed by the following differential equations:   There are 24 nodes and 40 reactions in our autocatalytic biosynthesis model. All flux functions of 1210 the 40 reactions are scalable (see Table S1); their formulas are derived in the previous section. 1211

Calculation of Poincare section and bifurcation analysis
The 24 nodes, including 12 metabolites (red) and 12 polymers (green) are shown in the figure  1212 below. The external nutrients 0 , \ , ] belong to the environment and are not part of the system 1213 nodes. The 40 fluxes are indicated by the arrows in the figure. Note that many fluxes are coupled with 1218 ATP → ADP reactions and hence have ATP as an upstream node and ADP as a downstream 1219 node. For clarity, all ATP-related arrows are not drawn; instead, the ATP coupling is indicated by 1220 the "yellow coin" symbol. The full stoichiometry of the reaction is listed in Table S2.  1221  1222 In addition, many reactions require enzyme catalysis, indicated by blue letters next to the 1223 reactions. These enzymes are required for flux, but are not consumed by the reaction. Hence, 1224 they do not belong to the upstream or downstream nodes of these reactions. The dependence 1225 of these enzymes is shown by the flux function in Table S2. 1226 1227 (c) Parameter range of the autocatalytic biosynthesis model 1228 1229 We tested a wide range of parameter values and various initial conditions for the system. Most 1230 of the trajectories ( ) converged to a fixed point or limit cycle, resulting in convergence of the 1231 long-term growth rate. We chose a "standard parameter set" (see Table S3) in which the system 1232 trajectory ( ) converges to a steady state. Due to the simplicity of our toy model, our 1233 parameter set cannot be directly compared to realistic metabolite concentrations. However, we 1234 can compare the relative fraction of biomass of each polymer and metabolite between 1235 simulations and the cell. Note that all parameters follow the unit system of biomass vector ( ), 1236 which is unitless. 1237 1238 To compare with classical biochemical rate constants, one can use Eq. (5) presented in the main 1239 text. We found that under our standard parameter set (Table S3) For each reaction and each slot of a maintenance set, a random node was drawn from 1286 according to probability 4l ( ) and assigned to the maintenance node slot. This specified the downstream nodes (using sampling with replacement). The probability of assigning a node to a 1322 maintenance or downstream set followed a power law exponent 4l and T+ . 1323 1324 For each condition, 500 random networks were generated and tested for their autocatalytic 1325 ability using the algorithm shown in the previous section. The autocatalytic probability 5SlØ was 1326 calculated as follows: 1327 1328  ( ) converges to a steady-state, with the three species coexisting. C, Same as B except that \ = 1.1, ] = 1.2. In the long-term, ( ) converges to a steady-state, with a single species dominating. D, Same as B except that \ = 0.9, ] = 1.2. In the long-term, ( ) converges to a heteroclinic cycle and exhibits oscillations with increasing period length. E, Same as D except that a scalable noise term ℎ ∘ ? is added for =1,2,3 with noise level ℎ = 0.01. This shows that the heteroclinic cycle (in C) is robust under noise perturbation. x 1 x 2 x 3 x 1 x 2 x 3 x 2 x 3 x 3 x 2 x 1 x 1  The red line indicates the contour of the optimal catabolic flux fraction. C, Same as A, except that the ATP production stoichiometry is set to 20. D, Same as B, except that the ATP production stoichiometry is set to 20. E, Same as A, except that the ATP production stoichiometry is set to 30. F, Same as B, except that the ATP production stoichiometry is set to 30. G, Comparison between the optimal growth rate obtained by simulation and the optimal growth rate inferred from the objective function. This plot contains six different types of optimization conditions in which ATP production stoichiometry is 30, 20 and 10 per amino acid ( í'09 ), respectively.      Table S4) 0.005 (Fig. 3F, Fig. S4