## Abstract

A crucial question within the fields of origins of life and metabolic networks is whether or not a self-replicating chemical reaction system is able to persist in the presence of side reactions. Due to the strong nonlinear effects involved in such systems, they are often difficult to study analytically. There are however certain conditions that allow for a wide range of these reaction systems to be well described by a set of linear ordinary differential equations. In this article, we elucidate these conditions and present a method to construct and solve such equations. For those linear self-replicating systems, we quantitatively find that the growth rate of the system is simply proportional to the sum of all the rate constants of the reactions that constitute the system (but is nontrivially determined by the relative values). We also give quantitative descriptions of how strongly side reactions need to be coupled with the system in order to completely disrupt the system.

## 1 Introduction

To study origins of life, it is a crucial question whether self-replicating systems, or autocatalytic networks, are able to emerge and persist in a chemical reaction system [7, 9]. Reactions in an autocatalytic network are always coupled in an intricate way. Thus, when studying the behavior of an autocatalytic network, simplifications are always necessary. As the behavior of a reaction network is governed by the topology of the network (including whether there is coupling between reactions, and how strong the coupling is) and the kinetics [19], the maximum simplification is to only consider whether the coupling exists.

This technique has been successfully applied to the field of origins of life, for example, in investigating whether autocatalytic networks or cycles exist in a random catalytic reaction network [46]. Besides, it has also been used to investigate the robustness of the huge cell-cycle network in yeast [2, 10], and the small functional motifs in cellular control networks [18, 19], although in the latter case, three types of interactions were considered: activation, inhibition, and none.

Nevertheless, the coupling strength and the kinetics are critical indeed [3, 11, 13, 20, 22, 23]. For example, recent studies showed that although autocatalytic networks or cycles can easily emerge from graphical analysis, most cycles are incapable of fixating dynamically if the kinetics is taken into account [3, 22, 23]. Empirical work on the Azoarcus ribozyme system has also suggested that different reaction rates make certain cycles win over others, although both are autocatalytic cycles graphically [13, 20]. A question thus arises: How does the growth pattern of the whole autocatalytic cycle depend on the different reaction rates of the system?

Another aspect of autocatalytic networks or self-replicating systems is that they are always immersed in a larger network. There are studies describing how to graphically discern an autocatalytic network or motif from a large network [5, 19], but an autocatalytic motif might be disrupted if coupled with side reactions [8, 17]. For example, an empirical study showed that networks behave significantly differently in isolation from when they are exposed to another molecular network (since some pathways cannot be realized due to competition for shared reactants from side reactions or other networks) [1]. For another example, the reverse citric acid cycle was proposed to be the candidate of autocatalytic networks on the primitive Earth, but it is still an open question whether or not this cycle can be free of disruption from the many side reactions in the primitive atmosphere [14].

Intuitively, if the reaction rate of the side reaction is extremely small, the dynamics should go back to the scenario with no side reaction. If the side reaction, on the other hand, is fast enough (i.e., the coupling is strong enough), it may completely disrupt the autocatalytic network. Whether side reactions could disrupt a network essentially depends on how strongly they are coupled, rather than merely whether they exist or not. A question thus arises: What is the threshold of such a coupling strength, if it ever exists?

Although many have simulated the growth pattern of autocatalytic networks numerically [3, 12, 15], there are very few studies analytically investigating how different rates of the reactions in the autocatalytic network affect the growth pattern of the whole network ([21] is one of the few examples). This is probably due to the strong nonlinear effects that are often involved in autocatalytic networks. However, there is a wide range of autocatalytic networks, or self-replicating systems, that can be reduced to linear systems under certain conditions. In that case, the dynamics of the system can be approximately described by a set of linear ordinary differential equations (ODEs), which can be solved analytically. The “certain conditions,” which we shall see later in detail, either refer to the situation where there is an infinite amount of resources, or equivalently refer to the initial stage of self-replication when the amount of resources is much more than the amount of other reacting molecules involved.

In this article, we employ the artificial chemistry model presented in a previous one [12], which is a convenient framework to construct various chemical reaction systems. In Section 2, we introduce the method to construct and solve the ODEs of the self-replicating systems that can be reduced to linear systems, and then apply this method to investigate the effects of different reaction rates on the growth pattern of those self-replicating systems. In Section 3, we investigate how strongly side reactions need to be coupled with the self-replicating system in order to disrupt it partially or completely. Finally, conclusions are drawn in Section 4.

## 2 Method

Here we briefly recap the artificial chemistry model presented in [12]: (1) Each type of molecule is defined by its integer mass i, and denoted $i¯$; (2) only the synthesis reaction $i¯$ + $j¯$$i+j¯$ and the decomposition reaction $i+j¯$$i¯$ + $j¯$ are possible; (3) it is assumed that every reaction is elementary (mass action kinetics thus applies), all molecules are ideally gaseous, and the whole system is kept at constant pressure and temperature. Therefore, the reaction rate for a synthesis reaction $i¯$ + $j¯$$i+j¯$ can be written as (refer to Supplementary Information S1 in [12] for detailed derivations)
$γ+ij=ω+ij⋅Ni⋅Nj/N$
(1)
where ω+ij is the reaction rate constant (in units of s−1), Ni is the number of molecule $i¯$ inside the whole system, and N = ∑iNi is the total number of molecules. For a decomposition reaction $i+j¯$$i¯$ + $j¯$, the reaction rate is given by
$γ−ij=ω−ij⋅Ni+j$
(2)
Note that instead of using molecule concentrations, here we use molecule counts. Consequently, the reaction rates (γ+ij and γij) are in units of s−1, rather than mol · L−1 · s−1 (concentration per second) as commonly used.

As we know, in mass action kinetics, elementary synthesis reactions are second-order while elementary decomposition reactions are first-order. If we investigated the dynamics in terms of molecule concentrations, the units of the rate constants would be different: The former is in mol−1 · L · s−1, while the latter is in s−1. We thus cannot compare them directly. But the comparison between rate constants is the key to investigating the question we are interested in, namely, how strongly side reactions need to be coupled with the self-replicating system in order to disrupt it partially or completely. Instead, by using molecule counts, the unit of the reaction rate constant for the synthesis reaction (ω+ij) and that for the decomposition reaction (ωij) can be calculated by the same formula β exp (−κΔG) (thus having the same units), where β and κ are physical constants, and ΔG is the activation energy of the corresponding reaction. The comparison will then be valid, and correspond to the comparison between intrinsic properties of the reactions.

Now, we first analytically investigate the self-replicating system
$1¯+2¯→3¯1¯+3¯→4¯4¯→2¯+2¯$
(3)
given resource molecule $1¯$. It corresponds to the formose reaction that involves the formation of sugars from formaldehyde [12]. Specifically, molecule $1¯$ is formaldehyde, $2¯$ is glycolaldehyde, $3¯$ is glyceraldehyde, and $4¯$ is tetrose.

### 2.1 When Reaction Rate Constants are Identical

Here we assume that all of the reaction rate constants are identical, that is, ω+12 = ω+13 = ω−22 = ω. The dynamics of the system are hence governed by the following nonlinear ODEs:
$N˙2=ω−N1N2N+2N4N˙3=ω−N1N3N+N1N2NN˙4=ω−N4+N1N3N$
Note that Ni(t) and N(t) are functions of t, but for simplicity, we omit writing the t. Now we assume that
• (i)

there is an infinite reservoir of the resource molecule $1¯$, so that N1 always keeps constant (i.e., N1Q), and

• (ii)

Q is much larger than the number of all other molecules.

Then we have the linearization condition:
$N1N≡QQ+N2+N3+N4≈QQ=1$
(4)
Therefore, the ODE system above becomes linear, namely,
$N˙2=ω−N2+2N4N˙3=ω−N3+N2N˙4=ω−N4+N3$
which can be written in the matrix form
$N˙2N˙3N˙4=ω−1021−1001−1N2N3N4≔AN2N3N4$
(5)
There is a standard way to solve Equation 5 [16]. Briefly, the asymptotic dynamics of Ni are governed by the matrix A's maximum positive eigenvalue λ1, that is,
$Nit≐cieλ1t$
where ci is a constant, decided by the initial condition (it is interesting to notice that N2, N3, and N4 all increase exponentially at rate λ1). The characteristic equation of A is
$λ3+3ωλ2+3ω2λ−ω3=0$
and the eigenvalues correspond to its solutions
$λ1=(23−1)ωλ2≈−1.63−1.09jωλ3≈−1.63+1.09jω$
where j is the unit imaginary number.

### 2.2 When Reaction Rate Constants are Different

Now we analytically investigate the case where the reaction rate constants of the system (3) are different. We denote the sum of all the three constants as Ω := ω+12 + ω+13 + ω−22 . Then we can write
$ω+12=aΩω+13=bΩω−22=cΩ$
(6)
where naturally a + b + c = 1 and 0 ≤ a, b, c ≤ 1. We assume that the condition (4) still holds, so the ODE system in this case can be written as
$N˙2N˙3N˙4=Ω−a02ca−b00b−cN2N3N4≔ΩBN2N3N4$
(7)
Then the characteristic equation of the matrix ΩB is
$λ3+Ωλ2+ab+ac+bcΩ2λ−abcΩ3=0$
(8)
Note that a = b = c = 1/3 (consequently, ω = Ω/3) corresponds to the case where all rate constants are identical.
In the Supplementary Information (https://www.wuyichen.org/si-yuliu-artificial-life-2020), we show that the maximum positive eigenvalue, λ1, of the matrix ΩB is linearly proportional to Ω, that is,
$λ1=Ωλ1*$
(9)
where $λ1*$ is the maximum positive eigenvalue of the matrix B. That is to say, the sum of the rate constants (Ω) only affects the growth rate (λ1) of the self-replicating system linearly; and only the relative rate constants (characterized by the matrix B) govern the growth rate in a nontrivial way. Therefore, to capture the behavior of Equation 7, we only need to study the matrix B, which has only two free parameters, a and b. To illustrate this, we draw a heatmap of the maximum (actually the only) positive eigenvalue $λ1*$ along with a and b, as shown in Figure 1.
Figure 1.

Heatmap of the maximum positive eigenvalue $λ1*$ along with a and b, for the matrix B in Equation 7.

Figure 1.

Heatmap of the maximum positive eigenvalue $λ1*$ along with a and b, for the matrix B in Equation 7.

We see that $λ1*$ varies with a and b, but the maximum is located at a = b = 1/3 (consequently, c = 1/3 also). That is, the maximum growth rate is achieved when each reaction in this linear self-replicating system has the same rate constant. When the rate constants are different, $λ1*$ still exists (i.e., there is always one eigenvalue that is positive). Therefore, no matter how different the rate constants are, the whole system is always self-replicating, although with smaller growth rates. But when the rate constants are different in orders of magnitude, corresponding to the edges of the colored triangle in Figure 1 (namely a ≈ 0, b ≈ 0, or a + b ≈ 1), the growth rate might be too small to detect in practice.

### 2.3 This Method can be Applied to a Wide Range of Systems

If any chemical reaction system (i) is constructed under the artificial chemistry framework presented in [12] (of which the main points have been described at the beginning of Section 2 of that article), and (ii) satisfies the criterion that for every reaction, one and only one type of its reactants comes from the products of other reactions in the system, this system can be linearized and analyzed by the method presented above.

In [12, Section 2], we applied this method to one of the chemical reaction systems that satisfy an extra criterion, namely, that there are some types of intermediate molecules and the number of times they appear on the product side is larger than that on the reactant side (note that an intermediate molecule is defined to be any molecule that appears on both the reactant side and the product side of the whole system). This type of chemical reaction system is guaranteed to be self-replicating, that is, the numbers of some molecules can grow exponentially in proper conditions [12]. We used this method to investigate the effects of different reaction rates on the growth pattern of self-replicating systems.

In Section 3, we will apply this method to self-replicating systems with side reactions, to investigate how strongly side reactions need to be coupled with the self-replicating system in order to disrupt it partially or completely. Note that in our context, side reactions are defined to be those extra reactions that consume intermediate molecules of a self-replicating system.

For those chemical reaction systems that can be linearized (the system in Equation 3 is one such example), the linear ODE system can be constructed in the same way as in Equation 5, as long as the linearization condition holds as given in Equation 4. If we assume that (i) the number of the resource molecule always keeps constant, that is, there is an infinite reservoir of resources, and (ii) this constant is much larger than the total number of all other molecules, then this linearization condition is guaranteed. The assumption of the infinite resource is unrealistic, though. Nevertheless, during the initial stage of self-replication, the resource molecules are much more numerous than other reacting molecules indeed (e.g., QN2 + N3 + N4 for the system in Equation 3). The infinite-resource assumption is thus valid during this stage.

To summarize, a linear ODE system such as Equation 5 is able to describe the mean-field dynamics of a chemical reaction system (i) when there is an infinite amount of resources, (ii) when the number of resource molecules is much greater than the number of other reacting molecules involved, or (iii) during the initial stage of self-replication.

## 3 Effects of Side Reactions

In this section, we investigate how strongly side reactions need to be coupled with the self-replicating system in order to disrupt the whole network. Many side reactions might be considered, but here we only consider those that do not violate the linearization condition mentioned above.

### 3.1 First Example System

Now consider the formose reaction in Equation 3 coupled with a side reaction, described by the following equations:
$1¯+2¯→3¯1¯+3¯→4¯4¯→2¯+2¯2¯+5¯→7¯$
(10)
We still denote Ω = ω+12 + ω+13 + ω−22. For the fourth reaction, we denote ω+25 = κΩ, where κ ≥ 0.
Here we assume infinite reservoirs of resource molecules $1¯$ and $5¯$, and N1N5QN2 + N3 + N4 + N7. The reason why we need a reservoir of $5¯$ is that the side reaction $2¯$ + $5¯$$7¯$ can never occur if there is no $5¯$. Now, we have
$N1N=N1∑i=1,2,3,4,5,7Ni≈N1N1+N5=12$
(11)
and similarly, N5/N ≈ 1/2. By this approximation, we can write down the linear ODE system
$N˙2N˙3N˙4=Ω−κ/2−a/202ca/2−b/200b/2−cN2N3N4≔ΩDN2N3N4$
We need to find the maximum positive eigenvalue $λ1*$ of the matrix D. For different values of κ, we can obtain different heatmaps for $λ1*$. We put them together, and then obtain a three-dimensional heatmap, shown in Figure 2.
Figure 2.

The three-dimensional heatmap of the maximum positive eigenvalue $λ1*$ along with a, b, and κ, for the system (10). The upper triangular region is empty because of the constraint c = 1 − ab ≥ 0. There are also empty parts in the lower triangular region. Those empty parts represent the fact that for those (a, b), none of the three eigenvalues of the matrix D has a positive real part, that is, the self-replicating system is completely disrupted.

Figure 2.

The three-dimensional heatmap of the maximum positive eigenvalue $λ1*$ along with a, b, and κ, for the system (10). The upper triangular region is empty because of the constraint c = 1 − ab ≥ 0. There are also empty parts in the lower triangular region. Those empty parts represent the fact that for those (a, b), none of the three eigenvalues of the matrix D has a positive real part, that is, the self-replicating system is completely disrupted.

When κ = 0, the system in Equation 10 should go back to the system in Equation 3, but because we have two resource molecules here (making the approximations in Equation 4 and Equation 11 slightly different), there are small differences in $λ1*$. One major difference is that here the maximum $λ1*$ is approximately at (a, b, c) ≐ (0.37, 0.37, 0.26), instead of (a, b, c) = (1/3, 1/3, 1/3) as in Figure 1. But generally speaking, $λ1*$ of the system in Equation 10 with κ = 0 has almost the same pattern as that of the system in Equation 3.

When κ is larger, making the side reaction more strongly coupled, there are three phenomena to be noticed. (1) There are smaller regions of (a, b) having positive real eigenvalues, that is, smaller regions where self-replication can occur. For a certain κ, the region corresponds to the colored triangle (for convenience of later discussions, we denote the area of the region as S(κ)). (2) There exists a critical value κ = 1 beyond which the whole system in Equation 10 is completely unable to self-replicate, no matter what the value of (a, b) is. (3) For a fixed value of (a, b), $λ1*$ gets smaller (if it exists), meaning that the self-replication gets slower.

For convenience, we define the relative area s(κ) = S(κ)/S(0). From Figure 2 we see qualitatively that s(κ) shrinks with increasing κ. We further found that s(κ) shrinks quadratically, that is, the square root of s(κ) decreases linearly.

Now, we are going to use numerical simulations to show strong evidence of this statement (unfortunately, our attempt to prove it analytically failed). We use p(κ) to denote, for a certain κ, the probability of a randomly chosen set of (a, b) (under the constraint 0 ≤ a, b ≤ 1 and c = 1 − ab ≥ 0) that allows the whole system to self-replicate (equivalently, a positive eigenvalue exists). Evidently, p(κ) = s(κ), and by simulations, p(κ) can be worked out. Figure 3 (the blue solid line) shows $pκ$, which turns out to be a straight line. Specifically,
$pκ=sκ=1−κ2if0≤κ≤10ifκ>1$
Clearly, κ = 1 is the critical value beyond which the whole system is not able to self-replicate at all, no matter what the value of (a, b) is. We do not have analytical explanations why it is a declining straight line, but we will give more examples and some exceptions later.
Figure 3.

$pκm$ for different systems. $pκ$ for the system (10) corresponds to the blue solid line. $pκ3$ for the system (12) happens to correspond to the blue solid line also. $pκ3$ for the system (13) corresponds to the red dashed line.

Figure 3.

$pκm$ for different systems. $pκ$ for the system (10) corresponds to the blue solid line. $pκ3$ for the system (12) happens to correspond to the blue solid line also. $pκ3$ for the system (13) corresponds to the red dashed line.

### 3.2 Second Example System

Now we consider the following system, where the first four reactions constitute a self-replicating system (as shown in [12]) while the fifth one is a side reaction, given the resource molecules $2¯$ and $7¯$, and N2N7QN1 + N3 + N4 + N5 + N6 + N11:
$5¯→1¯+4¯2¯+3¯→5¯2¯+4¯→6¯6¯→3¯+3¯4¯+7¯→11¯$
(12)
By the same approach, we can write down the linear ODEs,
$N˙3N˙4N˙5N˙6=Ω−b/2002d0−κ/2−c/2a0b/20−a00c/20−dN3N4N5N6$
where Ω = ω−14 + ω+23 + ω+24 + ω−33, ω−14 = aΩ, ω+23 = bΩ, ω+24 = cΩ, ω33 = dΩ, and ω+47 = κΩ. Now we have three free parameters a, b, c (as d = 1 − abc) to determine $λ1*$ for a certain κ, which becomes hard to visualize.

Nevertheless, we have three similar observations to those for the first example system in Equation 10: (1) When the rate constant of the side reaction κ = 0, the maximum $λ1*$ is approximately at (a, b, c, d) ≐ (0.21, 0.29, 0.29, 0.21). This is because N2/N ≈ 1/2. But if N2/N ≈ 1 (that is, there is only one type of resource molecule $2¯$), the maximum $λ1*$ would be at (a, b, c, d) = (1/4, 1/4, 1/4, 1/4). (2) $pκ3$ is also a declining straight line, shown in Figure 3 (namely the blue solid line, which coincides with that for the system in Equation 10). So, κ = 1 is again the critical value. It is interesting to notice that for this system, $pκ3$ is a straight line, while for system in Equation 10, $pκ$ is a straight line. That is, in both systems, the root index is equal to the number of free parameters. (3) For fixed (a, b, c), when κ gets larger, $λ1*$ gets smaller, if it exists.

### 3.3 Third Example System

However, κ = 1 does not necessarily correspond to the critical value. The following system is one such example:
$2¯→1¯+1¯3¯→1¯+2¯1¯+4¯→5¯5¯→2¯+3¯4¯+5¯→9¯$
(13)
Here the first four reactions in Equation 13 constitute a self-replicating system (as shown in [11]), while the fifth one is a side reaction that consumes the intermediate molecule $5¯$. We further assume that the resource molecule is $4¯$ and N4QN1 + N2 + N3 + N5 + N9. We can write down the linear ODEs,
$N˙1N˙2N˙3N˙5=Ω−c2ab00−abd00−bdc00−κ−dN1N2N3N5$
(14)
where Ω = ω11 + ω12 + ω+14 + ω−23, ω−11 = aΩ, ω−12 = bΩ, ω+14 = cΩ, ω−23 = dΩ, and ω+45 = κΩ.

Similarly, we have three observations: (1) When κ = 0, the maximum $λ1*$ is approximately at (a, b, c, d) ≐ (0.24, 0.18, 0.29, 0.29), even though in this case there is only one resource molecule $4¯$. (2) $pκ3$ is a declining straight line, shown in Figure 3 (the red dashed line), and 3 is the number of free parameters. But in this case, κ = 4 is the critical value. (3) For fixed (a, b, c), when κ gets larger, $λ1*$ gets smaller, if it exists.

### 3.4 Fourth Example System

It is also possible that there is no critical value for κ. Consider the system consisting of the first four reactions in Equation 13 and another side reaction $2¯$ + $4¯$$6¯$, given N4QN1 + N2 + N3 + N5 + N6. Similarly, we can write down the linear ODEs,
$N˙1N˙2N˙3N˙5=Ω−c2ab00−κ−abd00−bdc00−dN1N2N3N5$
(15)
with similar notation to that in Equation 14, but ω+24 = κΩ.

For this system, no matter what κ is, $λ1*$ always has positive eigenvalues for any (a, b, c, d), that is, p(κ) ≡ 1. So there is no critical value for κ. Nevertheless, it is observed, as in the previous examples, that for fixed (a, b, c, d), when κ gets larger, $λ1*$ gets smaller. That is to say, this side reaction $2¯$ + $4¯$$6¯$ does reduce the growth rate of the self-replicating system, but no matter how strongly this reaction is coupled, there is always a positive growth rate, although it might be extremely small.

## 4 Conclusion

For a wide range of chemical reaction networks under certain conditions, their mean-field dynamics can be described by a linear ODE system, which can be solved analytically. We have elucidated the method to construct such an ODE system and the conditions under which the method is valid. We used this method to analyze linear self-replicating systems. But for those self-replicating systems that cannot be linearized, other approaches, such as numerical methods and stochastic simulation, have to be used, which are beyond the scope of the current article.

Using this method, we found that for all linear self-replicating systems or autocatalytic networks, every type of molecules (except the resource molecules that are consumed) grows exponentially. More surprisingly, all the exponential growth rates are identical. Furthermore, the growth rate is directly proportional to the sum of reaction rate constants, as seen in Equation 9. That is, only the relative rate constants determine the growth rate in a nontrivial way.

The growth rate of a linear self-replicating system reaches a maximum when the reaction rate constants are compatible with each other: In some cases, they are equal (e.g., the systems in Equations 3 and 12), while in other cases, they are slightly different (e.g., the system in Equation 13). Nevertheless, the maximum growth rate also depends on external conditions, for example, if there are one or two resources available, as in the systems in Equations 10 and 12.

Side reactions (those consuming intermediate molecules) always disrupt self-replicating systems. It is generic that, the larger the side reaction’s rate constant κ is, the slower the whole system grows, under otherwise identical conditions. The second effect of side reactions is less generic: The larger κ is, the fewer scenarios there are where the whole system still self-replicates; and there may be a critical value for κ such that the whole system is completely disrupted (that is, whatever the scenarios are, the system is not able to self-replicate, as in the systems in Equations 10, 12, and 13). The second effect does not always exist, though; for example, for the system in Equation 15, no matter what κ is, the whole system is always able to self-replicate, although the growth rate might be extremely small if κ is very large.

Another interesting phenomena deserves to be noticed, which we do not have analytical explanations for now. For some reaction systems (for example, those in Equations 10, 12, and 13), $pκm$ is a declining straight line where the root index m is equal to the number of free parameters. This observation is not generic for all linear self-replicating systems, but interesting enough to be further investigated.

To summarize, whether side reactions disrupt a self-replicating chemical reaction system depends not only on whether side reactions exist, but also on how strongly side reactions are coupled with the system. The method we presented in this article gives quantitative answers to the “how strongly” question, at least for a wide range of reaction networks. It sheds light on the problems caused by side reactions in studies of origins of life [9, 14, 22] and metabolic networks [1, 8, 20].

## Acknowledgments

This article is based upon work conducted while Y. L. and T. L. were in residence at Institut Mittag-Leffler in Djursholm, Sweden, during the autumn semester of 2018.

The authors declare no conflict of interest.

## References

1
Ashkenasy
,
G.
,
Jagasia
,
R.
,
,
M.
, &
,
M. R.
(
2004
).
Design of a directed molecular network
.
Proceedings of the National Academy of Sciences of the U.S.A.
,
101
(
30
),
10872
10877
.
2
Davidich
,
M. I.
, &
Bornholdt
,
S.
(
2013
).
Boolean network model predicts knockout mutant phenotypes of fission yeast
.
PLOS ONE
,
8
(
9
),
e71786
.
3
Filisetti
,
A.
,
Graudenzi
,
A.
,
Serra
,
R.
,
Villani
,
M.
,
Füchslin
,
R. M.
,
Packard
,
N.
,
Kauffman
,
S. A.
, &
Poli
,
I.
(
2012
).
A stochastic model of autocatalytic reaction networks
.
Theory in Biosciences
,
131
(
2
),
85
93
.
4
Hordijk
,
W.
,
Smith
,
J. I.
, &
Steel
,
M.
(
2015
).
Algorithms for detecting and analysing autocatalytic sets
.
Algorithms for Molecular Biology
,
10
,
15
.
5
Hordijk
,
W.
, &
Steel
,
M.
(
2004
).
Detecting autocatalytic, self-sustaining sets in chemical reaction systems
.
Journal of Theoretical Biology
,
227
,
451
461
.
6
Kauffman
,
S. A.
(
1986
).
Autocatalytic sets of proteins
.
Journal of Theoretical Biology
,
119
,
1
24
.
7
Kauffman
,
S. A.
(
2011
).
Approaches to the origin of life on earth
.
Life
,
1
,
34
48
.
8
Kreyssig
,
P.
,
Escuela
,
G.
,
Reynaert
,
B.
,
Veloz
,
T.
,
Ibrahim
,
B.
, &
Dittrich
,
P.
(
2012
).
Cycles and the qualitative evolution of chemical systems
.
PLOS ONE
,
7
(
10
),
e45772
.
9
Lancet
,
D.
,
Zidovetzki
,
R.
, &
Markovitch
,
O.
(
2018
).
Systems protobiology: Origin of life in lipid catalytic networks
.
Journal of Royal Society of Interface
,
15
,
20180159
.
10
Li
,
F.
,
Long
,
T.
,
Lu
,
Y.
,
Ouyang
,
Q.
, &
Tang
,
C.
(
2004
).
The yeast cell-cycle network is robustly designed
.
Proceedings of the National Academy of Sciences of the U.S.A.
,
101
(
14
),
4781
4786
.
11
Liu
,
Y.
, &
Sumpter
,
D. J. T.
(
2018
).
Is the golden ratio a universal constant for self-replication?
PLOS ONE
,
13
(
7
),
1
18
.
12
Liu
,
Y.
, &
Sumpter
,
D. J. T.
(
2018
).
Mathematical modeling reveals spontaneous emergence of self-replication in chemical reaction systems
.
Journal of Biological Chemistry
,
293
(
49
),
18854
18863
.
13
Mathis
,
C.
,
,
S.
,
Walker
,
S.
, &
Lehman
,
N.
(
2017
).
Prebiotic RNA network formation: A taxonomy of molecular cooperation
.
Life
,
7
(
4
),
38
.
14
Orgel
,
L. E.
(
2008
).
The implausibility of metabolic cycles on the prebiotic earth
.
PLoS Biology
,
6
(
1
),
0005
0013
.
15
Piedrafita
,
G.
,
Montero
,
F.
,
Morán
,
F.
,
Cárdenas
,
M. L.
, &
Cornish-Bowden
,
A.
(
2010
).
A simple self-maintaining metabolic system: Robustness, autocatalysis, bistability
.
PLoS Computational Biology
,
6
(
8
),
e1000872
.
16
Råde
,
L.
, &
Westergren
,
B.
(
2004
).
Mathematics Handbook for Science and Engineering
(5th ed.).
Berlin-Heidelberg
:
Springer-Verlag
.
17
Szathmáry
,
E.
(
2000
).
The evolution of replicators
.
Philosophical Transactions of the Royal Society B: Biological Sciences
,
355
(
1403
),
1669
1676
.
18
Tyson
,
J. J.
,
Chen
,
K. C.
, &
Novák
,
B.
(
2003
).
Sniffers, buzzers, toggles and blinkers: Dynamics of regulatory and signaling pathways in the cell
.
Current Opinion in Cell Biology
,
15
(
2
),
221
231
.
19
Tyson
,
J. J.
, &
Novák
,
B.
(
2010
).
Functional motifs in biochemical reaction networks
.
Annual Review of Physical Chemistry
,
61
(
1
),
219
240
.
20
Vaidya
,
N.
,
Manapat
,
M. L.
,
Chen
,
I. A.
,
Xulvi-Brunet
,
R.
,
Hayden
,
E. J.
, &
Lehman
,
N.
(
2012
).
Spontaneous network formation among cooperative RNA replicators
.
Nature
,
491
(
7422
),
72
77
.
21
Virgo
,
N.
,
Ikegami
,
T.
, &
McGregor
,
S.
(
2016
).
Complex autocatalysis in simple chemistries
.
Artificial Life
,
22
(
2
),
138
152
.
22
Walker
,
S. I.
, &
Mathis
,
C.
(
2018
).
Network theory in prebiotic evolution
. In
C.
Menor-Salván
(Ed.),
Prebiotic chemistry and chemical evolution of nucleic acids
(pp.
263
291
).
Cham
:
Springer International Publishing
.
23
Wynveen
,
A.
,
Fedorov
,
I.
, &
Halley
,
J. W.
(
2014
).
Nonequilibrium steady states in a model for prebiotic evolution
.
Physical Review E
,
89
,
022725
.

## Author notes

**

Y. L. and D. H. contributed equally to this work.