## 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 [4–6]. 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

*i*, and denoted $i\xaf$; (2) only the synthesis reaction $i\xaf$ + $j\xaf$ → $i+j\xaf$ and the decomposition reaction $i+j\xaf$ → $i\xaf$ + $j\xaf$ 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\xaf$ + $j\xaf$ → $i+j\xaf$ can be written as (refer to Supplementary Information S1 in [12] for detailed derivations)

*ω*

_{+ij}is the reaction rate constant (in units of s

^{−1}),

*N*

_{i}is the number of molecule $i\xaf$ inside the whole system, and

*N*= ∑

_{i}

*N*

_{i}is the total number of molecules. For a decomposition reaction $i+j\xaf$ → $i\xaf$ + $j\xaf$, the reaction rate is given by

*γ*

_{+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.

### 2.1 When Reaction Rate Constants are Identical

*ω*

_{+12}=

*ω*

_{+13}=

*ω*

_{−22}=

*ω*. The dynamics of the system are hence governed by the following nonlinear ODEs:

*N*

_{i}(

*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\xaf$, so that

*N*_{1}always keeps constant (i.e.,*N*_{1}≡*Q*), and - (ii)
*Q*is much larger than the number of all other molecules.

*N*

_{i}are governed by the matrix

*A*'s maximum positive eigenvalue

*λ*

_{1}, that is,

*c*

_{i}is a constant, decided by the initial condition (it is interesting to notice that

*N*

_{2},

*N*

_{3}, and

*N*

_{4}all increase exponentially at rate

*λ*

_{1}). The characteristic equation of

*A*is

*j*is the unit imaginary number.

### 2.2 When Reaction Rate Constants are Different

*ω*

_{+12}+

*ω*

_{+13}+

*ω*

_{−22}. Then we can write

*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

*B*is

*a*=

*b*=

*c*= 1/3 (consequently,

*ω*= Ω/3) corresponds to the case where all rate constants are identical.

*λ*

_{1}, of the matrix Ω

*B*is linearly proportional to Ω, that is,

*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 $\lambda 1*$ along with

*a*and

*b*, as shown in Figure 1.

We see that $\lambda 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, $\lambda 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., *Q* ≫ *N*_{2} + *N*_{3} + *N*_{4} 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

*ω*

_{+12}+

*ω*

_{+13}+

*ω*

_{−22}. For the fourth reaction, we denote

*ω*

_{+25}=

*κ*Ω, where

*κ*≥ 0.

*N*

_{1}≡

*N*

_{5}≡

*Q*≫

*N*

_{2}+

*N*

_{3}+

*N*

_{4}+

*N*

_{7}. The reason why we need a reservoir of $5\xaf$ is that the side reaction $2\xaf$ + $5\xaf$ → $7\xaf$ can never occur if there is no $5\xaf$. Now, we have

*N*

_{5}/

*N*≈ 1/2. By this approximation, we can write down the linear ODE system

*D*. For different values of

*κ*, we can obtain different heatmaps for $\lambda 1*$. We put them together, and then obtain a three-dimensional heatmap, shown in Figure 2.

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 $\lambda 1*$. One major difference is that here the maximum $\lambda 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, $\lambda 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*), $\lambda 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.

*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 −

*a*−

*b*≥ 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\kappa $, which turns out to be a straight line. Specifically,

*κ*= 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.

### 3.2 Second Example System

*N*

_{2}≡

*N*

_{7}≡

*Q*≫

*N*

_{1}+

*N*

_{3}+

*N*

_{4}+

*N*

_{5}+

*N*

_{6}+

*N*

_{11}:

*ω*

_{−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 −

*a*−

*b*−

*c*) to determine $\lambda 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 $\lambda 1*$ is approximately at (*a*, *b*, *c*, *d*) ≐ (0.21, 0.29, 0.29, 0.21). This is because *N*_{2}/*N* ≈ 1/2. But if *N*_{2}/*N* ≈ 1 (that is, there is only one type of resource molecule $2\xaf$), the maximum $\lambda 1*$ would be at (*a*, *b*, *c*, *d*) = (1/4, 1/4, 1/4, 1/4). (2) $p\kappa 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\kappa 3$ is a straight line, while for system in Equation 10, $p\kappa $ 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, $\lambda 1*$ gets smaller, if it exists.

### 3.3 Third Example System

*κ*= 1 does not necessarily correspond to the critical value. The following system is one such example:

*N*

_{4}≡

*Q*≫

*N*

_{1}+

*N*

_{2}+

*N*

_{3}+

*N*

_{5}+

*N*

_{9}. We can write down the linear ODEs,

*ω*

_{−11}+

*ω*

_{−12}+

*ω*

_{+14}+

*ω*

_{−23},

*ω*

_{−11}=

*a*Ω,

*ω*

_{−12}=

*b*Ω,

*ω*

_{+14}=

*c*Ω,

*ω*

_{−23}=

*d*Ω, and

*ω*

_{+45}=

*κ*Ω.

Similarly, we have three observations: (1) When *κ* = 0, the maximum $\lambda 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\xaf$. (2) $p\kappa 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, $\lambda 1*$ gets smaller, if it exists.

### 3.4 Fourth Example System

*κ*. Consider the system consisting of the first four reactions in Equation 13 and another side reaction $2\xaf$ + $4\xaf$ → $6\xaf$, given

*N*

_{4}≡

*Q*≫

*N*

_{1}+

*N*

_{2}+

*N*

_{3}+

*N*

_{5}+

*N*

_{6}. Similarly, we can write down the linear ODEs,

*ω*

_{+24}=

*κ*Ω.

For this system, no matter what *κ* is, $\lambda 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, $\lambda 1*$ gets smaller. That is to say, this side reaction $2\xaf$ + $4\xaf$ → $6\xaf$ 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\kappa 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

## Author notes

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