Exponential growth for self-reproduction in a catalytic reaction network: relevance of a minority molecular species and crowdedness

Explanation of exponential growth in self-reproduction is an important step toward elucidation of the origins of life because optimization of the growth potential across rounds of selection is necessary for Darwinian evolution. To produce another copy with approximately the same composition, the exponential growth rates for all components have to be equal. How such balanced growth is achieved, however, is not a trivial question, because this kind of growth requires orchestrated replication of the components in stochastic and nonlinear catalytic reactions. By considering a mutually catalyzing reaction in two- and three-dimensional lattices, as represented by a cellular automaton model, we show that self-reproduction with exponential growth is possible only when the replication and degradation of one molecular species is much slower than those of the others, i.e., when there is a minority molecule. Here, the synergetic effect of molecular discreteness and crowding is necessary to produce the exponential growth. Otherwise, the growth curves show superexponential growth because of nonlinearity of the catalytic reactions or subexponential growth due to replication inhibition by overcrowding of molecules. Our study emphasizes that the minority molecular species in a catalytic reaction network is necessary to acquire evolvability at the primitive stage of life.


Introduction
Self-reproduction is a crucial step for the transitions from nonliving to living matter [1,2,3]. External resources are converted into internal components, the abundance of internal components increases, and they reproduce themselves as a set of components. How such reproduction of components is achieved has been the focus of research interest for decades.
Here, for reproduction, amounts of all components should increase roughly at the same rate. Then achievement of exponential growth is essential for stable reproduction. First, if the amounts of all components increase as a set while they help each other to replicate, then their replication is proportional to their amounts, so that the exponential growth is expected. Second, because of the attainment of exponential growth, Darwinian selection is achieved across different units of growth (e.g., protocells) [4,5]. Third, during exponential growth, the amount of each component i increases as x i (0) exp(λ i t) with time t, so that if growth rate λ i is identical across all components i, then initial composition {x i (0)} is preserved across components i. By contrast, in other types of growth, preservation of the initial concentration would be difficult.
Hence, to achieve self-reproduction while maintaining the concentration, molecular mechanisms of replication of each component in the entity are not sufficient for selfreplication. The replication of each component should be orchestrated in a controlled manner. In other words, to produce a copy of a cell with nearly the same composition, yields of the components should be proportional to the amount of the composition in this entity, with the same proportion coefficient.
Note that biological reactions are facilitated by catalysts, whereas all the catalysts are synthesized as a result of catalytic reactions. If a set of catalysts is synthesized as a result of these catalytic reactions and the growth rate of each component is proportional to the concentration of the catalyst, a step toward exponential growth may be made.
Indeed, it has long been said that catalytic reactions are essential for reproduction of cells at the primitive stage of life [2,6,7,8,9,10]. In particular, a hypercycle model provides a basic scheme of catalytic reactions, in which different molecular species catalyze the replication of each other, to prevent functional molecules from being lost by inevitable replication errors [4,11,12,13,14,15]. Dynamics of catalyst concentrations, however, are inherently nonlinear, and consequently, it is a nontrivial question whether the exponential growth synchronized over all components is achieved in a mutually catalytic reaction network.
Besides the nonlinearity in chemical reaction dynamics, smallness of the molecule number is often important. When the number of each species of molecule is large, the rate equations for continuous concentration variables will be applicable to study the reaction dynamics. At the primitive stage of cells and in modern cells, however, the number of some molecules is sometimes small, so that the discreteness in numbers 0, 1, 2,... has to be considered seriously [16,17,18,19]. Furthermore, each molecule for replication is a large polymer, and when such molecules are replicated, they soon reach crowded conditions [20,21]. As replicated molecules reach crowdedness, the standard reaction-diffusion equations may not be valid, where large excluded volumes can introduce non-negligible effects. Nevertheless, these two factors, discreteness and crowdedness in biochemical reactions (which are also important in modern cells), have not been seriously taken into account in the context of the research on the origin of life or protocells.
Recently, we considered a hypercycle model with two mutually catalyzing chemicals, to demonstrate that reproduction of a protocell via a growth-division process emerges when the replication and degradation rates of one chemical are slower than those of the other, and such a chemical is a minority, thus leading to discreteness, whereas replication causes molecular crowding [22]. In addition, the protocell divides after the minority molecule is replicated at a slow synthesis rate. Thus, synchrony between cellular divisions and molecular replication is achieved. The results were first demonstrated by means of Langevin dynamics in a continuous space [22] and later via a cellular automaton model in a square lattice [23] to clarify the effect of excluded volumes.
Here, we clarify the relevance of minority molecules to exponential growth of the system. By simulating the cellular automaton model in two-(2D) and threedimensional(3D) square lattices, we show that exponential growth occurs naturally because of the minority molecule, i.e., when the reaction rate of one chemical is considerably slower than that of the other, and replication causes molecular crowding. When the reaction rates of both chemicals are slow, in contrast, superexponential growth in the molecule number is observed because of the nonlinearity of catalytic reactions, whereas when both reaction rates are fast enough, subexponential growth is obtained because of the inhibition by overcrowding of molecules. In addition, we compare our data with the results of the reaction-diffusion model to demonstrate that discreteness is essential for the exponential growth, to achieve orchestrated self-replication of each molecular species. This paper is organized as follows. In section 2, we briefly review three types of growth laws and their consequences for selection. In section 3, the cellular automaton model is introduced, incorporating the mutually catalyzing reaction of molecules. On the basis of simulation of the model, we classify the growth curves in section 4 to show that the minority status of one molecular species and crowdedness indeed lead to exponential growth. In section 5, we summarize and discuss our results.

The simple growth law and its selection consequences
Here, we discuss simple model equations for exponential, sub-, and superexponential growth and show that Darwinian selection holds only during exponential growth, according to the arguments presented by Eigen and Schuster [4].

Exponential growth and survival of the fittest
Simple exponential growth assumes that the process of reproduction is characterized by the stoichiometric equation where A indicates the unit of growth such as replication of a molecule, S denotes the consumed resource, and k is the growth rate constant. The rate equation for this process is written as where x and s are the density of A and S, respectively, andẋ represents the time derivative of x.
When the density of S is assumed to be stationary and is incorporated into k, the solution of eq. (2) is well known: x(t) = x(0) exp(kt). When two populations with different rate constants compete, the one with the higher growth rate wins the competition, and the opponent is competitively excluded.
To simply illustrate this situation, we consider the growth dynamics of two populations (labeled as 1 and 2) in a flow reactor aṡ where φ = k 1 x 1 +k 2 x 2 and x 1 and x 2 are the densities of population 1 and 2, respectively.

Subexponential growth and survival of everyone
The pioneering studies on template-directed replication [24,25,26,27] are known to comply with a parabolic growth law (owing to product inhibition) and lead automatically to coexistence under conditions allowing for subexponential growth without competitive exclusion. The growth of a template follows the rate equatioṅ where l is the rate constant for spontaneous template formation, and 0 ≤ p < 1. To show that the subexponential growth results in coexistence, we consider the growth dynamics in a flow reactor aṡ Here, x 1 + x 2 = 1. When one of the species is dominant and the other is rare, we can show that the rare species can grow in number. By writing densities of the dominant and rare species as x 1 = 1 − δ and x 2 = δ, for small δ, we geṫ Accordingly,ẋ 2 > 0 holds for sufficiently small δ. This means that the rare species can always invade the dominant species. On the other hand, the solution of eq. (6) is given by It can be shown straightforwardly that this fixed point x * 1 , x * 2 is stable against perturbations by linear stability analysis. Hence, the subexponential growth implies coexistence, i.e., survival of everyone.

Superexponential growth and survival of the first
When two or more individuals are necessary to replicate another molecule, the simplest scheme accounting for the replication is where n ≥ 2. Assuming again that the density of resource S is stationary and incorporated into rate constant k, then for n = 2, the growth equation iṡ and the solution reads approaching infinity as t → 1/[x(0)k]. The importance of hyperbolic growth lies in the consequences for selection. When two competitors with different k values grow together, the one with greater x(0)k will "blow up" first. This means that the outcome of selection depends on the initial condition x(0), and thus, hyperbolic growth implies survival of the first comer, i.e., the advantage of greater x(0). The same results are obtained for the growth equationẋ = kx n with any exponent n > 1.
Here, we can again adopt the competitive growth dynamics in a flow reactor aṡ where φ = k 1 x n 1 + k 2 x n 2 . In this case, stability of the solution x * i = 1 and x * j = 0 for others is given byδ = (nk j x * j − 1)δ = −δ for the perturbation x j = 0 + δ and x i = 1 − δ. Hence, both fixed points x * 1 = 1 and x * 2 = 1 are stable. In other words, if one species is dominant initially, the other species cannot invade even if its k value is greater than that, i.e., survival of the first.
Although we illustrated the behaviors by the case of two species, it is quite straightforward to show that for any number of species, the exponential growth, i.e., x ∝ x leads to survival of the fittest, subexponential growth, i.e.,ẋ ∝ x p with p < 1 leads to the survival (coexistence) of everyone, and superexponential growth, i.e.,ẋ ∝ x n with n > 1 leads to survival of the first.

Model
Now, to study a stage in the origins of life, we consider mutually catalyzing molecules that replicate themselves. As the simplest case, we assume that two molecular species X and Y mutually catalyze the replication of each other as X + Y → 2X + Y and Y + X → 2Y + X. According to the discussion in section 2, the dynamics of the reactions obeying the simple model equations essentially follow the superexponential growth in which two individuals are necessary to replicate each other. In this case, the concentrations of X and Y generally show hyperbolic growth, and exponential growth is not attained.
Here, the simple model equations ignore the spatial distribution of molecules and imply that the system is homogeneous. In reality, however, molecules undergo reactions in space, and it is not obvious whether the assumption is always valid, i.e., spatial structure may change the growth dynamics. In addition, catalytic molecules are generally giant polymers; therefore, effects of crowding may be relevant due to the excluded volume.
Accordingly, we investigate the dynamics with effects of spatial distribution and excluded volumes of molecules to identify a possible mechanism underlying the exponential growth. Here, we employ a simple cellular automaton model to ensure the two effects, while we also analyze the conventional reaction-diffusion model without the effects of discreteness in molecule number and crowding of molecules in the Appendix.
We consider square lattices in 2D space with dimensions L x × L y , and 3D space with dimensions L x × L y × L z and periodic boundary conditions are imposed. Each site is empty or occupied by one molecule to ensure appropriate representation of the excluded volume of molecules. The species identity is represented by a "color" of the molecule, and replication of molecules proceeds based on the catalytic relation between the species as described below. We use discrete simulation steps and update the system at each step by applying three processes to each molecule in a random order.
The first process is molecular replication, as outlined in figures 1(a) for 2D space and 2(a) for 3D space. In this paper, we consider the simplest hypercycle, in which two molecular species X and Y mutually catalyze the replication of each other as When a molecule is located at a site neighboring its catalyzing molecule, replication occurs with a given probability, p, and a new molecule is added at a randomly chosen site (among the six sites for 2D space and ten sites for 3D space) neighboring the reaction pair. The examples in figures 1(a) and 2(a) illustrate a case where replication of a green molecule (X) is catalyzed by a red molecule (Y). If the chosen neighboring site is occupied, then the molecule is not replicated. The added new molecule becomes X or Y with probabilities γ X and γ Y = 1 − γ X , respectively. Note that in the above reaction, some substrate (or resource) chemicals S X , S Y are presumed, as in the reactions X + Y + S X → 2X + Y and Y + X + S Y → 2Y + X. Here, these substrate molecules are assumed to be supplied sufficiently.
The second process is molecular degradation. At each step, degradation takes place via removal of each molecule at a given probability representing the degradation rate of the system (figures 1(b) and 2(b)). Degradation proceeds as X → φ and Y → φ with probabilities a X and a Y , respectively.
The third process is diffusion. At each step, every molecule moves to one of the nearest neighboring sites if the destination site is empty [figures 1(c) and 2(c)]. Hence the time steps and lattice sizes are chosen so that random walk takes place at each step. Degradation: Each molecule is removed from the system with fixed probabilities a X and a Y for molecular species X and Y , respectively. (c) Diffusion: every molecule moves to one of its four nearest sites (highlighted in gray) if the selected site is empty.
The destination site is chosen randomly from the neighboring sites, or the molecule remains at the original site if the site is occupied.

The growth curve for the symmetric case
First, we show the growth curve for 2D space in the case where the probabilities of replication of X and Y are equal, i.e., γ X = γ Y . Simulations were carried out from the initial condition where the molecules are located at coordinates L ini × L ini , and the molecules randomly assigned to X or Y with equal probabilities. Here, the values of L ini are set to 10. Figure 3 illustrates how the numbers of molecules X and Y increase with simulation steps at different p values. Here, we set a X = a Y = 0. When p is small, the growth curves are downward convex on a log-log scale and apparently "blow up" at some point If molecules X and Y are located next to each other, then replication can occur with probability p, and a new molecule is added at one of its ten neighboring sites if the site is empty. (b) Degradation: Each molecule is removed from the system with fixed probabilities a X and a Y for molecular species X and Y , respectively. (c) Diffusion: every molecule moves to one of its six nearest sites if the site is empty.
where all the sites of the system are occupied. As shown in section 2.3, when two individuals are necessary to replicate another molecule, the growth curve is written as eq. (10), i.e., superexponential growth. The equation can be applied to the present case by assuming that the numbers of X and Y are equal. To confirm that the solution of eq. (10) is suitable for small p, in figure 4, we present evolution of the inverse of the numbers of X and Y . The equation indicates that the inverse decreases with the function −kt + C, in agreement with the numerical results. In addition, values of k are proportional to those of p in the simulation.
As p increases, the curves in figure 3 gradually shift to a power law and approximately follow t 2 . This slowing down of the shape of curves indicates the effect of spatial structure. When p is large, the replication of molecules is fast enough so that the molecules maintain a single cluster, and replication occurs mainly at the periphery of the cluster (figure 5). When the replication takes place mostly at the periphery, the increase in molecule numbers at each step is proportional to x  of molecules, and d denotes spatial dimensionality. Then the increase in the molecule number ∆x is expressed as ∆x = x 1/2 and x 2/3 at d = 2 and 3, and the solution of the equation leads to x(t) ≈ t 2 and t 3 , respectively. Actually, such power law behavior is observed for the 3D model as shown in figure 6, which fits power law t 3 when p is large.
The dependence on degradation rate a = a X = a Y is plotted in figure 7, where the growth in the molecule number occurs as degradation is smaller than a given threshold. When the growth occurs, the growth curve essentially conforms to the power law at large p (figure 7a) and is superexponential when p is small (figure 7b), independently of the degradation rates. This result suggests that attainment of exponential growth in symmetric cases, even if possible, requires fine-tuning of the parameters and is not practically feasible.
Qualitatively similar behavior is also observed in a reaction-diffusion system with site capacity as shown in the Appendix.

Exponential growth
As explained in the previous section, the growth curves follow either the subexponential function where replication is inhibited by overcrowding of molecules, or the superexponential function owing to the nonlinearity in the catalytic reactions. Then a question arises, How is self-replication via exponential growth possible?
Here, we find that exponential growth is possible when the replication and degradation rates of one species, say Y , are much slower than those of the other species, X. As for the parameter dependence, our previous study showed that the system  In the Extinction region at large a Y , replication of Y is slower than its degradation, and therefore all the molecules degrade and become extinct. In the Explosion region at large p Y = pγ Y , molecules X and Y increase in number in a mixture of X and Y . In the Division region, the numbers of molecules X and Y increase with the spatial structures in which each slowly replicating molecule Y is surrounded by a group of fast-replicating  molecules X. To see how the number of molecules increases, we demonstrate evolution of the number of X and Y in figure 8 against different values of parameters p Y for a given value of a Y as in the dotted box. Here, the initial condition for the simulation is a single Y located within a group of X molecules with dimensions L ini × L ini , where the values of L ini are fixed at 10. For the parameters corresponding to the Division region, the numbers of both molecules X and Y increase linearly on the semi-log scale [see the right-hand panels of figure 8(b)], implying the exponential growth of molecule numbers. In addition, the growth constant in the exponential function is approximately the same for both molecules X and Y , which is proportional to p Y = pγ Y .
For the values of p at which the localized structure occurs, the replication of molecules is fast enough as compared with the rate of diffusion. Then, as p Y = pγ Y increases to the Explosion region, the replication of molecules is inhibited by the crowded conditions, and the growth curves are upward convex, indicating subexponential growth [left-hand panels in figure 8(b)].
The exponential growth via the minority molecule is also implemented in the 3D model when γ Y ≪ γ X (figure 9), where the number of molecules grows via formation of localized structures ( figure 10). Note, however, that the exponential growth is not observed in the reaction-diffusion approach (see Appendix). This finding suggests that exponential growth is achieved due to the synergetic effect of discreteness and crowdedness of molecules.

Summary and discussion
In summary, we demonstrated that self-replication via exponential growth occurs only when the rates of replication and degradation of one molecular species are much slower than those of the other in 2D and 3D cellular automaton models of the hypercycle. Herein, the synergetic effect of molecular discreteness and crowding is necessary to ensure the exponential growth. In contrast, if the reaction rates are not so different between the chemicals, then the growth curves reflect superexponential growth because of the nonlinearity of catalytic reactions, or subexponential growth owing to replication inhibition by overcrowding of molecules.
As emphasized in section 1, exponential growth is essential for orchestrated replication of all the components and for Darwinian selection of self-reproducing units. Our results suggest that the minority molecules with slower replication and degradation rates represent a possible mechanism behind exponential growth, by naturally preventing overcrowding of molecules. Therefore, for a balanced yield of components even in stochastic chemical reactions, the replication of each component is orchestrated in a controlled manner.
In biological systems, polymer replication is commonly adopted as a basic process, while catalytic polymer molecules are usually gigantic. Then, as replication reactions progress, crowding becomes inevitable and typical in live cells. Achievement of such a high local concentration is sometimes advantageous because effective concentration of macromolecules increases, and therefore reactions can occur more frequently. On the other hand, overcrowding may inhibit reactions because it can deplete resources coming from the external environment. The exponential growth controlled by a slowly replicating minority molecular species can take advantage of molecular crowding and simultaneously prevent overcrowding.
Although the present study does not include a selection process, our results indicate that the localized structure with minority molecules may be optimized after several rounds of selection during the exponential growth, because this structure can work as a unit for Darwinian selection. In our previous study [28], by presuming a compartment and its division, i.e., exponential growth at the cell level, it was suggested that a state involving minority molecules is relevant to enhanced evolvability because changes   Red and green dots denote molecules Y and X, respectively. For visual clarity of the presentation, radii of X and Y are shown 0.3 shorter and twice longer, respectively, than those of actual sizes.
in minority molecules influence catalytic activity of the protocell more strongly than those in majority molecules. Our results provide further insights into the advantages of minority molecules for achieving evolvability because exponential growth, being essential for Darwinian selection, is naturally shaped only by a minority molecule. One of the remaining issues is to determine whether such minority molecules with slow dynamics emerge at the primitive stage of the catalytic reaction network for reproduction. Note that without such molecules, superexponential growth becomes possible, and faster growth may be plausible at first glance. This superexponential growth, however, may be replaced by subexponential (power law) growth owing to overcrowding, as molecules with a higher replication speed emerge, and hence this kind of growth cannot be sustained. Furthermore, the explosion in the molecule number in the superexponential-growth case may cause a serious lack of resource molecules for the replication, and the growth is likely to stop. (Note that in the simple models that we simulated here, resources are supplied very rapidly, but in reality, there may be resource shortages as molecules are replicated).
In contrast, during the exponential growth orchestrated by a minority molecule, a spatially localized structure emerges that alleviates the lack in resources, and consequently, the growth is sustained. The localized structure works as a unit for selection, and with the mutational change in the minority molecule, its catalytic activity may increase, which may further enhance the difference in replication speeds between molecules X and Y in the model, thereby stabilizing the difference in the rate of replication. This reaction rate separation, on the other hand, will enhance the evolvability because the change in such a slow minority molecule will influence the synthesis of other components, leading to a directional change [29,30].
As the growth continues, resources for replication of the molecules may become limited even if the localized structure is formed around the minority molecule. Such resource limitation will lead to diversification of components, as recently demonstrated on models of the protocell [31,32]. Such diversification may yield a localized structure of diverse components synthesized with the aid of a minority molecule, as seen in modern cells.

Acknowledgments
This research was partially supported by a Grant-in-Aid for Scientific Research (S) (15H05746) from the Japan Society for the Promotion of Science (JSPS). The authors declare that they have no conflicts of interest.

Appendix A. The reaction-diffusion model
To confirm that discreteness of molecules is essential for exponential growth shown in the main text, we investigated a 2D reaction-diffusion model. We consider square lattices in 2D space with dimensions 100 × 100 and a periodic boundary condition. For each site at position (i, j)(1 ≤ i, j ≤ 100), we define densities of X and Y respectively denoted by C ij X and C ij Y , and they follow the reaction-diffusion equations, Here, in the right-hand side of each equation, the first and second terms represent replication and degradation of the molecules, and the third term denotes diffusion to the neighboring sites. In the replication term, we introduce the capacity of a site, C max , so that the replication is inhibited as total densities of X and Y approach the capacity.
An initial condition is imposed such that the densities of X and Y are respectively fixed at 0.1 for a small region (10 × 10) and for the other sites fixed at zero.
For the symmetric case, γ X = γ Y , the growth curve is shown in figure A1 for different p values. The curves are qualitatively similar to those in the cellular automaton model: At small p, the number of molecules blows up at some point according to function 1/(A−Bpt), and at large p, the increase follows a power law. For large p, the evolution of spatial distribution of X and Y shows behavior similar to that of the cellular automaton model ( figure A2).
On the other hand, during the growth in the asymmetric case, γ X ≫ γ Y , the increases in X and Y numbers do not show synchronous growth with the exponential curve as in the cellular automaton model, and spatial distributions of X and Y are rather homogeneous (figure A4).   Figure A3: Evolution of the amounts of X and Y in the 2D reaction-diffusion model for the asymmetric case γ X > γ Y . Here, C max = 10, p = 1, a X = 0.01, and a Y = 0.