A formal model of autocatalytic sets emerging in an RNA replicator system

Background: The idea that autocatalytic sets played an important role in the origin of life is not new. However, the likelihood of autocatalytic sets emerging spontaneously has long been debated. Recently, progress has been made along two different lines. Experimental results have shown that autocatalytic sets can indeed emerge in real chemical systems, and theoretical work has shown that the existence of such self-sustaining sets is highly likely in formal models of chemical systems. Here, we take a first step towards merging these two lines of work by constructing and investigating a formal model of a real chemical system of RNA replicators exhibiting autocatalytic sets. Results: We show that the formal model accurately reproduces recent experimental results on an RNA replicator system, in particular how the system goes through a sequence of larger and larger autocatalytic sets, and how a cooperative (autocatalytic) system can outcompete an equivalent selfish system. Moreover, the model provides additional insights that could not be obtained from experiments alone, and it suggests several experimentally testable hypotheses. Conclusions: Given these additional insights and predictions, the modeling framework provides a better and more detailed understanding of the nature of chemical systems in general and the emergence of autocatalytic sets in particular. This provides an important first step in combining experimental and theoretical work on autocatalytic sets in the context of the orgin of life.


Background
Recently, significant new experimental results on spontaneous network formation among cooperative RNA replicators were reported [1]. These results continue and strengthen a line of ongoing work on creating autocatalytic sets in real chemical systems [2][3][4][5]. Moreover, they show the plausibility and viability of the idea of autocatalytic sets, especially in the context of the origin of life, as already developed in various forms several decades ago [6][7][8][9][10][11][12].
However, such chemical experiments, important as they are, are difficult, costly, and time-consuming to perform. In contrast, in our own work we have developed a theoretical framework of autocatalytic sets, which we consider a necessary condition for the origin of life, that can be studied computationally and mathematically [13][14][15][16][17][18][19][20]. This framework has provided many important insights into the emergence and structure of autocatalytic sets in its own right, and is considered to provide theoretical support for experimental observations [1,21].
The formal framework has, so far, mostly been applied to an abstract model of a chemical reaction system known as the binary polymer model. Even though this model already has a fair amount of chemical realism (for example, some recent experiments are an almost literal chemical implementation of the binary polymer model [5]), a direct application to a real chemical system was still lacking. Here, we construct and analyze a formal model version of the recent experimental RNA replicator system [1], and show how this results in: • an accurate reproduction of experimental results, • the correction of a misinterpretation of some of the original results, • additional results and insights that could not be obtained from the experiments alone, and • testable predictions about the behavior of the chemical system.
With this, we claim that the formal framework can be applied directly and meaningfully to real chemical systems, generating additional insights and predictions that are hard to obtain from experiments alone, thus providing a better understanding of real (bio)chemical systems. This represents an important first step towards merging experimental and theoretical work on autocatalytic sets.

Chemical reaction systems and autocatalytic sets
We begin by briefly reviewing the relevant definitions and main results of the formal framework. A chemical reaction system (CRS) is defined as a tuple Q = {X, R, C} consisting of a set of molecule types X, a set of chemical reactions R, and a catalysis set C indicating which molecule types catalyze which reactions. We also consider the notion of a food set F ⊂ X, which is a subset of molecule types that are assumed to be freely available from the environment (i.e., they do not necessarily have to be produced by any of the reactions). Informally, an autocatalytic set (or RAF set) is now defined as a subset R ⊆ R of reactions (and associated molecule types) which is: 1. reflexively autocatalytic (RA): each reaction r ∈ R is catalyzed by at least one molecule type involved in R , and 2. food-generated (F): all reactants in R can be created from the food set F by using a series of reactions only from R itself.
A more formal (mathematical) definition of RAF sets is provided in [14,17], including an efficient algorithm for finding RAF sets in a general CRS. It was shown (using the binary polymer model) that RAF sets are highly likely to exist, even for very moderate levels of catalysis (between one and two reactions catalyzed per molecule, on average) [14][15][16], and that this result still holds when a more realistic "template-based" form of catalysis is used [17,18].
The RAF sets that are found by the RAF algorithm are called maximal RAF sets (maxRAFs). However, it was shown that a maxRAF can often be decomposed into several smaller subsets which themselves are RAF sets (subRAFs) [19]. If such a subRAF cannot be reduced any further without losing the RAF property, it is referred to it as an irreducible RAF (irrRAF). The existence of multiple autocatalytic subsets can actually give rise to an evolutionary process [22], and the emergence of larger and larger autocatalytic sets over time [19].

A model of the RNA replicator system
The main idea behind the RNA replicator system reported recently [1] is the assembly (ligation) of ribozymes (catalytic RNA molecules) from two smaller RNA fragments. These ribozymes can then catalyze the assembly of other ribozymes (or in some cases their own assembly). Which ribozymes catalyze which assembly reactions is determined by one specific nucleotide in the "guide sequence" of the (potential) catalyst and one other specific nucleotide in the "target sequence" of a reactant. If these two nucleotides are each other's base pair complement, then the given ribozyme can catalyze the given reaction. Starting with RNA fragments with different nucleotides in these guide and target sequences, a mixture of auto-and cross-catalytic RNA replicators evolves over time (by means of the assembly reactions taking place), in which cooperative networks (autocatalytic sets) form spontaneously [1].
The ribozymes in this system are labeled MjN, where M denotes the specific nucleotide (A, C, G, or U) in the guide sequence of an RNA molecule, and N denotes the specific nucleotide in its target sequence [1].
The value of j denotes the specific location where the ribozyme was assembled from two smaller RNA fragments (there are three possible locations in the original experiments). However, as in some of the experimental results [1], for the purposes of looking for autocatalytic sets, this value can be ignored. So, in total there are 4 × 4 = 16 possible ribozymes MjN. For example, GjA and AjC are two such ribozymes.
The reaction set X in the formal model of this RNA replicator system contains these 16 MjN ribozymes plus the smaller RNA fragments from which they are assembled. In fact, these RNA fragments are the subset of molecules that forms the food set F . In the model, we simply lump these fragments together into just two food molecule types (i.e., each ribozyme is the assembly of two "generic" fragments). For most of the results shown here this suffices, although a refinement will be made later on in one of the computational experiments (below).
The reaction set R in the model consists of the 16 assembly reactions that create the ribozymes from food molecules. For simplicity, the notation MjN will denote both a molecule type (ribozyme) as well as the reaction that created it. Finally, the catalysis set C consists of all molecule/reaction pairs where the nucleotide M in the guide sequence of a molecule is the base-pair complement of the nucleotide N in the target sequence of a reaction. For example, molecule GjA can catalyze the reaction that produces molecule AjC, since G and C are complementary. The catalysis set C thus consists of all such matching pairs The full mathematical definition of the CRS Q = {X, R, C} of the RNA replicator system is thus as follows: Note that since the reactants in all 16 reactions in R are food molecules, every subset R ⊆ R is automatically food-generated (F). Therefore, identifying possible RAF sets in this system only requires checking the (RA) part of the definition. We now analyze this model in detail and compare the results with those from the original experiments [1].

Results and discussion
The existence of RAF sets Applying the RAF algorithm to the reaction set R (as defined above), returns the set R itself. In other words, the system as a whole (all 16 reactions) forms a maximal RAF set. Moreover, this would be expected, given the extent of catalysis, even if the actual assignment of catalysis was randomized. In fact, under such a (random) model the probability that R forms an RAF is approximately ( This can be derived as follows. For molecule x and reaction r, let I(x, r) be the indicator function defined by: Now suppose that Q = (X, R, C) is a catalytic reaction system. Then under the random model [13][14][15], the collection of indicator functions (I(x, r) : x ∈ X, r ∈ R) are independent and identically distributed Bernoulli random variables. Thus, if we let p = P[I(x, r) = 1], then the probability q r that a given reaction r ∈ R is catalyzed by at least one molecule from X is: Moreover, if we set the expected amount of catalysis in the random model to the value observed in Q then: Thus, Eqn. (1) becomes: Now, if R is F-generated, and if X is the set of molecules involved (as products or reactants) in R, then R is an RAF precisely if every reaction in R is catalyzed by at least one molecule from X. This latter probability, under the random model, is: which, for |C| = 64 and |R| = 16 (as in the RNA replicator model) evaluates to (1 − e −4 ) 16 , as claimed.
Notice that Eqn. (4) is essentially independent of |X| (the number of molecules) provided this is not too small.
Even though the entire reaction set R is a maximal RAF set in itself, RAF sets can often be decomposed into smaller RAF subsets [19]. Two such subsets were presented in the original experimental results ( [1] Fig. 4). However, a comparison of those diagrams with the formal model indicates that they are incomplete, as several catalysis arrows are missing. Fig. 1 (here) shows the corrected diagrams, with all catalysis arrows included. The first corrected diagram (Fig. 1, left) contains an RAF set of size seven (enclosed within the box): each node (reaction) is catalyzed by at least one of the members of the set, which satisfies the (RA) part (and, as noted above, the (F) part is automatically satisfied). Thus, the conclusion from the original experiments that "at the 1h time point, no closed network was possible" [1] is a misinterpretation due to the omission of several catalysis arrows. As the corrected diagram shows, at the 1h time point an RAF set of size seven already exists (Fig. 1, left), which has grown to size 11 at the 8h time point (Fig. 1, right, which forms an 11-reaction RAF set).

The structure of RAF sets
As mentioned, RAF sets can often be decomposed into smaller subRAFs. Indeed, even the 7-reaction RAF subset in Fig. 1 (left) itself is composed of many smaller subRAFs. Previously we introduced a formal method to identify and classify RAF subsets [19], resulting in a so-called Hasse diagram that visualizes the partially ordered set of all possible subRAFs and their subset relations. Applying this method to the 7-reaction RAF set in Fig. 1 (left) results in the Hasse diagram shown in Fig. 2.
This Hasse diagram contains 68 nodes, i.e., there are 68 possible subRAFs in the given 7-reaction RAF set.
Note that there are 2 7 = 128 possible subsets of a set of seven elements. So, more than half of these actually form RAF sets themselves, which shows how "rich" and diverse the RNA replicator system really is. In fact, if the same ratio (about half) holds for the full 16-reaction maximal RAF set, then we can expect more than 2 16 /2 ≈ 32, 000 nodes (subRAFs) in the Hasse diagram of the maxRAF, i.e., far too many to visualize in a meaningful way.
The edges in the Hasse diagram of an RAF set show the many possible ways in which the full RAF set can be built up from smaller subsets, or, in other words, how RAF sets can emerge and evolve [19,22]. Which of these trajectories is actually followed depends on, for example, initial conditions and stochastic events such as "spontaneous" (uncatalyzed) reactions. As shown above, in the original experiment the system went through a stage of a particular 7-reaction RAF set, which then grew to an 11-reaction RAF set.
However, what the Hasse diagram suggests, given the many possible subRAFs and ways of combining them into larger RAFs, is that when the experiment is repeated, most likely a different trajectory will be followed. This is a testable prediction that follows directly from the formal model and its results.

The emergence of RAF sets
Performing the actual chemical experiments in a laboratory is costly and time-consuming. However, we can use the formal model to simulate molecular flow on the reaction network [20]. Using the well-known Gillespie algorithm [23,24], we performed such simulations on the full 16-reaction model, starting with an initial supply of food molecules (RNA fragments) only. Fig. 3 shows the result of one such simulation (time and concentration are in arbitrary units). For simplicity, we set all reaction rates equal, but with a factor c difference between catalyzed and uncatalyzed ("spontaneous") reactions. We tried various values for this factor c (anywhere from c = 2 to c = 100), but qualitatively there is little difference in the overall dynamics, except that everything happens at longer or shorter time-scales, depending on the exact value of c.  Fig. 4). Indeed, the system has gone through a different trajectory in the simulation run compared to the chemical experiment.
With these results we do not intend to claim that the simulation is an exact reproduction of the original experiment, where at regular intervals 10% of the current solution was transferred to a new solution of fresh RNA fragments [1]. Although we can include such "transfer" steps in our simulations as well, we mainly wish to show that the overall process of going through larger and larger subRAFs is accurately reproduced by the model, and to confirm that in each repetition of the experiment (simulation), a different trajectory is indeed followed, as suggested by the Hasse diagram.

The advantage of RAF sets
As a further demonstration of the value of the modeling approach, we consider the "cooperation versus selfishness" question [1]. Using a variant of the formal model, we investigate the system of three cooperative molecules and three equivalent selfish ones, competing for the same food resources. The reaction graph of this system is shown in Fig. 5 (inset). Ribozymes E1 and S1 are assemblies of one pair of food molecules (f 1 and f 2 ), molecules E2 and S2 are assemblies of another pair of molecules (f 3 and f 4 ), and molecules E3 and S3 of yet another pair (f 5 and f 6 ). Molecule E1 catalyzes the assembly of E2, E2 that of E3, and E3 that of E1, in a cooperative way (forming an RAF set of size three). Molecules S1, S2, and S3 each catalyze their own assembly. So, the cooperative (green) part of the system (RAF set) competes with the selfish (red) part of the system for the same food molecules, as in the original experiment [1].  Fig. 2a, mixed system). In Fig. 5, the difference between cooperation and selfishness is even larger, although the difference is not always this large in each simulation run. Fig. 6 (top) shows the results of another simulation on the same reaction network where the difference is much smaller. In fact, occasionally the selfish (red) system actually outperforms the cooperative (green) one, an example of which is shown in Fig. 6 (bottom).
There is a simple explanation for the difference between cooperation and selfishness in this model. To get the cooperative system (RAF set) going, only one of the three (green) reactions has to happen "spontaneously", i.e., uncatalyzed; this is always possible, but at a lower rate (by a factor c) than a ribozyme-catalyzed reaction. This happens around time point 0.14 in the simulation shown in Fig. 5.
However, in the selfish system, all three (red) reactions have to happen spontaneously at least once to get the full system going. This results in an almost three times longer waiting time on average. In the simulation shown in Fig. 5, the third spontaneous (uncatalyzed) red reaction happens around time point 0.35. By that time, however, the (green) RAF set has already built up enough "momentum" to outcompete the selfish (red) system, i.e., enough molecules of types E1, E2, and E3 are already around to increasingly catalyze each other's assembly. However, due to the stochastic nature of the system (having to wait for uncatalyzed reactions to happen), occasionally the selfish (red) system gets a head-start (with low probability) and outcompetes the cooperative (green) system.
Finally, there is also an advantage for the RAF set in terms of robustness. Suppose that at some point all molecules of type E3 and S3 are removed from the system, and a new supply of their food molecules (f 5 and f 6 ) is provided. This happens around time point 2 in Fig. 5, at the sudden dip in the concentrations.
Since there are still E1 and E2 molecules present, the (green) RAF set quickly recovers without any delay.
However, the (red) selfish system again has to wait until reaction S3 happens uncatalyzed at least once (which still has not happened by time point 3). Clearly, this suggests that RAF sets are more robust than selfish systems against perturbations, another prediction that can be tested with real chemical experiments.

Conclusions
We have taken the chemical RNA replicator system described recently [1] and formalized and investigated it within our mathematical RAF framework. As the experimental results showed, autocatalytic sets seem to form spontaneously in this chemical system [1]. The existence of RAF sets in this system is verified by the formal model. In fact, many of the experimental results are accurately reproduced by the model, such as the emergence of larger and larger RAF sets over time, and the advantage of cooperative systems over selfish systems when they compete for the same resources. Of course there are many refinements and additional details that can be added to the current model, but even in its most basic form it already captures the main structural and dynamical properties of the real chemical system.
Moreover, the modeling approach provides several additional results and insights. First, it allows for a correction in the misinterpretation of some of the original experimental results. Second, it shows the "richness" of the RNA system in terms of the many subRAFs and ways they can be combined into larger subRAFs (as visualized by the Hasse diagram). Third, this "richness" suggests the testable prediction that each repetition of the experiment will, most likely, follow a different trajectory towards a realization of the maximal RAF set, which is confirmed computationally by the molecular flow simulations presented here.
Lastly, it provides a simple explanation for the advantage of RAF sets over selfish systems, together with another testable prediction that RAF sets are more robust than selfish systems against perturbations.
In conclusion, we have shown that the formal RAF framework can be directly and meaningfully applied to real chemical systems, generating additional insights that are hard or even impossible to obtain from experiments alone, and thus provides a more detailed understanding of the nature of chemical systems in general and the emergence of autocatalytic sets in particular. This forms an important and much needed first step towards merging experimental and theoretical lines of work on autocatalytic sets in the context of the origin of life.

Competing interests
The authors declare that they have no competing interests. The bottom row represents subRAFs of size one (the three autocatalytic reactions in Fig.1 (left)). Each next row up represents subRAFs of a larger size, up to the full 7-reaction RAF set at the top. Due to space constraints, only the four irreducible RAFs in this diagram (the four nodes that do not have an incoming edge from a lower level) and the full 7-reaction node at the top are labeled.  Top: an example where the difference between the two systems is only minimal. Bottom: an example where the selfish system actually outcompetes the cooperative one.