Mathematical Analysis of a Prototypical Autocatalytic Reaction Network

Network autocatalysis, which is autocatalysis whereby a catalyst is not directly produced in a catalytic cycle, is likely to be more common in chemistry than direct autocatalysis is. Nevertheless, the kinetics of autocatalytic networks often does not exactly follow simple quadratic or cubic rate laws and largely depends on the structure of the network. In this article, we analyzed one of the simplest and most chemically plausible autocatalytic networks where a catalytic cycle is coupled to an ancillary reaction that produces the catalyst. We analytically analyzed deviations in the kinetics of this network from its exponential growth and numerically studied the competition between two networks for common substrates. Our results showed that when quasi-steady-state approximation is applicable for at least one of the components, the deviation from the exponential growth is small. Numerical simulations showed that competition between networks results in the mutual exclusion of autocatalysts; however, the presence of a substantial noncatalytic conversion of substrates will create broad regions where autocatalysts can coexist. Thus, we should avoid the accumulation of intermediates and the noncatalytic conversion of the substrate when designing experimental systems that need autocatalysis as a source of positive feedback or as a source of evolutionary pressure.

The importance of autocatalysis for determining the origin of life is twofold [19]. At the initial stages of prebiotic evolution, where organic building blocks accumulated, autocatalysis is a possible solution for the "mixtures" problem [20][21][22], which involves a low abundance of any particular reaction product from diverse starting materials because of the statistical distribution of the reaction products. Autocatalysis would accelerate the formation of a limited set of products and consume starting materials for forming these products, thus avoiding the formation of the complex mixture. At later stages of prebiotic evolution, autocatalysis serves as a driving force for the natural selection of information carriers [2,12,14,23]. It amplifies the fittest carriers against others. Many variations of information carriers have been proposed such as rybozimes [24], autocatalytic sets of polypeptides [1,2,12,14,22], Life 2019, 9, 42 2 of 10 or vesicles carrying compositional information [23], but autocatalysis is always a driving force behind the evolution of these species.
Direct autocatalysis, which is a process whereby a product directly catalyzes its own production through formation of only short-living intermediates, is often differentiated from network autocatalysis, where multiple stable products cooperatively accelerate their own production [25][26][27][28][29]. Because of its high mechanistic versatility compared with direct autocatalysis, network autocatalysis has become a basis for a variety of models for determining the origin of life [1,2,11,12,14,23,30]. Nevertheless, the functional properties (e.g., the ability to evolve) of a particular network depend on their kinetic behavior which [31,32], consequently, depend on the structure of the network [33]. Therefore, it is important to analyze the kinetics of specific and chemically plausible autocatalytic reaction networks [34][35][36].
Semenov, Whitesides, and coworkers have recently published two autocatalytic reactions that can be reduced to a catalytic cycle, followed by the noncatalytic conversion of one product of this catalytic cycle to the catalyst itself ( Figure 1) [37,38].
Life 2019, 9, x FOR PEER REVIEW 2 of 10 information carriers have been proposed such as rybozimes [24], autocatalytic sets of polypeptides [1,2,12,14,22], or vesicles carrying compositional information [23], but autocatalysis is always a driving force behind the evolution of these species. Direct autocatalysis, which is a process whereby a product directly catalyzes its own production through formation of only short-living intermediates, is often differentiated from network autocatalysis, where multiple stable products cooperatively accelerate their own production [25][26][27][28][29]. Because of its high mechanistic versatility compared with direct autocatalysis, network autocatalysis has become a basis for a variety of models for determining the origin of life [1,2,11,12,14,23,30]. Nevertheless, the functional properties (e.g., the ability to evolve) of a particular network depend on their kinetic behavior which [31,32], consequently, depend on the structure of the network [33]. Therefore, it is important to analyze the kinetics of specific and chemically plausible autocatalytic reaction networks [34][35][36].
Semenov, Whitesides, and coworkers have recently published two autocatalytic reactions that can be reduced to a catalytic cycle, followed by the noncatalytic conversion of one product of this catalytic cycle to the catalyst itself ( Figure 1) [37,38].  In the first example, we can separate a catalytic cycle of cysteamine (CSH), which catalyzed the acylation of cystamine (CSSC) by a thioester; this is followed by converting one of its byproducts, ethanethiol, into CSH ( Figure 1a) [37]. In the second example, the catalytic cycle of the coppercatalyzed azide-alkyne cycloaddition reactions is followed by the formation of a catalytically active complex from triazole derivatives (Figure 1b) [38]. Because of the abundance of catalytic reactions in Figure 1. Schemes that represent (a) the thiol-based autocatalytic reaction and (b) the copper-catalyzed azide-alkine cycloaddition-based autocatalytic reaction as a catalytic cycle coupled to a non-catalytic reaction that converts one of the products of catalytic transformation into the catalyst itself. Examples of autocatalysts are cysteamine (CSH) and a copper complex (Cu(I)TATZ).
In the first example, we can separate a catalytic cycle of cysteamine (CSH), which catalyzed the acylation of cystamine (CSSC) by a thioester; this is followed by converting one of its byproducts, ethanethiol, into CSH ( Figure 1a) [37]. In the second example, the catalytic cycle of the copper-catalyzed azide-alkyne cycloaddition reactions is followed by the formation of a catalytically active complex from triazole derivatives (Figure 1b) [38]. Because of the abundance of catalytic reactions in nature, we speculate that these motifs might be among the most common autocatalytic motifs in simple (i.e., non-biological) reaction networks.
Here, we would like to analyze the kinetic behavior of one type of these motifs where the Michaelis-Menten-type catalysis is coupled to one additional irreversible reaction ( Figure 2).
We decided to analyze this particular motif because the Michaelis-Menten scheme can be applied to many catalytic reactions and possibility to reduce extended motifs to this motif by considering only rate-limiting steps. We would like to determine under which conditions this reaction network can be used in place of quadratic autocatalysis in experiments investigating chemical evolution and chemical systems with nonlinear kinetics and whether the kinetics of this network will always cause mutual exclusion of competing replicators. Schemes that represent (a) the thiol-based autocatalytic reaction and (b) the coppercatalyzed azide-alkine cycloaddition-based autocatalytic reaction as a catalytic cycle coupled to a noncatalytic reaction that converts one of the products of catalytic transformation into the catalyst itself. Examples of autocatalysts are cysteamine (CSH) and a copper complex (Cu(I)TATZ). In the first example, we can separate a catalytic cycle of cysteamine (CSH), which catalyzed the acylation of cystamine (CSSC) by a thioester; this is followed by converting one of its byproducts, ethanethiol, into CSH ( Figure 1a) [37]. In the second example, the catalytic cycle of the coppercatalyzed azide-alkyne cycloaddition reactions is followed by the formation of a catalytically active complex from triazole derivatives ( Figure 1b) [38]. Because of the abundance of catalytic reactions in

Analysis of Kinetics for a Network with an Infinite Supply of Substrates
Equations (1)-(3) with constant S 1 and S 2 describe the network in this approximation: Let us first consider two of the simplest cases: (i) when quasi-steady-state approximation (QSSA) can be applied to ( dI dt = 0) and P ( dP dt = 0) and (ii) when QSSA can be applied to A ( dA dt = 0) and P ( dP dt = 0). Simple calculations (see the Appendix A) show that in the first case the reaction is perfectly autocatalytic with the effective rate constant k e = k 1 S 1 k 2 /(k 2 + k −1 ): Obviously, in a situation with a limited amount of S 1 , the reaction will behave as quadratic autocatalysis with rate constant k e .
In the second case, the reaction is entirely autocatalytic for I with k 2 as the autocatalytic constant (see the Appendix A): Let us next look at a situation where we apply QSSA only to I. This situation is important for experimental systems because it has been shown that many catalytic reactions follow Michaelis -Menten kinetics, which implies QSSA for I. The equations (1)-(3) can be reduced to the second order Equation (6) on P: From 6, with initial conditions P(0) = 0 and A(0) = A 0 , we can derive an expression for A (see the Appendix A for details).
Life 2019, 9, 42 4 of 10 Note that for physical (positive) rate constants and concentrations, λ 1 is always positive and λ 2 is always negative. The first term in Equation (7) has a positive exponent (λ 1 ) and has the higher coefficient in front of the exponent other than the second term, and therefore, it will dominate from the beginning. Chemically, this equation means that if intermediates of the catalytic cycle do not accumulate in significant amounts, the reaction will behave as exponential autocatalysis from the beginning of the experiment and any deviations from exponential growth will decay over time.
Next, we will examine a situation where we do not have an accumulated product of the catalytic cycle because of the substantially high rate of its conversion to the autocatalyst. This situation is described by QSSA with dP dt = 0 and cannot be reduced to the second order differential equation; instead, we need to solve a system involving two equations: For the initial conditions, A(0) = A 0 and I(0) = 0; this system has a solution: Both the C 1 and C 2 coefficients are positive and C 1 + C 2 = 1 (for an exact expression of C 1 and C 2 , see the Appendix A). The value λ 1 is always positive and λ 2 is always negative for any positive rate constants. Therefore, the value of the decaying term, A 0 C 2 e λ 2 t , will not exceed that of A 0 . These calculations indicate that if in an experimental system the product P does not accumulate and A 0 is much smaller than the concentration of the substrate, S 1 , this system will behave as an almost perfect quadratic autocatalysis.
Finally, we will briefly examine a situation where we do not apply any QSSA. For an autocatalyst A, the solution has a general form: If any part of λ 1-3 is positive, then A will grow exponentially after an initial lag period (the term with a positive exponent will dominate in Equation (11)). As we show in the Appendix A, because rate constants are positive, one of the λ 1-3 values must be positive. Thus, independently of rate constants and after some lag period, this reaction network will produce exponential growth. Experimental systems are limited by the amounts of the substrates S 1 and S 2 and might not have time to reach an exponential phase, especially with a high initial concentration of A.

Competition of the Autocatalysts of Two Different Autocatalytic Networks for Common Substrates
To describe a practically interesting situation, we analyzed the competition between two autocatalytic networks in a continuously stirred tank reactor (CSTR). Network 1 consists of A 1 , I 1 , and P 1 and network 2 consists of A 2 , I 2 , and P 2 ; they compete for common substrates S 1 and S 2 . Here, the system is defined by: Life 2019, 9, 42 5 of 10 where S 10 and S 20 are the concentrations at which S 1 and S 2 are supplied to the reactor. We should mention at this point that if networks compete only for substrate S 1 and S 2 is considered to be in an unlimited supply (i.e., we consider Equations (12)-(18) with S 2 being constant), autocatalysts A 1 and A 2 cannot coexist at a steady state if any difference between the rate constants of the reactions of the networks exists (see the Appendix A). The model with constant S 2 describes experimental systems where S 2 is in a big excess in relation to S 1 or where conversion of P to A is a monomolecular reaction.
If S 2 is variable, we cannot draw a simple conclusion about coexistence and need to use numerical analysis. We analyzed 12-19 using Mathematica script (see the Appendix A). Figure 3a shows concentrations of A 1 and A 2 at t = 2000 where, in most cases, the system reaches a steady state. We set all constants to unity, the initial concentrations of A 1 and A 2 to 0.001, k 0 to 0.1, and varied k 1 and k 3 , which characterize reactions with substrates. The graph has two characteristic futures: (i) autocatalysts A 1 and A 2 do not coexist, if one has a nonzero concentration and another falls to zero; (ii) competition is not sensitive to k 3 as soon as k 3 is sufficiently high. These features mean that these networks can undergo Darwinian evolution and that this evolution will be more sensitive to improvements in k 1 than in k 3 .
An important difference between model 12-19 and plausible chemical systems is the presence of the noncatalytic conversion of substrates to autocatalysts in many experimental systems. To consider these reactions, we have to add k 4 S 1 and k 4 S 1 to Equations (12) and (13) correspondingly, and subtract both these terms from Equation (18). We performed a numerical analysis of the modified equations with the same parameters as in Figure 3a and with k 4 and k 4 set to 0.01 (Figure 3b). The graph shows that autocatalysts can coexist and the region of coexistence is higher for low k 3 values. Thus, we should avoid accumulating P in experimental systems with noncatalytic background reactions. We also explored how k 2 and k −1 influence the competition between two replicators (Figure 3c,d). The plots demonstrate that k 2 , which represents k cat in a classical Michaelis-Menten scheme, has a stronger influence on the competition than k 3 does, but it is less important than k 1 . Interestingly, simultaneously varying k 2 and k −1 produces almost a symmetrical plot, which indicates their equal contribution to the competition between replicators (Figure 3d). Overall, the data indicate that k 1 , which is responsible for initiating the catalytic cycle, has the greatest effect on the competition, whereas k 3 , which is responsible for an axillary reaction that produces an extra molecule in the catalyst, has a minimum effect.
contribution to the competition between replicators (Figure 3d). Overall, the data indicate that k1, which is responsible for initiating the catalytic cycle, has the greatest effect on the competition, whereas k3, which is responsible for an axillary reaction that produces an extra molecule in the catalyst, has a minimum effect.

Figure 3. Competition between replicators A1 and A2 in continuously stirred tank reactor (CSTR). (a)
The system is based on equations 12-19. All reaction rate constants except k1 and k3 are set to 1; k0 is set to 0.1, S10 and S20 are set to 1, and k1 and k3 are varied from 0.1 to 2. (b) When considering the noncatalytic formation of A1 and A2 with k4 = k4' = 0.01, other parameters are identical to a. (c) The system is based on equations 12-19. All reaction rate constants except k1 and k2 are set to 1; k0 is set to 0.1, S10 and S20 are set to 1, and k1 and k2 are varied from 0.1 to 2. (d) The system is based on equations 12-19. All reaction rate constants except k−1 and k2 are set to 1; k0 is set to 0.1, S10 and S20 are set to 1, and k−1 and k2 are varied from 0.1 to 2. In all plots, the concentrations of replicators A1 and A2 are plotted at t = 2000, A1(0) = 0.001, A2(0) = 0.001.

Conclusions
In this work, we determined, as precisely as possible, where and how to use reactions in experimental systems, which are based on the scheme shown in Figure 2. The results provide two main conclusions: (i) As soon as an intermediate of catalytic I or an intermediate product P does not accumulate in the reaction (at least one of them can be described by QSSA), the reaction kinetics does not deviate from exponential autocatalysis after a short lag period. (ii) The competition between these networks results in the mutual exclusion of autocatalysts if the noncatalytic formation of autocatalysts is negligible. Therefore, although these networks are not direct autocatalysts from a mechanistic perspective, they will behave as simple quadratic autocatalysts in most experimental systems.
Finally, if some variable information in autocatalyst A is transferred to product P and then retained during conversion of P to A, the network fulfills the conditions for Darwinian evolution. Interestingly, an experimental system from Otto's group is [39,40], in a way, already following this mechanism. The growth of the supramolecular stacks, which is catalyzed by the terminus of the stack, (a) The system is based on Equations (12)- (19). All reaction rate constants except k 1 and k 3 are set to 1; k 0 is set to 0.1, S 10 and S 20 are set to 1, and k 1 and k 3 are varied from 0.1 to 2. (b) When considering the noncatalytic formation of A 1 and A 2 with k 4 = k 4 = 0.01, other parameters are identical to a. (c) The system is based on Equations (12)- (19). All reaction rate constants except k 1 and k 2 are set to 1; k 0 is set to 0.1, S 10 and S 20 are set to 1, and k 1 and k 2 are varied from 0.1 to 2. (d) The system is based on Equations (12)- (19). All reaction rate constants except k −1 and k 2 are set to 1; k 0 is set to 0.1, S 10 and S 20 are set to 1, and k− 1 and k 2 are varied from 0.1 to 2. In all plots, the concentrations of replicators A 1 and A 2 are plotted at t = 2000, A 1 (0) = 0.001, A 2 (0) = 0.001.

Conclusions
In this work, we determined, as precisely as possible, where and how to use reactions in experimental systems, which are based on the scheme shown in Figure 2. The results provide two main conclusions: (i) As soon as an intermediate of catalytic I or an intermediate product P does not accumulate in the reaction (at least one of them can be described by QSSA), the reaction kinetics does not deviate from exponential autocatalysis after a short lag period. (ii) The competition between these networks results in the mutual exclusion of autocatalysts if the noncatalytic formation of autocatalysts is negligible. Therefore, although these networks are not direct autocatalysts from a mechanistic perspective, they will behave as simple quadratic autocatalysts in most experimental systems.
Finally, if some variable information in autocatalyst A is transferred to product P and then retained during conversion of P to A, the network fulfills the conditions for Darwinian evolution. Interestingly, an experimental system from Otto's group is [39,40], in a way, already following this mechanism. The growth of the supramolecular stacks, which is catalyzed by the terminus of the stack, is the catalytic step with information transfer; breaking the stack, which generates an extra terminus, is a P to A step that retains the information. Acknowledgments: E.V.S. acknowledges the Russian Government support Grant 08-08.

Conflicts of Interest:
There are no conflicts of interest.
Steady-state conditions for 12-19 are as follows: If we consider S 2 to be constant and the reaction rates are different for two networks, this system of equations does not have solutions where both A 1 and A 2 are nonzero. VI.