Excess Mutual Catalysis Is Required for Effective Evolvability

The fundamental question of how primitive life emerged on the prebiotic Earth has drawn considerable scientific attention throughout the centuries [2, 5, 14, 15, 22, 42, 59, 64]. The path from organic mixtures (i.e., the primeval soup) to reproducing lifelike protocells has traditionally been dominated by two different views: the genetic, or replicator-first, approach, and the metabolism-first approach [2, 42]. Both acknowledge the need for reliable information storage and transfer, assisted by self-replication. The replicator-first approach suggests that life began with a single self-perpetuating biopolymer (e.g., RNA) [14, 15, 18, 19, 22, 37, 42, 64], which later evolved into multimolecular networks under the replicator's control. Orgel [41] has highlighted the relationship between molecular replication and the concept of autocatalysis or self-catalysis. The metabolism-first approach suggests that the very first life precursors must have been relatively complex molecular networks arising via spontaneous accretion of simpler organic molecules [3, 9, 24, 25, 34, 48, 51, 53, 60]. In this scenario, it is further proposed that faithful reproduction directly stems from certain network attributes. Therefore, one should better understand the network properties of the implicated molecular assemblies [1, 47, 57, 66] if one can merge the two seemingly conflicting scenarios for prebiotic evolution.

One embodiment of the metabolism-first view is the lipid world scenario, which considers non-covalent assemblies of amphiphiles, such as micelles and vesicles formed by lipids [8, 39, 48, 50, 53, 69]. These are assumed to store information in the form of nonrandom molecular compositions, and pass it to progeny via homeostatic growth accompanied by fission [49]. The graded autocatalysis replication domain (GARD) kinetic model for prebiotic evolution quantitatively describes the details of such a process. It elaborates some of its evolution-related attributes [10–12, 27, 30, 44, 50, 58, 62, 67, 68], with an implied route to minimal protocells [8, 45, 49, 63, 65]. The model is based on a catalytic network, usually presented in the form of a matrix β with autocatalysis (self-catalysis) and cross (mutual) catalysis terms. Importantly, the system is kept away from thermodynamic equilibrium by assembly fission [49]. Key in GARD dynamics are compotypes—clusters of replication-prone quasi-stationary states (composomes, a term derived from the notion of compositional genomes [49]), proposed to play a crucial role in the GARD's evolutionary behavior. Introducing substantial inhibition in β is expected to result in net catalysis because an inhibitor of an inhibitor is an activator [20].

Catalysis, the enhancement of reaction rate by an external chemical component, was recognized as early as 1836 by Berzelius, and Ostwald applied the term autocatalysis in 1890 to reactions that gain speed as they proceed [26, 44]. In the genetic approach to life's origin, researchers invoke one or several autocatalytic molecules as the core of a prebiotic entity. This is exemplified by the hypercycle, a set of self-replicating polynucleotides, coding for and acted upon by enzymes [10, 30, 58]. In the metabolism-first domain, autopoiesis [67] and the chemoton model [12] are examples of collective autocatalysis [25].

Collectively autocatalytic systems feature a central role not only for self-catalysis, but also for mutual catalysis. In this, they arguably resemble present-day living cells, which harbor self-catalytic polynucleotides as well as a plethora of mutual catalysts that constitute metabolic pathways. Here we utilize a metabolism-first simulator to examine the relative importance of the two catalytic modes (self- and mutual catalysis). Previously [11], an abstract chemistry model has been used to demonstrate that self-maintaining organizations arise only once self-catalysis is completely inhibited [11, 62]. We attempt to extend such results in the realm of the GARD kinetic model, asking what features of the β network contribute to the evolution of the ensuing compositional assemblies. It is shown that excess mutual catalysis is a necessary, though not sufficient, condition for displaying several evolutionlike characteristics, including a high number of composome types, higher evolvability scores, and a significant response to selection.

Recently, it has been argued that collectively autocatalytic metabolic networks, such as the GARD, do not allow for fitter compositional genomes to be maintained by selection. Vasas et al. [68] compared the frequency ranking of random GARD compositional assemblies before and after selection, and found that the relative ranks changed only slightly. This was taken as evidence for an inherent evolutionary limitation of metabolism-first scenarios. Here it is demonstrated, based on a large number of simulations, that when quasi-stationary composomes rather than arbitrary compositions serve as selection targets, GARD networks are capable of a significant response to selection. Importantly, this can happen chiefly when a high proportion of mutual catalysis is present in a GARD network. The results highlight the potentially important role of mutual catalysis, as compared to self-catalysis, in the emergence of early lifelike systems.

The matrix β represents a network of self-catalytic (diagonal elements) and mutually catalytic (off-diagonal elements) catalytic interactions (Figure 1), with self-catalysis represented by the case j = i ( Appendix A.1, Equation 13). The matrix elements are randomly drawn from a lognormal distribution ( Appendix A.1 and Equation 12) [49].

The model is subjected to a kinetic Monte Carlo simulation based on Gillespie's algorithm [16, 17, 51] using parameter values similar to those employed in previous studies (Table 1). Simulations are run using MATLAB versions 7.6–7.10 (the GARD10 code is available upon request). A set of 10,000 GARD simulations is generated, all with the same parameters, and each with a different matrix β generated by the MATLAB pseudorandom number generator with seeds 1–10,000. The validity of the conclusions drawn here is ascertained by repeating the simulations with smaller data sets, with seeds 1–2,000 and 2,001–4,000, striving to verify that the entire 10,000-strong data set adequately represents the GARD simulation space. The random sampling of β values may be perceived as representing different possible GARD environmental chemistries.

All the compositions belonging to composomes in the entire simulation undergo k-means clustering [56, 61], and the centers of mass of the resulting clusters are defined as compotypes.

An initially random population of a fixed number of assemblies, L_pop, is allowed to simultaneously grow according to Equations 1 and 2 and its idiosyncratic composition. When one of the assemblies reaches the limiting size N_max, it divides by random fission, and a randomly chosen assembly from among the other L_pop − 1 assemblies is removed, thus keeping the population size constant. This is repeated for GEN splits (Table 1). This protocol is based on the classical Moran process [36, 68, 70], and to some degree reflects an earlier attempt to simulate GARD populations [38].

The frequency of the target in each population is defined as the number of assemblies that are highly similar (H ≥ 0.9) to the target compotype taken from regular GARD for the same β network (Figure 13 in  Appendix A.4). Selection is exerted by performing a simulation with the same parameters, biasing the growth of assemblies toward a target compotype as for regular GARD (Equations 8 and 9). The selection excess is defined as in Equation 10, where f_T− and f_T are respectively the fractions of assemblies within the population belonging to the target compotype before and after selection.

We used GARD simulations to ask what is the selection response of compositional assemblies. A value for the selection excess was obtained for each of 10,000 simulations, using a modest value of the fitness gain, s = 1.1, in line with previous work [68]. Figure 2a shows the correlation between the frequencies of the target compotype with and without selection (examples of regular GARD carpets before and after selection are given in Figure 14 in  Appendix A.4). An overall skew is seen here toward positive selection. The figure also demonstrates that significant positive selection, as well as negative, occurs over most of the range of f_T.

Figure 2b shows the distribution of selection excess values for the entire data set (Equation 10). Importantly, a considerable percentage of the simulations (33%) show positive selection, with a mean selection excess of 1.38 for selection excess >1.05, and as much as 10% shows selection excess >1.5. Interestingly, 31% of the cases showed negative selection, with a mean selection excess of 0.775 for selection excess <0.95, and about 36% were neutral to the selection pressure. Similar to the skewness in Figure 2a, there is a slight bias in favor of positive selection, as indicated by an overall mean selection excess equal to 1.05. Notably, higher mean selection values positively correlate with the number of other compotypes coexisting with the target compotype in a given system (Figure 15 in  Appendix A.5).

GARD simulations are used to see how attributes of the catalytic network embodied in the matrix β govern the evolution-related dynamics of compositional assemblies. It is asked how the mutually catalytic power p_mc (Equation 4) influences the selection response. A clear trend appears here, whereby strong positive or negative selection is found almost entirely for p_mc higher than 1 (Figure 3b).

The main trends appear also at lower simulation counts, barring small-number fluctuations at high p_mc (Figure 3a). For example, for the range of p_mc > 100, a meaningful p-value with 5% significance level is achieved only after performing more than 2,500 simulations (Table 3 in  Appendix A.6). The other two evolution-related parameters withstand similar scrutiny (below).

The foregoing simulations of the regular GARD model addressed the case in which at each time point only one GARD assembly is considered. To enhance the capacity to draw conclusions about selection in GARD, 1,000 simulations were performed, each for a population of 1,000 assemblies, under the constant population conditions. Figure 4 shows an example of the dynamics for one of the networks. Starting from a population of random assemblies, the population frequency of the target compotype gradually grows over the first 10,000 split events, reaching a plateau with fluctuations, signifying the compositional preference imposed by the matrix β towards this compotype. When selection toward this compotype is applied (Equation 9), this general behavior is retained, with a faster growth and a higher plateau, that is, showing positive response to selection. Similar to Figure 3b, strong positive or negative selection is much more prevalent for p_mc values higher than 1 (Figure 5b).

The effect of selection pressure on the frequency of the target compotype for all 1,000 networks is presented in Figure 5a. Similarly to Figure 2a, an overall skew toward positive selection is seen (about 50% of cases), with some cases of negative (about 15% of cases) or no response to selection, and with a mean selection excess of 1.254 ± 0.804. Significantly, the ratio of the number of simulations showing positive selection to that showing negative selection increased more than threefold, from 1.06 in the regular GARD to 3.33 in the population-GARD. In line with previous work [68], the growth bonus was calculated when the assembly size was N_min (Equation 8). When the bonus was calculated for all time points between N_min and N_max (for a smaller set of 100 population-GARD simulations), the overall selection response seems to become even more positive (70% of cases), with a higher selection excess value of 1.399 ± 0.997.

The influence of p_mc on one of the attributes of GARD diversity, the mean number of different compotypes appearing in a simulation, is now analyzed. It is found that as p_mc increases, so does the mean number of compotypes, reaching a maximal value of nearly 3 at p_mc = 100 (Figure 6a). Furthermore, in the realm of excess self-catalysis (p_mc < 0.5), one compotype appears in an overwhelming majority of the cases (91%) (Figure 6b). In contrast, compotype counts between 2 and 6 are almost entirely confined to the domain of excess mutual catalysis (p_mc > 2). Curiously, even among the ∼5,000 simulations that show only one compotype, a large majority have p_mc > 2, suggesting that high mutual catalysis is a necessary but not sufficient condition for a high number of compotypes.

The similarity autocorrelation function (Equation 6) and its derived parameters (Equation 7) are employed to obtain information on the evolutionlike dynamics of GARD assemblies. One possible interpretation of the value of τ is a depiction of the whole-simulation average of the assembly compositional lifetime. Longer τ may be taken to represent better average maintenance of compositional similarity between consecutive GARD generations, symbolizing better reproduction fidelity. Likewise, 1/τ may be thought of as related to the compositional mutation rate. Indeed, effective compositional preservation is implicated by the most frequent number of generations, τ ≈ 3, with a non-negligible probability for τ ≥ 10 (Figure 7a). Note that τ does not represent the composomal lifetime. In fact, the most probable target compotype lifetime (taking for simplicity the maximal time from each simulation) is 30, and the average is 434 generations (Figure 7c). The other similarity autocorrelation parameter, H₀, is interpreted here as showing the residual compositional similarity among assemblies along many generations in the entire simulation. Thus, 1 − H₀ is taken as proportional to the overall compositional diversity of assemblies across the entire simulation. Note that H₀ is not strongly correlated with the compotype count (Figure 16 in  Appendix A.5, correlation coefficient −0.049, r² = 0.89) and therefore constitutes a rather independent diversity assessment attribute. The most probable H₀ value is ∼0.5, with a smaller probability peak at H₀ ≈ 1. The latter stems from simulations in which a single compotype tends to dominate.

The effect of the assembly pre-fission size N_max on GARD's evolutionlike behavior was studied by performing two additional sets of 10,000 simulations, each with the same parameters as in Table 1 except for N_max = N_G/2 and 2N_G (Figure 9; Figure 17 and Table 3 in  Appendix A.6). While for the smaller N_max value, the percentage of beneficial outcomes seems to be even higher than for the nominal N_max = N_G, a larger N_max value appears to have a disruptive effect because the system is nearing the equilibrium steady state [23]. This is especially seen for the compotype count and the evolvability score.

One of the dominant concepts in prebiotic evolution research is the replicator-first scenario [10, 32, 40]. Based on the concept that molecular replication is related to self-catalysis [41], such views may be perceived as related to the RNA-first scenario, positing that life began with a unique self-replicating polyribonucleotide. In this realm, it is argued that more complex interaction networks have arisen only at later stages, as when precursors for the autocatalytic molecule have been exhausted [31]. Our simulation results demonstrate an advantage for a network-first scenario, in which a large number of molecular components mutually interact. While arising from a metabolism-related framework, such results may be taken as relevant to the question of whether life's early precursors were a set of replicators or a metabolic network. Note that the present work makes a direct comparison between a metabolic network with frequent self-catalytic interactions and a metabolic network with frequent mutually catalytic interactions, and therefore has only indirect relevance to the question of the validity of replicator models. It is conceivable that future work incorporating templating biopolymers together with mutually catalytic networks will better resolve this issue.

A widespread argument against metabolism-like entities being the first seed of life is the assertion that metabolic networks cannot store and propagate information. The GARD model may be viewed as a counterexample, as it is endowed with a (limited) capacity to store and propagate compositional information. This has implications for a set of previously proposed models involving networks of molecular interactions. Two of the earliest relevant concepts are Gánti's chemoton [12, 13, 63] and Maturana and Varella's autopoietic systems [35, 67]. Autopoiesis characterizes a spatially confined network of molecular components, whose mutual interactions continuously regenerate the network itself. The chemoton is described as a system of three subnetworks: metabolite generation, template copying, and membrane synthesis. We prudently suggest that GARD may be viewed as a special case of autopoietic-chemoton-like models, where template copying and compartmentation are embodied in one entity, and a continuous supply of metabolites is afforded by the spontaneous accretion of lipids from the buffered environment.

An important result of this work is that networks within a certain range of kinetic parameters, namely those that exhibit excess mutual catalysis, lead to enhanced diversity and evolvability of GARD compotypes. The compotype count is a direct indication of the degree of composomal diversity. This result is related to an important aspect of early evolution: Self-catalysts tend to propagate their own identity and suppress processes essential for the increasing complexity necessary for transitions from early seeds of life toward systems resembling present-day life. The presently demonstrated importance of mutual catalysis echoes the notion of systems prebiology [21, 57], whereby it is suggested that life began its trajectory from complex chemical mixtures obeying network behavior similar to that of metabolism in present-day cells.

One of the unique corollaries of the GARD model is the emergence of composomes, dynamic states of compositional assemblies that embody both metabolism-like characteristics and a rudimentary capacity to store and propagate molecular information [49]. Composomes may be considered as forming bridges between seemingly disparate views of the early seeds of life: metabolism first and replicators first. Compotypes are further defined as centers of mass of composome clusters, which may be regarded as analogous to species or quasi species [6]. This is due to the fact that a compotype is a distinct entity, with distinct physical properties and hence fitness encoded in its compositional information, different from those of other compotypes but still harboring considerable internal variability of constituents. Therefore, compotypes are considered as natural targets of selection, as compared to randomly chosen compositions, as previously pursued [68]. Note that here we have a measure of selection inherently present in the GARD model even in the absence of external selective pressure, due to the fact that different composomes have different average growth rates. This is seen in the present population GARD simulations, which are seeded with a random population, but show a gradual increase of the population frequency of a specific compotype even in the absence of externally imposed selection. This increase comes at the expense of other compositions because of the constant population condition.

The present results show that GARD assemblies can exhibit positive or negative selection toward a compotype target, as well as no selection at all. While in regular GARD the overall average selection excess is merely 1.05, it is noteworthy that as many as 10% of the simulations show high selection excess, >1.5. Importantly, these general results are borne out both in simulations of the regular model and in simulations involving populations of assemblies. Previously, GARD population dynamics has been studied by addressing various emergent properties, including a comparison of finite and infinite chemical environments [38]. Another study [70] showed that compositional inheritance also emerges in the GARD model variants involving assembly populations and spatial proximity interaction effects, and that it emerges in both a thermodynamic and a kinetic interaction regimen.

Analyzing GARD, both positive and negative selections can be observed in practice only when the underlying network exhibits mutual catalysis excess. This conclusion is strengthened by its demonstration in two different simulation modes: in the regular model and in populations. Notably, positive selection is observed appreciably more often in population GARD simulations, perhaps reflecting the advantage of addressing populations of competing entities with different reproductive rates. Furthermore, this selection response tends to be augmented as the number of coexisting compotypes increases in a given simulation, which may indicate a capacity of selective forces to provide an edge to the target compotype in inter-compotype competition. Further in-depth analyses (currently underway) of the ultrastructure of the β network, as well as subnetworks (quasi compartments [68]), could lead to a better understanding of the influence of p_mc and the compotype count on selection.

The present method for biasing the growth rate of a GARD target composition is in principle similar to that used previously [68]. In both cases, modifications are in effect introduced to β matrix elements. However, the previous analysis utilizes an interim formalism, the Eigen equation, for replication-mutation dynamics [10], and the selection-related modification is exerted by multiplying the growth rate by fH, defined in the same way as in Equation 8. The method utilized here involves direct modification (Equation 9), a possible explanation for the discrepant results obtained by the two reports. There are, however, additional significant differences between the two studies: (a) a pre-fission value N_max = 100 used here, as compared to N_max = 6 used previously, an obligatory small value required for the realistic application of Eigen's formalism with the available computing power; (b) a large difference in repertoire size (N_G = 100 here versus N_G = 10 in the earlier study); (c) the performance here of 10,000 random simulations, considered essential for proper statistical rigor, as compared to only a single simulation done previously. Both points (b) and (c) provide a significant edge to the present simulations in sampling the β interaction space, which allows drawing conclusions with higher certainty. In the future it will be interesting to consider additional methodologies to exert external selection. One could be a variant of the presently used method, whereby the β network will be biased by a constant factor and not employing target similarity-oriented bias. Another could be biasing the environmental concentration ρ_i (Equation 2) by a constant factor based on the molecules that are contained in the target compotype.

The GARD model embodies the inheritance of compositional information in the realm of a lipid world scenario for early evolution [20, 21, 23, 27, 48, 49, 51, 55–57]. The GARD has recently been pursued in several additional publications [20, 39, 68, 70] and has been chosen as an archetypal metabolism-first realization [68]. This suggests that despite being a simulated toy model, the GARD has sufficient complexity to shed light on some important questions in the field of prebiotic origins. In the present work an attempt is made to shed further light on some of the GARD's evolutionary features. It is expected that the present insights will become instrumental in further efforts to extend the GARD beyond the monomer world [54], as has been preliminarily explored [55]. This might be necessary to reveal the capacity of the GARD model to capture the much-needed open-ended attributes of natural selection and evolution.

We thank Raphael Zidovetzki and Tsviya Olender for enlightening comments, Aron Inger for assistance in coding, and Ágnes Tóth Petróczy, Azhar Ali Shah, Hugues Bersini, Leong Ting Lui, and Natalio Krasnogor for discussions. This work is partly supported by EU-FP7 project MATCHIT and by the Crown Human Genetics Center at the Weizmann Institute of Science.

The fitting procedure is as follows: (1) Calculate H₀ as the mean of c(Δt) in the interval [GEN/4, GEN/2]. (2) Guess τ^‡ as the first instance c(Δt) drops below H₀. (3) Smooth the c(Δt) tail by forcing: c(Δt > τ^‡) = H₀. (4) Fit an exponential (Equation 7) to the smoothed c(Δt), using nonlinear least squares with a tolerance of 10^-5.

Examples are given in Figure 12 in  Appendix A.4.

See Table 2.

See Figures 10,¹¹¹²¹³–14.

See Figures 15–16.

See Figures 17 and Table 3.