Mathematical Analysis of a Prototypical Autocatalytic Reaction Network

Abstract

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.

1. Introduction

Autocatalytic reactions are of profound importance for at least three problems in the chemical sciences: (i) chemical evolution and the origin of life, [1,2,3,4,5,6] (ii) dissipative chemical systems [7,8,9], and (iii) chemical signaling and information processing systems [10]. Autocatalytic self-replication is an element needed for Darwinian selection. Biological evolution is driven by cellular division and life most likely originates from prebiotic reactions that involve some form of self-replication and autocatalysis [1,11,12,13,14]. Dissipative systems, such as chemical and biological oscillators and Turing structures [15,16,17,18], often require autocatalytic reactions as a source of positive feedback. Finally, biochemical signaling systems need autocatalysis for signal amplification [10].

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], 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 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.

2. Results and Discussion

2.1. Analysis of Kinetics for a Network with an Infinite Supply of Substrates

Equations (1)–(3) with constant S₁ and S₂ describe the network in this approximation:

$$\frac{dA}{dt}=-k_1{S}_{1}A+(k_2+k_{-1})I+k_3{S}_{2}P$$

(1)

$$\frac{dI}{dt}=k_1{S}_{1}A-(k_2+k_{-1})I$$

(2)

$$\frac{dP}{dt}=k_{2}I-k_3{S}_{2}P$$

(3)

Let us first consider two of the simplest cases: (i) when quasi-steady-state approximation (QSSA) can be applied to ($\raisebox{1ex}{$dI$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=0$) and P ($\raisebox{1ex}{$dP$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=0$) and (ii) when QSSA can be applied to A ($\raisebox{1ex}{$dA$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=0$) and P ($\raisebox{1ex}{$dP$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=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₁S₁k₂/(k₂ + k₋₁):

$$\frac{dA}{dt}=\frac{k_1{S}_1{k}_2}{(k_2+k_{-1})}A$$

(4)

Obviously, in a situation with a limited amount of S₁, 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₂ as the autocatalytic constant (see the Appendix A):

$$\frac{dI}{dt}=k_{2}I$$

(5)

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:

$$\frac{d^{2}P}{d{t}^2}+k_3{S}_2\frac{dP}{dt}-\frac{k_3{S}_2{k}_1{S}_1{k}_2}{(k_2+k_{-1})}P=0$$

(6)

From 6, with initial conditions P(0) = 0 and A(0) = A₀, we can derive an expression for A (see the Appendix A for details).

$$\begin{array}{l}A(t)=\frac{A_0}{\sqrt{k_3{S}_2(k_3{S}_2+4k_1{S}_1{k}_2/(k_2+k_{-1}))}}((k_3{S}_2+{\lambda }_1)e^{{\lambda }_{1}t}-(k_3{S}_2+{\lambda }_2)e^{{\lambda }_{2}t})\\ {\lambda }_{1/2}=\frac{-k_3{S}_2\pm \sqrt{k_3{S}_2(k_3{S}_2+4k_1{S}_1{k}_2/(k_2+k_{-1}))}}{2}\end{array}$$

(7)

Note that for physical (positive) rate constants and concentrations, λ₁ is always positive and λ₂ is always negative. The first term in Equation (7) has a positive exponent (λ₁) 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 $\raisebox{1ex}{$dP$}\!\left/ \!\raisebox{-1ex}{$dt$}\right.=0$ and cannot be reduced to the second order differential equation; instead, we need to solve a system involving two equations:

$$\frac{dA}{dt}=-k_1{S}_{1}A+(2k_2+k_{-1})I$$

(8)

$$\frac{dI}{dt}=k_1{S}_{1}A-(k_2+k_{-1})I$$

(9)

For the initial conditions, A(0) = A₀ and I(0) = 0; this system has a solution:

$$\begin{array}{l}A(t)=A_0(C_1{e}^{{\lambda }_{1}t}+C_2{e}^{{\lambda }_{2}t})\\ {\lambda }_{1/2}=\frac{-k_1{S}_1-k_2-k_{-1}\pm \sqrt{{(k_1{S}_1+k_2+k_{-1})}^2+4k_{1}S{k}_2}}{2}\end{array}$$

(10)

Both the C₁ and C₂ coefficients are positive and C₁ + C₂ = 1 (for an exact expression of C₁ and C₂, see the Appendix A). The value λ₁ is always positive and λ₂ is always negative for any positive rate constants. Therefore, the value of the decaying term, $A_0{C}_2{e}^{{\lambda }_{2}t}$, will not exceed that of A₀. These calculations indicate that if in an experimental system the product P does not accumulate and A₀ is much smaller than the concentration of the substrate, S₁, 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:

$$A(t)=C_1{e}^{{\lambda }_{1}t}+C_2{e}^{{\lambda }_{2}t}+C_3{e}^{{\lambda }_{3}t}$$

(11)

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₁ and S₂ and might not have time to reach an exponential phase, especially with a high initial concentration of A.

2.2. 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₁, I₁, and P₁ and network 2 consists of A₂, I₂, and P₂; they compete for common substrates S₁ and S₂. Here, the system is defined by:

$$\frac{d{A}_1}{dt}=-k_1{S}_1{A}_1+(k_2+k_{-1})I_1+k_3{S}_2{P}_1-k_0{A}_1$$

(12)

$$\frac{d{I}_1}{dt}=k_1{S}_1{A}_1-(k_2+k_{-1}+k_0)I_1$$

(13)

$$\frac{d{P}_1}{dt}=k_2{I}_1-(k_3{S}_2+k_0)P_1$$

(14)

$$\frac{d{A}_2}{dt}=-k_1^{\prime }{S}_1{A}_2+(k_2^{\prime }+k_{-1}^{\prime })I_2+k_3^{\prime }{S}_2{P}_2-k_0{A}_2$$

(15)

$$\frac{d{I}_2}{dt}=k_1^{\prime }{S}_1{A}_2-(k_2^{\prime }+k_{-1}^{\prime }+k_0)I_2$$

(16)

$$\frac{d{P}_2}{dt}=k_2^{\prime }{I}_1-(k_3^{\prime }{S}_2+k_0)P_2$$

(17)

$$\frac{d{S}_1}{dt}=-(k_1{A}_1+k_1^{\prime }{A}_2+k_0)S_1+k_{-1}{I}_1+k_{-1}^{\prime }{I}_2+k_0{S}_{10}$$

(18)

$$\frac{d{S}_2}{dt}=-(k_3{P}_1+k_3^{\prime }{P}_2+k_0)S_2+k_0{S}_{20},$$

(19)

where S₁₀ and S₂₀ are the concentrations at which S₁ and S₂ are supplied to the reactor.

We should mention at this point that if networks compete only for substrate S₁ and S₂ is considered to be in an unlimited supply (i.e., we consider Equations (12)–(18) with S₂ being constant), autocatalysts A₁ and A₂ 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₂ describes experimental systems where S₂ is in a big excess in relation to S₁ or where conversion of P to A is a monomolecular reaction.

If S₂ 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₁ and A₂ at t = 2000 where, in most cases, the system reaches a steady state. We set all constants to unity, the initial concentrations of A₁ and A₂ to 0.001, k₀ to 0.1, and varied k₁ and k₃, which characterize reactions with substrates. The graph has two characteristic futures: (i) autocatalysts A₁ and A₂ do not coexist, if one has a nonzero concentration and another falls to zero; (ii) competition is not sensitive to k₃ as soon as k₃ 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₁ than in k₃.

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^{\prime }{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^{\prime }$ set to 0.01 (Figure 3b). The graph shows that autocatalysts can coexist and the region of coexistence is higher for low k₃ values. Thus, we should avoid accumulating P in experimental systems with noncatalytic background reactions. We also explored how k₂ and k₋₁ influence the competition between two replicators (Figure 3c,d). The plots demonstrate that k₂, which represents k_cat in a classical Michaelis-Menten scheme, has a stronger influence on the competition than k₃ does, but it is less important than k₁. Interestingly, simultaneously varying k₂ and k₋₁ produces almost a symmetrical plot, which indicates their equal contribution to the competition between replicators (Figure 3d). Overall, the data indicate that k₁, which is responsible for initiating the catalytic cycle, has the greatest effect on the competition, whereas k₃, which is responsible for an axillary reaction that produces an extra molecule in the catalyst, has a minimum effect.

3. 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.