Alex Blokhuis

David Lacoste

Philippe Nghe

Universal motifs and the diversity of autocatalytic systems

Autocatalysis is essential for the origin of life and chemical evolution. However, the lack of a unified framework so far prevents a systematic study of autocatalysis. Here, we derive general stoichiometric conditions for catalysis and autocatalysis in chemical reaction networks from basic principles. This allows for a classification of minimal autocatalytic motifs called cores. While all known autocatalytic systems indeed contain minimal motifs, the classification also reveals hitherto unidentified motifs. We further examine conditions for kinetic viability of such networks, which depends on the autocatalytic motifs they contain and is notably increased by internal catalytic cycles. Finally, we show how this framework extends the range of conceivable autocatalytic systems, by applying our stoichiometric and kinetic analysis to autocatalysis emerging from coupled compartments. The unified approach to autocatalysis presented in this work lays a foundation towards the building of a systems-level theory of chemical evolution.

EXAMPLES, DEFINITIONS AND CONVENTIONS

Catalysis and autocatalysis

The following reactions have the same net mass balance but a different status regarding catalysis: A (I) ----B, A + E (II) ----B + E, A + B (III) ------2 B.

(S1) Since no species is both a reactant and product in reac-Supp

The capacity of living systems to replicate themselves is rooted in a chemistry that makes more of itself, i.e. an autocatalytic system. Autocatalysis appears to be ubiquitous in living systems from molecules to ecosystems [1]. It is also likely to have been continually present since the beginning of life and is invoked as a key element in prebiotic scenarios [2][3][4][5]. Surprisingly, autocatalysis is considered to be a rarity in chemistry [6]. Developments in systems chemistry are changing this view, with an increasing number of autocatalytic systems synthesized de novo [7][8][9]. Chemical replicators have been endowed with biomimetic properties such as protein-like folding [10] and parasitism [11]. Autocatalysis has also found technological applications, e.g. enantiomer enrichment and acid amplification [12][13][14].

Understanding autocatalysis represents a primary challenge for theory. Models based on autocatalysis were first built to explain a diversity of dynamical behaviors in so called dissipative structures, such as bistable reactions [15], oscillating reactions, and chemical waves [16]. Autocatalysis then became a central topic in the study of self-replication dynamics in biological and prebiotic systems [3,[17][18][19] (see [20][21][22] for recent reviews).

Despite this history, a unified theory of autocatalysis is still lacking. Such a theory is needed to understand the origins, diversity and plausibility of autocatalysis. It would also provide design principles for artificial autocatalytic systems. Here, we present a framework that unifies the different descriptions of autocatalysis and is based on reaction network stoichiometry [23][24][25][26][27].

Let us start from basic definitions in chemistry as established by IUPAC [28] (see SI Sec. I for full definitions), where autocatalysis is a particular form of catalysis: A substance that increases the rate of a reaction without modifying the overall standard Gibbs energy change (∆G • ) in the reaction; the process is called catalysis. The catalyst is both a reactant and product of the reaction. Catalysis brought about by one of the products of a (net) reaction is called autocatalysis.

From this definition, we derive conditions to determine whether a subnetwork embedded in a larger chemical network, can be catalytic or autocatalytic. These conditions provide a mathematical basis to identify minimal motifs, called autocatalytic cores. We found that cores have five fundamental categories of motifs. They allow classification of all previously described forms of autocatalysis, and also reveal hitherto unidentified autocatalytic schemes. We then study the kinetic conditions, which we call viability conditions, under which autocatalytic networks can appear and be maintained on long times. We find that networks have different viabilities depending on their core structure, and notably that viability is increased by internal catalytic cycles. Finally, we expand the repertoire of autocatalytic systems, by demonstrating a general mechanism for its emergence in multicompartment systems (e.g. porous media, vesicles, multiphasic systems). This mechanism strongly relaxes chemical requirements for autocatalysis, making the phenomenon much more diverse than previously thought. tion (I), it should be regarded as uncatalyzed. Reactions (II) and (III) instead contain species which are both a reactant and a product, species E in reaction (II) and species B in reaction (III) and following the definition above, these species can be considered as catalysts. In reaction (II), the amount of species E remains unchanged, in contrast to the case of reaction (III), where the species B experiences a net production. For this reason, reaction (III) represents genuine autocatalysis. Although reaction (II) is usually referred to as simply catalyzed in the chemistry literature, we propose to call it an example of allocatalysis to contrast it with the case of autocatalysis, catalysis being common to both.

We emphasize that stoichiometric considerations are necessary but not sufficient to characterize catalysis, which according to the definition should also accelerate the rate of the net reaction. In the following, we will first generalize the stoichiometric conditions, then examine kinetic ones.

Stoichiometric matrix and reaction vectors

Reaction networks are represented as a stoichiometric matrix ν ν ν [23,26], in which columns correspond to reactions and rows to species. The entries in a column are the stoichiometric coefficients of the species participating in that reaction, the coefficient is negative for every species consumed and positive for every species produced. A reaction vector g g g = [g 1 , .., g r ] T results in a change of species numbers ∆n n n = ν ν ν • g g g. The support of g g g, denoted supp(g g g), is the set of its non-zero coordinates. A reaction cycle is a non-zero reaction vector c c c such that no net species number change occurs : ν ν ν • c c c = 0 0 0, or equivalently, c c c belongs to the right null space of ν ν ν. Vectors b b b T belonging to the left null space of ν ν ν induce conservation laws, because in that case b • n b • n b • n represents a conserved quantity. The case of all coefficients b k nonnegative is referred to as a masslike conservation law. For example in Fig. 1a, conserved quantities are n E + n EA (catalysts) and n A + n EA + n B (total compounds).

Lastly, catalyzed reactions may not always be distinguished from uncatalyzed one in the stoichiometric matrix. For instance, in reactions (II-III), catalysts cancel on each side leading to the same column vector as for (I). This is avoided by describing catalysis through a sequence of reactions steps from which it emerges, so that a participating species is either a reactant or a product:

A + E IIa ----EA IIb ----E + B, (S2) A + B IIIa ------AB IIIb ------2 B . ( S3
)
We call this convention non-ambiguity and assume henceforth that it is respected.

EA A EA A B E EA -1 0 0 -1 1 1 1 -1 AB B A B AB -1 -1 1 -1 2 0 E E a) b) c) f ) e) d) B + EA E + AB A B + AB 2B 1'' 2'' 1' 2' 1' 2' 1' 2' 1'' 2'' 1'' 2''

CATALYSIS AND AUTOCATALYSIS IN STOICHIOMETRIC MATRICES

In this section, we will consider any possible submatrix ν ν ν of ν ν ν, the stoichiometric matrix of a reaction network, and ask whether the stoichiometry of the corresponding subnetwork, called a motif, is compatible with the definitions of allocatalysis or autocatalysis. Note that such identification neither makes a priori assumptions on the values and signs of reaction vector coefficients, nor on kinetics, or on which species are catalytic or not. A matrix ν ν ν is a restriction of ν ν ν to certain rows and columns, which respectively correspond to the species and reactions of the motif under consideration.

The restriction of the rows means that the species of ν ν ν are separated into internal species of the motif (rows of ν ν ν) and external species (remaining rows of ν ν ν). These external species could be, in some cases, chemostatted [26], and represent feedstock compounds, also called the food set [29], and waste from the point of view of internal species of the motif. In Fig. 1, external species have been colored in blue, while stoichiometric submatrices have been boxed in yellow. Restriction of columns separates reactions which are part of the motif and those which occur outside of it. A motif such that each of its reactions has at least one reactant and at least one product is called autonomous. This means that every column of ν ν ν contains a positive and a negative coefficient. Below, we pose autonomy as a condition for catalysis. Indeed, it ensures that the production of any species of the motif is conditional on the presence of other chemical species of the motif. Otherwise, rate acceleration would be allowed unconditional on an already present substance, in opposition to the definition of catalysis. Autonomy is less restrictive than former conditions for autocatalysis [24], and is similar to the siphon concept in Petri Nets [30], but without assumption on reaction signs (see SI Appendix, Sec. 1). Note that it does not forbid that reactions outside of the motif produce species of the motif.

Criterion for allocatalysis

By definition, allocatalysis is an ensemble of reactions by which a set of species remain conserved in number (the catalysts) while other external species undergo a turnover which changes their numbers. This leads to the following conditions:

There exists a set of species S S S, a submatrix ν ν ν of ν ν ν restricted to S S S, and a non-zero reaction vector c c c such that

: i) ν ν ν is autonomous; ii) supp(c c c) is included in the columns of ν ν ν; iii) c c c is a reaction cycle of ν ν ν (ν ν ν • c c c = 0 0 0)
, and; iv) ν ν ν •c c c = 0. The members of S S S which participate in c c c (i.e. that are consumed and produced) are called allocatalysts, c c c an allocatalytic cycle and ν ν ν an allocatalytic matrix. Condition (i) has been discussed above. Condition (ii) expresses the involvement of the catalysts in the reactions c c c, where all columns of ν ν ν are non-zero due to (i), so that all reactions of c c c involve catalysts. Condition (iii) expresses the conservation of catalysts and (iv) the net reaction. Since the reaction cycle c c c is a cycle of the reduced matrix but not of the original matrix, some authors have qualified it as emergent and shown that it can establish a nonequilibrium steady state driven by the turnover of the external species [26]. Note that being allocatalytic is not a property of the sub-matrix ν ν ν alone but involves the larger matrix ν ν ν as imposed by condition (iv).

Criterion for autocatalysis

By definition, autocatalysis is the process by which a combination of reactions involves a set of species which all increase in number conditional on species in the set itself (the autocatalysts), while other species undergo a turnover. This leads to the following conditions:

There exists a set of species S S S, a submatrix ν ν ν of ν ν ν restricted to S S S, and a reaction vector g g g such that: i) ν ν ν is autonomous, ii) all coordinates of ∆n n n = ν ν ν • g g g are strictly positive, or equivalently, ν ν ν has no mass-like conservation laws. The members of S S S consumed (and produced) by g g g are called autocatalysts, g g g an autocatalytic mode and ν ν ν an autocatalytic matrix.

Condition (i) ensures the conditionality of the reactions on autocatalysts, as it forbids cases where species of S S S are produced from external reactants only, thus playing the role of conditions (i) and (ii) in the definition of allocatalysis. Condition (ii) expresses the increase in autocatalyst number. The equivalence between the two formulations of condition (ii) is an immediate consequence of Gordan's theorem [31]. Importantly, the second formulation of (ii) does not involve an autocatalytic mode g g g, so that (i) and (ii) can be expressed as properties of a matrix itself, in contrast with allocatalysis. This allows us to look for minimal autocatalytic motifs, which we do next. Note that external species must feed the autocatalytic system in order to guarantee the net mass increase imposed by condition (ii).

Autocatalytic cores

An autocatalytic core is an autocatalytic motif which is minimal because it does not contain any smaller autocatalytic motif. Consequently, an autocatalytic system is either a core, or it contains one or several cores. The stoichiometric conditions show that characterizing cores is equivalent to finding all autonomous matrices whose image contains vectors with only strictly positive components. This well-posed formulation allowed us to show that the stoichiometric matrix ν ν ν of an autocatalytic core must verify a number of non-obvious properties reported below and demonstrated in the SI Appendix Sections 2 and 3.

First, ν ν ν must be square (the number of species equals the number of reactions) and invertible. The inverse has a chemical interpretation. By definition of the inverse, the k-th column of ν ν ν -1 is a reaction vector such that species k increases by one unit, making it an elementary mode of production. Likewise, the reaction vector obtained by summing the columns of ν ν ν -1 leads to a net increase by one unit of every autocatalyst, which thus represents an elementary mode of autocatalysis. This shows how stoichiometry informs on fundamental modes of autocatalysis [27].

Second, every forward reaction of a core involves only one core species as a reactant. While this excludes reactions between two different core species, a single core species may react with itself. As ν ν ν is square, this also implies that every species of a core is consumed (none is only produced), thus is an autocatalyst. Furthermore, every species is the reactant of a single reaction. Overall, every species is uniquely associated with a reaction as being its reactant, so that ν ν ν admits a representation with a negative diagonal and zero or positive coefficients elsewhere, at least one coefficient of each column being strictly positive to ensure autonomy.

These properties are constraining enough to allow an exhaustive enumeration of reaction graphs that are cores. Autocatalytic cores are found to belong to five categories, denoted as Type I to Type V. Fig. 2a represents typical members of each category as reaction hypergraphs (see SI Appendix Fig. S1 for general cases). As can be seen in these graphs, all minimal motifs contain a fork, which is a reaction with a single reactant and two products. The presence of at least one fork is consistent with the intuition that autocatalysis requires reaction steps that amplify the amount of autocatalysts. In Type I cores, the fork ends with two copies of the same compound, whereas in Types II to V, forks end with different compounds. Note that Type I cores can be seen as a particular case of Type II if one formally allows different nodes to represent a same chemical species. It is obligatory for Types II and III to contain one fork with two distinct products, for Type IV to contain two such forks, and for Type V to contain three. In Fig. 2a, the orange squares on the links between the nodes indicate that these links could contain further nodes and reactions in series, provided certain rules on cycles given in the next paragraph.

The five types differ in their number of graph cycles[32] and the way these cycles overlap. Type I consists of a single graph cycle that is weight-asymmetric, defined as the product of the stoichiometric coefficients of its reaction products being different than that of its reactants. Types II-V can be described as two, respectively three, overlapping graph cycles, where any such graph cycle involving a strict subset of the core species must be an allocatalytic cycle, i.e. weight-symmetric (it would otherwise be of Type I, contradicting minimality).

Unification of autocatalytic schemes

The stoichiometric characterization of autocatalysis provides a unified approach to autocatalytic networks reported in the literature. The examples below are further detailed in SI Appendix Sec. 4. The formose reaction is a classic example of autocatalysis known to contain many autocatalytic cycles [34]. Fig. 2b andc show Type I and III cores both found in the formose reaction. Similarly, autocatalytic cores of Type I and III can be found in the Calvin cycle and reverse Krebs cycle (SI Appendix Fig. S4). Some reaction steps Fig. 2b may be catalyzed externally (e.g. by enzymes, base, ions), but external catalysis in general does not alter the core. By the same token, proposed examples of auto-induction introduced in [21,35] contain Type I and III cores (SI Appendix Fig. S3).

In the GARD (Graded Autocatalysis Replication Domain) model for self-enhancing growth of amphiphile assemblies [4,5], all underlying autocatalysis is described (SI Appendix Fig. S6) by Type I cycles with one fork and Type II cycles built up from sequential nonoverlapping allocatalytic cycles (cross-incorporation, such as N 3 in Fig. 3). More generally, when such catalytic cycles are compactly written as single reactions as in (S1), they can be treated in the RAF (Reflexively Autocatalytic and Food-generated) framework [29], where they form irreducible RAF-sets [36]. This formally establishes the recently suggested link [5,37] between these models.

Another reported form of autocatalysis is 'chemical amplification' due to cavitands [38]. The mechanism involves a reactive compound in a molecular cage, whose free counterpart can react to form two species that exchange with the caged species, thus amplifying its release. We find that this process can be described within our framework and corresponds to a Type III core (SI Appendix Fig. S5). Overall, previously described autocatalytic schemes comprise Types I, II and III. We have not yet found examples of Types IV and V.

VIABILITY OF AUTOCATALYTIC NETWORKS

Stoichiometric conditions do not guarantee that autocatalysts within motifs amplify. Whether an initial autocatalyst amplifies or degrades depends on kinetic considerations. To address this so-called fixation problem [17,22], we examined the probability P ex of extinction (or 1 -P ex of fixation) of species within autocatalytic motifs, as a function of transition probabilities of reaction steps. Considering a homogeneous system with a steady supply of reactants, several authors have noted that in the highly dilute autocatalyst regime, appreciable rates require first-order autocatalysis [17,22,39], i.e. each for-ward reaction step only involves one autocatalyst. Among first-order order networks, fixation models have so far focused on Type I networks (e.g. Fig. 2b), which have a single graph cycle containing n species. In a transition step, a given species may either proceed irreversibly to the next species or disappear as a result of degradation. King found that if every reaction step k among n steps of the cycle has a success probability Π + k (1 -Π + k being the degradation probability), fixation is possible for a doubling probability p 2 = n k=1 Π + k ≥ 1/2 [39]. This minimum value of p 2 above which fixation is possible is called the decay threshold [19,40]. Bagley et al. [17] used birth-death processes to derive P ex for an autocatalytic loop containing one species (n = 1). Schuster reported detailed time-dependent statistics for such networks in various contexts [22].

Here, we extend the treatment of the fixation problem so as to include reversible reactions and networks beyond Type I using the theory of branching processes [41]. In these stochastic processes, an autocatalytic species X s is, after a sequence of reaction steps in the network, replaced by k copies. Reaction sequences yielding k copies happen with a probability p k , such that

X s p0 -→ ∅, X s p1 -→ X s , ... X s p k -→ kX s , ... . (S4)
The probability P ex that X s goes extinct is then the probability that its k descendants independently go extinct:

P ex = p 0 + p 1 P ex + p 2 P 2 ex + ... = ∞ k=0 p k P k ex . (S5)
The main difficulty here is to derive p k from transition probabilities Π k . A procedure for this is given in SI Sec. 5, where branching processes are constructed from reaction networks. Below, we exemplify this method by generalizing known results for Type I networks, solutions for other networks being detailed in the SI Appendix Sec. 5. We then apply it to compare the P ex of autocatalytic motifs which differ in their core structures.

Reversible Type I cycles

Consider a Type I cycle consisting of n reaction steps, such as N 1 in Fig. 3b, and let us start at the first step with species X 1 (marked node). Ultimately, X 1 will either be successfully converted and yield 2X 1 or be degraded prematurely, which simplifies (S4) to

∅ p0 ← -X 1 p2 -→ 2X 1 , ( S6
)
with p 0 + p 2 = 1. The overall outcome described by (S6) corresponds to the simplest type of branching process: a birth-death process. (S5) then becomes a quadratic equation that yields

P ex = 1 p2 -1, p 2 ≥ 1 2 , 1, p 2 < 1 2 .
(S7)

This generalizes Bagley et al's observation for Type I networks to n > 1 and reversible reactions. For reversible reactions, p 2 is found by considering all possible sequences of forward and backward reactions along the cycle. From X k , let Π - k be the transition probability to revert to X k-1 , and Π + k to convert to X k+1 . We have

p 2 = n k=1 Π + k Γ k , (S8) Γ k+1 = ∞ s=0 (Π - k+1 Γ k Π + k ) s = 1 1 -Π - k+1 Γ k Π + k , (S9)
where Γ k recursively (Γ 1 = 1) counts the statistical weight of all back-and-forth trajectories from X k to itself, in terms of Π - k and Π + k . In the irreversible reaction limit Π - k → 0, Γ k → 1 King's expression for p 2 [39] is recovered.

b) a)

P ex P ex N 1 N 2 N 3 N 4 N 5

Viability of autocatalytic cores

To investigate how autocatalytic motif structure affects survival, we calculated P ex for five different cores[42] (N 1 to N 5 , Fig. 3): they are of equal size (6 reaction steps, 6 species), all reactions proceed irreversibly with the same success probability ζ, which plays a similar role as the transition probability Π + k in the example above and is sometimes called specificity [19,39,40]. Fig. 3 highlights how P ex depends on ζ for each core structure. The highest ζ for extinction (P ex = 1) is observed for the Type I cycle N 1 , and progressively lower values are found for N 2 to N 4 , which are all of Type II. Type V network N 5 tolerates the lowest specificity ζ before extinction, sustaining almost three times higher failure rates 1-ζ than N 1 . These differences can be qualitatively understood by counting the minimum number of steps needed to produce more autocatalysts. In respective order, networks N 1 to N 5 in Fig. 3b do so in six, four, three, three and two steps. In particular, given their symmetries, the P ex of N 3 and N 5 have the same dependence on ζ as a 3 and 2-membered Type I cycle, respectively.

It has been suggested that large networks are disfavored in general [39]. The examples of Fig. 3 indicate that this can be counterbalanced by the presence of more allocatalytic cycles in the network. This is in particular the case for autocatalytic sets, as every net reaction must participate in an allocatalytic cycle.

EXTENSIONS: MULTICOMPARTMENT AUTOCATALYSIS

We finally show how stoichiometric criteria allow the identification of autocatalysis that emerges from compartments coupled via selective exchange, as found in systems comprising vesicles, pores, emulsions and complex coacervates [43]. The reaction network in Fig. 4a is incapable of autocatalysis as it does not contain any autocatalytic core. However, when we place this network in two compartments α and β coupled by a membrane permeable for A and A 2 B, a Type II core emerges (Fig. 4b).

The core identification indicates a possible setting for autocatalysis: U and V are chemostatted in α, and AB is chemostatted in β (Fig. 4c). The reaction involving U and V may in principle also take place in β, but it is not required for autocatalysis as it is not part of the type II core. In the present example, U and V are absent in β.

We now apply our viability analysis to this autocatalytic network in the presence of degradation reactions (r 6 and r 7 in Fig. 4c). Let us introduce a characteristic rate k for reactions r 2 , r 4 and r 5 , a degradation rate k d for reactions r 6 and r 7 and an exchange rate k ex for reactions r 1 and r 3 . Figure 4d shows the extinction probability as a function of the degradation rate k d and exchange rate k ex , both normalized by k. To overcome the degradation threshold (P ex < 1), the ratio k d /k must lie above a certain threshold (vertical black dotted line in Fig. 4d), and the rate of exchange k ex should outpace the rate of degradation (black slanting dotted line in Fig. 4d). Stochastic simulations (SI Appendix Sec. 8 and Fig. S9) confirm autocatalytic growth in these conditions.

+ + + + α β r 3 r 4 r 1 c) a) d) k k ex P ex r 2 r 5 r 6 r 7 A A AB AB A 2 B A 2 B α β b) Type II Type II U V V U U V A AB A 2 B
In this example of coupled compartments, compounds are no longer restricted to one role: AB is an autocatalyst in α and a feedstock in β. Chemical reactions are no longer restricted to one direction: The reaction used for reproduction in α is reused in β to provide the missing step to close the cycle. Such multicompartment autocatalysis is however more general. For instance, a single reaction A ----B + C can give rise to Type III motifs, given three compartments coupled by selective exchange as detailed in SI Appendix Sec. 9 and Fig. S9.

DISCUSSION

We presented a theoretical framework for autocatalysis based on stoichiometry, which allows a precise identification of the different forms of autocatalysis. Starting with a large stoichiometric matrix, we provide criteria for reaction network motifs that allow allocatalysis and autocatalysis. A detailed analysis of the graph structure contained in these reduced stoichiometric matrices reveals that they contain only five possible recurrent motifs, which are minimal in the sense that they do not contain smaller motifs. Fundamental modes of production of minimal autocatalytic cores are encoded in the column vectors of the inverse of the autocatalytic core submatrix. Autocatalytic cores are found to have a single reactant species for each reaction. This means that autocatalytic networks require the availability of certain chemical species in their cores to operate properly, but also implies that the proper functioning of an autocatalytic network will guarantee the stable supply of certain products, a definitive advantage when these products are key enzymes or metabolites.

We identified these minimal motifs in known examples of autocatalysis such as the formose reaction, central metabolic cycles, the GARD model and RAF sets. Autocatalytic cores also provide a basis for algorithms to identify these recurring autocatalytic motifs in large chemical networks [24,44,45], as has been done for gene regulatory networks [46]. In this way, we may be able to break the complexity of large chemical networks into smaller, more manageable structures [47]. Additionally, autocatalytic cores are the building block of evolution in prebiotic chemistries [36], thus their identification paves the way to a systematic exploration of the possible modes of chemical evolution [48,49].

Autocatalytic motifs provide different degrees of robustness, which we evaluated using the notion of viability. Viability can be computed as a survival probability in an appropriately defined branching process. This approach is generally applicable to autocatalytic models upon identification of their cores, highlighting the interest of a unified framework. Viability results from a competition between reactions that produce autocatalysts and side-reactions such as degradation. This is intimately related to the 'paradox of specificity' [19,40]: autocatalytic motifs are more likely to be found in large networks with many different chemical components engaging in many different reactions, but putting many components together favors side-reactions, leading to extinction.

Multicompartment autocatalysis introduced here offers a way around this problem: coupled compartments effectively enlarge the number of species without requiring new reactions. In multicompartment autocatalysis, cycles rely on the environmental coupling of reaction networks, which allows access to conditions unattainable in a single compartment. In this way, autocatalysis can emerge from reaction schemes as simple as a bimolecular reaction, provided certain semi-permeability conditions are met for the exchange of compounds between compartments. In the example shown here (Fig. 4), this allowed us to reuse the compounds and reactions to complete autocatalytic cycles. The principle is more general, however: autocatalysis may also emerge from coupling phases with physical-chemical conditions conducive to different reactions, as observed in liquid-solid [50], solid-gas [51] interfaces. Liquid-liquid interfaces in cellular organization and multiphase coacervates [43] are a promising places to further explore such principles.

Overall, our framework shows that autocatalysis comes in a diversity of forms and can emerge in unexpected ways, indicating that autocatalysis in chemistry must be more widespread than previously thought. This invites to search for further extensions of autocatalysis, which provides new vistas for understanding how chemistry may complexify towards life [52].

Theoretical methods and derivation of results are detailed in the appendices comprising the following sections: A) Terminology and definitions, B) derivation of autocatalytic cores from Graph theory, C) their chemical interpretation and D) application to formose, autoinduction, metabolic cycles, chemical amplification, RAF sets, GARD E) branching process derivation and determination of P ex . F) determination of P ex for Fig. 3. G) determination of P ex for Fig. 4d. H) stochastic simulations. I) autocatalysis from one bimolecular reaction and 3 compartments.

To promote the consistent use of terminology, IUPAC committees establish recommendations which serve as a basis for the Compendium of Chemical Terminology (The Gold Book), of which some relevant entries are reproduced here. Comments in italic have been added for clarity.

Chemical Reaction: A process that results in the interconversion of chemical species. Chemical reactions may be elementary reactions or stepwise reactions. It should be noted that this definition includes experimentally observable interconversions of conformers. Detectable chemical reactions normally involve sets of molecular entities as indicated by this definition, but it is often conceptually convenient to use the term also for changes involving single molecular entities (i.e. 'microscopic chemical events').

Catalyst: A substance that increases the rate of a reaction without modifying the overall standard Gibbs energy change in the (net) reaction; the process is called catalysis. The catalyst is both a reactant and product of the (catalyzed) reaction. The words catalyst and catalysis should not be used when the added substance reduces the rate of reaction (see inhibitor). Catalysis can be classified as homogeneous catalysis, in which only one phase is involved, and heterogeneous catalysis, in which the reaction occurs at or near an interface between phases. Catalysis brought about by one of the products of a (net) reaction is called autocatalysis. Catalysis brought about by a group on a reactant molecule itself is called intramolecular catalysis. The term catalysis is also often used when the substance is consumed in the (net) reaction (for example: base-catalysed hydrolysis of esters). Strictly, such a substance should be called an activator.

Autocatalytic Reaction: A (net) chemical reaction in which a product (or a reaction intermediate) also functions as a catalyst. In such a reaction the observed rate of reaction is often found to increase with time from its initial value.

Allocatalysis

We refer to allocatalysis as the form of catalysis in which, at the end of a catalytic cycle, the catalyst(s) have not changed in number. By their equal participation in either direction, allocatalysts will thus drop out of the net reaction. Some authors refer to autocatalysis as homocatalysis and allocatalysis as heterocatalysis, which is a lexicologically consistent choice of terms that express an opposition (same vs different). This opposition between same and different is e.g. found in the IUPAC terminology for a homogeneous catalysis (occurring in the same phase) and heterogeneous catalysis. For the IUPAC recommended terminology 'autocatalysis', a consistent choice that expresses this opposition is 'allocatalysis' (self vs other).

Remarks a. Directionality of autocatalysis

In an allocatalytic reaction, catalysts are produced and consumed in equal amounts, and the term does not distinguish between the direction in which the reaction proceeds (catalysis equally increases the rates of 'forward' and 'backward' reactions by the thermodynamic criterion). In contrast, the term autocatalysis only applies in one direction, due to the requirement of having catalysts be the product of the net reaction. Consequently, in writing the simplest reaction balance

A + B ----AB ----2 B, (S1)
the definition applies when the net reaction A ----B exhibits acceleration, owing to B being a catalyst and product of the net reaction. In the opposite sense, B ----A, B is no longer a product, but it can still accelerate the reaction. This case is typically referred to as 'reverse autocatalysis'. A simple experimental example is the catalytic disproportionation of water adsorbed on a copper surface

H 2 O(ads) [? ]. 2H 2 O(ads) 1 --H 2 O • OH(ads) + H(ads), (S2) H 2 O • OH(ads) 2 --H 2 O(ads) + OH(ads). (S3)
Performing r 1 and r 2 from left to right, the catalyst H 2 O is consumed, yielding reverse autocatalysis for the net reaction:

H 2 O(ads) ----H(ads) + OH(ads). (S4)
In the opposite sense, the stoichiometric conditions for autocatalysis are obtained:

H(ads) + OH(ads) + H 2 O(ads) ----2H 2 O(ads). (S5) b. Inhibition
A compound is a catalyst in the context of a particular experimental condition where it accelerates the rate of a reaction. A change of those conditions may change this label, e.g. the compound may become an inhibitor instead. Consider

A + B 1 --AB, (S6) A + E 2 --AE, (S7) AE + B 3 --AEB, (S8) AEB 4 --E + AB. (S9)
For a rate acceleration to occur for the net reaction A + B ----AB, steps 2,3,4 must together proceed faster than reaction 1. Template-assisted ligation of DNA or RNA provides an example where this may not be the case: below a critical annealing temperature the detachment of the product (here, AB in step 4) becomes increasingly slower for longer strands, and the template becomes an inhibitor.

c. Autonomy and siphons

The concept of autonomy, as defined in the main text, is closely related to the concept of a siphon in Chemical Reaction Networks (CRN) theory [? ]: A Siphon Σ is a subset Σ ⊂ S of all species S, for which each reaction that has a species in Σ as a product, has at least one of its reactants in Σ.

In this definition, a reaction can be irreversible in a mathematical sense: the reverse reaction does not exist. For a reversible CRN, a reverse reaction does exist, and the siphon definition must apply to the forward and backward direction, thus becomes equivalent to autonomy. A reaction must then imply both one or more products and one or more reactants from Σ (siphon reactions), or none at all (external reactions). Note that autonomy is less restrictive than the conditions posed in Barenholz et al. [24], where in addition species must be both reactant and product of a reaction. This last conditions is a proved as a consequence for minimal structures in our choice of formalism (see below).

Appendix B: Mathematical derivation of autocatalytic cores 1. Reaction graph definitions

Reaction graphs described below correspond to weighted directed hypergraphs without self-loops in the language of graph theory. In this section, the word cycle is used in the sense of graph theory (see below), not in the sense of reaction cycle used to denote right null-vectors of the stoichiometric matrix. In the following, letters used for scalars indicate positive numbers, and cycles and paths are understood as directed.

Definition 1. A reaction graph H is a triplet (S, R, M )

where {s 1 , ..., s n } is the species set, R = {r 1 , ..., r m } is reaction set, each r j being an ordered pair (X j , Y j ) of non-empty and non-intersecting subsets of S respectively called reactants and products, and M is the stoichiometric matrix with coefficients m ij . m ij = 0 when the species s i does not participate in the reaction r j , m ij < 0 when species s i is a reactant of r j , and m ij > 0 when species s i is a product of r j .

The stoichiometric matrix M contains all the information about the hypergraph, the column of M corresponding to reactions, and the rows to species.

Definition 2. A subgraph of H = (S, R, M ) is a triplet H = (S , R , M )
where S is a subset of S, R is a set of reactions whose reactants and products are in S and intersect the reactant set and product set of a reaction in R, with corresponding stoichiometric coefficients M . Definition 3. A reaction graph is square if it has the same number of reactions and species. Definition 4. A directed path is a sequence of alternating species and reactions, all reactions and species being distinct, where species which precede and succeed a reaction are respectively a reactant and a product of it. A path is a minimal path if it is not possible to form a path starting and ending at the same species using a strict subset of its reactions. A path is semi-open if it either starts or ends with an edge. Definition 5. In a directed path, an edge has a backbranch if one of its products is a species located upstream in the path, and this product is called a back-product. An edge has a forward-branch if one of its reactants is a species located downstream in the path, and this product is called a forward-reactant. Definition 6. A cycle has an identical definition as a path, except that the first and last species are the same species. A cycle is minimal if it is not possible to form a cycle with a subset of its species and reactions. Definition 7. A species S is the solitary reactant (product) of a reaction if S is the only reactant (resp. product) of this reaction, otherwise it is a co-reactant (resp. coproduct). Definition 10. In a square graph, a perfect matching is a bijection between species and reactions. It corresponds to all pairs (i, σ(i)) i=1...N , where i is the index of a species, σ is a permutation of [1...N ], and σ(i) is the index of a reaction. A perfect matching of a graph, or a subgraph, G is called a G-matching.

Remarks:

• Consider reactions E = ({a}, {b, c}) and F = ({b}, {c}); a -E -b -F -c is a path, but it is not minimal because it contains a -E -c. Simple paths are necessarily minimal (Fig. 5a).

• Consider E = ({a}, {b}) and F = ({b}, {a, c}); a -E -b -F -c is a minimal path and F has a backbranch (Fig. 5a). In particular, it contains a cycle a -E -b -F -a.

• Simple paths are exactly cycle-free minimal paths.

• More generally, a minimal path can have backproducts and forward-reactants (Fig. 5a).

• A hypergraph cycle can have reactions connecting several of its species, but a minimal cycle cannot. Thus, a minimal cycle only contains simple paths (it is identical to a cycle in a regular graph).

• Cycles are square.

Relationship with linear algebra

x 0 denotes a real vector with only strictly positive coordinates.

Definition 11.

A matrix M is productive if there exists a real vector γ such that M . γ 0. Equivalently, M intersects the strictly positive orthant IR n >0 . Definition 12. A matrix is autonomous if its columns all contain a strictly negative and a strictly positive coefficient. Definition 13. A minimal cycle is weight-symmetric if the product of the absolute values of the stoichiometric coefficients associated with its reactants equals the product of the stoichiometric coefficients associated with its products. Otherwise, the cycle is weight-asymmetric.

Remarks:

• Autonomous matrices are stoichiometric matrices of reactions systems such that every reaction has at least one reactant and at least one product.

• The Leibniz formula for the determinant shows that the non-zero terms of det(M ) exactly correspond to the products of stoichiometric coefficients m i,σ(i) determined by each possible perfect matching of the graph.

• A minimal cycle C has exactly two perfect matchings, which correspond to the matching of its reactions with their solitary reactants and solitary products respectively. det(C) = 0 if and only if C is weightsymmetric.

• Finding cores is equivalent to finding minimal matrices which are productive and autonomous. In this section, we denote M ∈ IR n,m the stoichiometric matrix of the graph H, where n is the number of rows (species) and m the number of columns (reactions). If M is productive, we denote γ a vector such that M.γ 0. Productivity of M is indifferent to the sign of its columns as the sign of the coefficients of γ is not constrained. Therefore, we choose the convention that any productive vector γ is positive, up to taking the opposite for some columns of M .

Remarks:

• M invertible implies M productive, as the image of M is then the full space IR n,m , which contains the strictly positive orthant.

• Below, species or reactions 'can be removed' is understood as 'can be removed while preserving productivity and autonomy'. Being able to remove a row (a species) or a column (a reaction) contradicts the minimality of a core.

• Removing columns (reactions) preserves autonomy, but not necessarily productivity.

• Removing rows (species) preserves productivity, but not necessarily autonomy.

• A row corresponding to a species that is always a co-reactant or a co-product can be removed without affecting autonomy (every column still contains positive and negative coefficients) and productivity.

Proposition 1. In a core, every species is both a reactant and a product.

Proof. Every species must be produced, otherwise it would not be possible to find a positive γ verifying M.γ 0. Now, suppose a species S is never a reactant and only a product. All reactions such that S is their only product can be removed without affecting the productivity of other species. In the resulting graph, either S is not the product of any reaction anymore, or S is only a co-product of the remaining reactions. In both cases, S can be removed without affecting autonomy. Proposition 2. A core is square, invertible, every species is the solitary reactant of a reaction, and is reactant for this reaction only.

Proof. Consider M ∈ IR n,m a productive stoichiometric matrix with n species and m reactions, with rank k. Obviously k ≤ m, n. We must also have m = k, otherwise a column could be removed while preserving the image of M , thus productivity. Additionally, every species must be the solitary reactant of at least one reaction, otherwise it could be removed. This implies m ≥ n. Overall, k = m ≥ n ≥ k, so that k = m = n, meaning that M is square and invertible. As every species S is the solitary reactant of at least one reaction R and m = n, the species is a reactant for only R.

Remark: Property 2 implies that M can be re-arranged in such a way that it has a strictly negative diagonal, and only non-negative off-diagonal coefficients. Proof. Consider a species x 0 of M . Below, we recursively construct sets D k = {x 1 , ..., x k } of increasing cardinal, such that every x i is downstream x 0 , until k = n -1, implying that for any species y = x, there exists a path from x to y. We denote R(S) the only reaction with reactant S, which is well defined by Property 2.

Step 1:

We take x 1 as a product of R(x 0 ). The definitions below are generalizations of the notion of ear decomposition in regular graphs (Fig. 5b).

Definition 15. A hyper-ear is a hypergraph comprising a minimal cycle C, called the base cycle, such that C has a reaction with a product x outside C, and a minimal path P starting at x, such that its last reaction, R, has a product in C, R being the only P-reaction to have a product in C. A proto-ear has a similar definition, but where C is a simple path (called the base path) instead of a minimal cycle.

Description of proto-ears and hyperears -In a hyper-ear or a proto-ear, any C-reaction can have x as a product, and any C-species can be the product of R, the last reaction of P (Fig. 5b). We denote by u a species which is the product of R, v any C-species which is the reactant of a C-reaction which produces x, '-' a simple path (including the empty path), o a simple path comprising at least one species, o being a subcase of the '-' category. By convention, a uv motif can correspond to a single same species which is both a R-product and a reactant for x production. Any proto-ear falls into a class described by a chain made of the symbols u, v, o, and '-', as soon as it contains at least a u and a v. Reciprocally, any such chain represents one or more proto-ears. Any hyper-ear is obtained by cyclic closure of the base path of a proto-ear. Such closure is denoted by the * symbol at the beginning and the end of the chain.

Theorem 1. Cores are of one of the following types:

• TYPE I: a weight-asymmetric minimal cycle;

• TYPE II: a cycle comprising all species as solitary reactants, and one or more weight-symmetric subcycles without intersection between them;

• TYPES III-IV-V: one of the three types of hyperears, any subcycle of which is weight-symmetric:

(III) * u -vo * ; (IV) * u -u -v * ; (V) * u -v -u -v * .
Proof. Core type are represented in Fig. 5c. H denotes the reaction graph and M its stoichiometric matrix.

SUFFICIENCY

Step 1: All types have a non-zero determinant. TYPE I: Invertibility is a direct consequence of the remark on the determinant of minimal cycles.

TYPE II: Consider a minimal weight-symmetric subcycle C of H. Any H-matching consistent with a Cmatching can be associated to another H-matching obtained by reversing the C-matching. Thus, det(H) has the form det(C).α + β, where the first term comprises all C-consistent H-matchings. As det(C) = 0, it suffices to show that β is a single non-zero term. By Property 3, C must have a reaction R 1 with a product x 1 outside of C. As by definition of TYPE II, subcycles are non-intersecting, R 1 must be the only C-reaction with a product outside C, and x 1 must be the only product of R 1 outside C. All non-zero terms in β require matching R 1 to x 1 , otherwise R 1 would be matched with a C-species, which would impose a C-matching (C being a minimal cycle), contradicting the definition of β. We now order the indexes k following the downstream order of the x k species along H, and show recursively that β corresponds to the H-matching where every R k matches x k : (step 1) By construction, R 1 matches x 1 . (step k) Suppose R k-1 matches x k-1 . Either R k is a simple reaction, and, as its reactant x k-1 is already matched, R k necessarily matches its only product x k . Or R k has a back-branch forming a minimal cycle. However, cycles are non-intersecting, so that the back-product x i of R k is necessarily downstream x k and upstream x k-1 , thus x i is already matched by the recursion hypothesis. Thus, R k can only match x k .

TYPES III-IV-V: Consider the stoichiometric matrix M with non-negative off-diagonal coefficients:

M =   -1 c e a -1 f b d -1   (S1)
We have det(M ) = -1 + df + ac + ade + bcf + be. Strict subcycles must be weight-symmetric, as otherwise, the minimality of M would be contradicted. Consequently, when their factors are both non-zero, the products df , ac, and be must be equal to 1. By contracting all simple reactions and multiplying the columns of M if necessary so that solitary reactant stoichiometries are all normalized to -1, TYPE III corresponds to d = f = 0 and a, b, c, e > 0, TYPE IV to f = 0 and a, b, c, d, e > 0, and TYPE V to all coefficients strictly positive. In all these cases, we have det(M ) > -1 + ac + be = 1.

Step 2: All types are minimal TYPE I: Removing any subset of species or reactions would result in an acyclic graph, contradicting Property 5.

TYPE II: Removing any subset of species or reactions either leads to a hypergraph where a subset of species is upstream of the rest, contradicting Property 6, or to a non-invertible minimal cycle.

TYPE III-V: Removal of any set of reactions (and a fortiori of species, given that they all are solitary reactants) leaves at most one cycle in the graph, the latter being non-invertible.

NECESSITY

By Property 5, there exists a minimal cycle C in H. Either C is weight-asymmetric, then H = C is of TYPE I. Or C is weight-symmetric. Then, by Property 3, there is a C-reaction with a product x outside C. By Property 4, there exists a path P from x to any species in C (Fig. 5d). We can take P minimal and such that only its last reaction has a product in C, thus forming a hyper-ear * mu -vm * where m symbols stand for any other hyperear motif. Below, we show that hyper-ears are either one of the types II-V, or contain a TYPE II, III or IV core as a subgraph, overall demonstrating that H is necessarily a hyper-ear of TYPE II-V.

Before continuing the proof of the theorem, we show an additional property based on the sufficiency of TYPE I and TYPE III.

Proposition 7. In a core, a minimal path can have backbranches but no forward-branch, and the cycles formed by back-branches are non-intersecting.

Proof. Consider a minimal path P in a core H. Forwardbranches require reactions to have multiple reactants, contradicting Property 1. Therefore, P is either a simple path or it has back-branches. Suppose there exists two back-products y and z of the respective reactions r y and r z , respectively forming intersecting cycles C y and C z (Fig. 5e). Without loss of generality, y and z can be taken closest, y upstream of z, and such that there is no other cycle nested within C y and C z than possibly themselves. C z is necessarily weight-symmetric, as it would otherwise contradict the minimality of H. There are two cases. Case 1: r z is upstream r y (C z is nested within C y ). Call P the path from the product of r z to r y then y then z. Then C z and P form a TYPE III core. Case 2: r y is upstream r z (C z and C y are entangled). Call P the path joining y to z. Then C z and P form a TYPE III core. Overall, we have shown that the existence intersecting cycles along a minimal contradicts the minimality of H.

We now go back to the proof of the theorem.

Proof. Consider the minimal path P of the hyper-ear, where P starts at x. By Property 7, P has no forwardbranch.

Case A (TYPE II): Suppose P has one or more back-branch. Consider a reaction R of P which has a back-branch, with back-product y and forming a cycle C . C is necessarily weight-symmetric by minimality of H. Call x the product of R outside C . Consider the path P starting at x , then going downstream of P until it reaches a shortest subpath u -v of C, then to x, and finally from x to y in P. By Property 7, all cycles along P are nonintersecting, and by construction, C is non-intersecting with the cycles of P. Consequently, P only has nonintersecting cycles, so that C and P form a TYPE II motif. Thus, the base cycle of the hyper-ear cannot have species outside u -v, as otherwise the core made of C and P would be a strict subset of H, contradicting its minimality. This shows that H is a * u -v * hyper-ear with a minimal path comprising non-intersecting cycles, the latter being necessarily weight-symmetric. Thus H is a TYPE II core.

Case B (TYPES III-V): Suppose P has no backbranch, in other words P is a simple path.

Subcase B.1 (TYPE III): Suppose H a hyper-ear which has only one u and one v symbol, thus a * u -v - * hyperear. The case * u -v * with a simple path is covered by the definition of TYPE II. If thesymbol does not represent an empty path, then it comprises at least one species, so that H is a * u -vo * hyper-ear with simple path P, corresponding to TYPE III. Subcase B.2 (TYPE IV): Suppose the hyper-ear comprises one u or v symbol in addition to the u -v motif. We first note that u -u -v and u -v -v proto-ears are isomorphic to * u -v - * , which falls into the categories of TYPE II or TYPE III cores. Therefore, in a core, a hyper-ear containing two successive u or v symbols is necessarily of the form * u -u -v * or * v -v -u * , or their cyclic permutations, as any additional symbol would allow to find a strict subgraph u -u -v or u -v -v forming a core, contradicting the minimality of H. Furthermore, * u -u -v * and * v -v -u * hyper-ears are isomorphic. Indeed, the matrices of their reduced forms correspond to the matrix shown in the SUFFICIENCY section of the theorem, where exactly one coefficient is set to zero.

Thus, these motifs as well as all their cyclic permutations, fall into the category of TYPE IV. Subcase B.3 (TYPE V): Subcase B.2 imposes that any motif containing four or more u or v symbols must alternate u and v symbols. Any motifs with five or more u or v symbols contains, up to permutation, a u -v -u -v proto-ear as a strict subgraph. However, the latter is isomorphic to TYPE IV cores, which implies that * uv -u -v * (TYPE V) is the only motif in this class which does not contradict the minimality of H.

Appendix C: Chemical interpretation of autocatalytic cores

We remind that the notion of 'graph cycle' (closed successions of nodes and edges) differs from the notion 'reaction cycle' (right null vectors of the stoichiometric matrix). We employ 'reaction cycle' and 'graph cycle' to distinguish these notions.

By the minimality of autocatalytic cores, an autocatalytic motif is either a core, or it contains one or several cores. If P is a property of cores, then for any autocatalytic motif, either it verifies P , or it contains an autocatalytic motif verifying P . In particular,

• by Property 1, every autocatalytic motif contains an invertible autocatalytic motif;

• by Property 2, every autocatalytic motif contains an autocatalytic motif such that every product is also a catalyst of the reaction, or equivalently, such that every species appears on both side of the total chemical equation.

The five categories of cores are schematically represented in Fig. 5c. The convention of representation are as follows. The edges of the graph correspond to reactions and the yellow nodes to species. Reaction have two sides, the reactant side and product side, and each side of the edge representing the reaction connects to one or several nodes representing the species, with a stoichiometry for each connection. In the Type I core, the edge from the top node to the bottom node ends with a fork, the two ends of which connect to a single bottom node (5c). This means that the bottom species is produced with a stoichiometry of 2 by this reaction. For simplicity, we have represented the connection between all other edges and nodes without fork. However, in all generality, any connection between an edge and a node could be a fork (e.g. could be of stoichiometry >1), as soon as the rules on graph cycle symmetry are respected, as explained below. Forks also appear in Types II-V, but where they connect to two distinct nodes, which are two distinct product species. The orange squares indicate that the reaction represented can be replaced by a chain of reactions, in which the first reaction in the chain has only one reactant and the last reaction in the chain a single product'.

The mathematical derivations are done without constraint on the number of reactants or products a reaction step can have. However, in a chemical system, elementary reactions are in principle either unimolecular or bimolecular. This restricts the range of possible core graphs. This restriction operates only at the level of the stoichiometry of each single reaction, not at the level of the overall structure of cores. Indeed, even without invoking the restriction on bimolecularity, the mathematical derivation results in cores such that every reaction step involves at most a single species as reactant and at most two species as products, so that bimolecularity can always be respected, provided rules on the stoichiometry of each connection.

In a chemical autocatalytic core, an edge with a fork connecting to two distinct nodes (e.g. a reaction with two distinct product species) must have a stoichiometry 1 for each connection in order to respect bimolecularity. If a reaction has only one product species, then its stoichiometry can be 1 or 2, as soon as it respects rules on cycle stoichiometry, which we detail now.

Consider a simple graph cycle C, simple meaning that every reaction has a single reactant species and a single product species. Note a i (resp. b i ) the stoichiometry of the reactant (resp. product) of reaction i.

If i |a i | = i |b i |,
we say that C is weight-asymmetric, otherwise it is weight-symmetric. Weight-asymmetric simple graph cycles are an example of a graph cycle which has no reaction cycle. Weight-symmetric graph cycle taken in isolation have a reaction cycle (their determinant is zero as shown in the derivation of cores above), thus they correspond to allocatalytic cycles in the context of a larger reaction graph where the reactions of the graph cycle consume and/or produce species outside of its own species set. The most classic example of an allocatalytic cycle is represented in Fig 5f, where reactant S is provided from the environment, binds to catalyst E to form complex ES converted into EP, finally dissociated into E which is recycled, and product P.

Type I cores are weight-asymmetric simple graph cycles, as C in 5c . Consequently, in types II-V, any simple graph cycle must be weight-symmetric (these graph cycles corresponding reaction cycles), otherwise the core would contain a Type I core, contradicting its minimality. For example in 5c: in Type II, C must be symmetric, but this does not apply to C because it is not even a simple cycle; in types III-V, C, C and C must be symmetric. Type II cores consist of a large graph cycle (C in Fig 5c ) comprising smaller graph cycles embedded within it (for example C in Fig 5c). Given the above, each of these smaller cycles must be weight-symmetric, and obeys the definition of an allocatalytic cycle. The Type II category includes circularly closed successions of such allocatalytic cycles, where the product of one allocatalytic cycles is the catalyst of the next allocatalytic cycle, which is a typical example of autocatalytic set. In addition however, Type II allows intermediate non catalyzed reaction steps.

Notably, every core follows the basic structure repre-sented on Fig 5d, comprising a base cycle C and a minimal path (in the sense of reaction hypergraphs, see former section) starting from a reaction fork and joining back to a node of C. This structure should enable a systematic algorithmic search for autocatalytic motifs in large stoichiometries.

Appendix D: Examples and application of the stoichiometric criteria

Example of autocatalytic submotif

We illustrate the concepts of autocatalytic submotif and the properties demonstrated above on a toy model for the Formose reaction, to which we have added one auxiliary reaction C 3 ----D 3 , so that we obtain a stoichiometric matrix

1 2 3 4 ν ν ν = C 1 C 2 C 3 D 3 C 4      -1 -1 0 0 -1 0 2 0 1 -1 0 -1 0 0 0 1 0 1 -1 0      C 1 + C 2 1 --C 3 C 1 + C 3 2 --C 4 C 4 3 --2 C 2 C 3 4 --D 3 (S1)
this full matrix has a left nullvector l l l = (

we have a mass-like conservation law

L = n C 1 + 2n C 2 + 3n C 3 + 3n D 3 + 4n C 4 . (S2)
Thus, this matrix does not correspond to an autocatalytic motif. Upon removing C 1 (then considered as a feedstock molecule), we obtain an autocatalytic matrix ν ν ν * 1 2 3 4

ν ν ν * = C 2 C 3 D 3 C 4    -1 0 2 0 1 -1 0 -1 0 0 0 1 0 1 -1 0    C 2 1 --C 3 C 3 2 --C 4 C 4 3 --2 C 2 C 3 4 --D 3 (S3)
This matrix is autonomous and invertible with inverse:

ν ν ν -1 * = 1 2 3 4    1 2 2 2 1 1 1 2 1 1 1 1 0 0 1 0    = (g g g C 2 , g g g C 3 , g g g D 3 , g g g C 4 ). ( S4
)
The columns of the inverse are reaction vectors associated to a given species, of which they produce one extra unit. We will refer to them as elementary production modes.

Among the production modes, g g g D 3 = (2, 1, 1, 1) T is the sole vector which uses the auxiliary reaction, and D 3 only ever occurs as a product

2 C 2 + 2 C 3 + C 4 g g g D 3 -------- -g g g D 3 2 C 2 + 2 C 3 + C 4 + D 3 . (S5)
If we now consider Γ Γ Γ = g g g C 2 +g g g C 3 +g g g D 3 +g g g C 4 , we obtain a reaction balance such that every species has a net increase:

Γ Γ Γ : 7 C 2 + 6 C 3 + 4 C 4 --→ 8 C 2 + 7 C 3 + 5 C 4 + D 3 .
(S6) ν ν ν * is not minimal as D 3 is only produced, and does not participate as a catalyst. Indeed, Property 2 implies the existence of an autocatalytic submotif without species which does not participate as a catalyst. An autocatalytic submotif is obtained by removing the reaction

C 3 ----D 3 , (S7)
and the species D 3 , to obtain the autonomous submatrix ν ν ν

1 2 3 ν ν ν = C 2 C 3 C 4   -1 0 2 1 -1 0 0 1 -1   C 2 1 --C 3 C 3 2 --C 4 C 4 3 --2 C 2 (S8)
The matrix ν ν ν is also invertible

ν ν ν-1 = 1 2 3   1 2 2 1 1 2 1 1 1   = (g g g C 2 , g g g C 3 , g g g C 4 ). ( S9
)
The columns of ν ν ν-1 correspond to reaction vectors, whose application yields one net copy of the corresponding molecule, e.g. ∆n k = (ν ν ν • g g g (k) ) k = 1. This is illustrated for C 2 in Fig. 6.

The individual replication cycles are

C 2 + C 3 + C 4 g g g C 2 -------- -g g g C 2 2 C 2 + C 3 + C 4 (S10) 2 C 2 + C 3 + C 4 g g g C 3 -------- -g g g C 3 2 C 2 + 2 C 3 + C 4 2 C 2 + 2 C 3 + C 4 g g g C 4 -------- -g g g C 4 2 C 2 + 2 C 3 + 2 C 4
We can construct Γ Γ Γ = g g g C 2 + g g g C 3 + g g g C 4 , which leads to the overall reaction

5 C 2 + 4 C 3 + 3 C 4 Γ Γ Γ ---- -Γ Γ Γ 6 C 2 + 5 C 3 + 4 C 4 . (S11)
We see here that every species has a net production and is on both sides of the balance, thus participates as a autocatalyst. Furthermore, ν ν ν is a Type I core and consequently does not contain any smaller autocatalytic submotifs.

Autoinduction

The concept of autoinduction was put forward by D. Blackmond [21], to distinguish between autocatalysis that is mediated by external catalysts (i.e. not part of the autocatalysts that are reproduced) called 'autoinduction', and autocatalysis that functions without the aid of external catalysts (i.e. all catalysts are autocatalysts). We thus obtain a hybrid of pure allocatalysis and autocatalysis, for which the simplest example would be

A + B + E ----2 A + E. (S12)
The IUPAC definitions impose that autoinduction qualifies as autocatalysis. It follows then from our framework that we can find autocatalytic cores in autoinduction networks, which is indeed confirmed in Fig. 7 for the two types of autoinduction that have been proposed [21].

The concept of 'autoinduction' addresses a notion of self-sufficiency (also encountered in RAF sets, Sec.D 5): external allocatalysts become essential to succesful autocatalysis, yet they are not reproduced. Depending on the context, they could be seen as part of the environment, in the same sense as essential feedstock species. Examples of autoinduction occur in autocatalytic metabolic networks (with locally allocatalytic enzymes) or e.g. the formose reaction which is often catalyzed by base and divalent metal ions.

The presence of external allocatalytic cycles does not add new cycles to the autocatalytic core. A practical consequence is that one can write catalyzed reactions very compactly for the core, while still maintaining nonambiguity, which we make use of in our analysis of metabolic cycles.

Metabolic cycles

At least two metabolic cycles are known to be autocatalytic. In our analysis of autoinduction, we pointed out that the core does not contain external allocatalysts (here: enzymes). Written purely in terms of autocatalysts, we find a Type II autocatalytic core for the reverse Krebbs cycle (Fig. 8a). For the Calvin cycle depicted in Fig. 8b, we identify three Type I cores (two structurally equivalent, differing only in the reaction chosen to link the same core species) and 4 of Type II.

Chemical amplification

Chen et al. [38? ] advanced a general strategy to achieve self-amplifying behavior, demonstrated by the encapsulation of the reactive compound DCC in a cavitand (indicated by an oval around the molecule in EDCC, Fig. 9), which they referred to as chemical amplification. This phenomenon is consistent with the IUPAC definition of autocatalysis, which we will now illustrate by finding its autocatalytic core.

Let us first assess a reaction network formed by two reactions proposed by Chen et al, in absence of the cavitand complexes. In that case, applying our formalism readily reveals that the network cannot exhibit autocatalysis: there are no autocatalytic submatrices.

1 2 ν ν ν = DCC X I1 Y DCU Z        -1 0 -1 0 1 -1 0 -1 0 1 0 1        DCC + X 1 --I1, I1 + Y 2 --DCU + Z. (S13)
Now, let us add the species EDCC and two exchange reactions with DCU and Z that liberate DCC, i.e.

DCC + X 1 --I1, (S14) I1 + Y 2 --DCU + Z, (S15) EDCC + DCU 3 --DCC + EDCU, (S16) EDCC + Z 4 --DCC + EZ, (S17)
for which we have the matrix

1 2 3 4 ν ν ν = DCC X I1 Y DCU Z EDCC EDCU EZ              -1 0 1 1 -1 0 0 0 1 -1 0 0 0 -1 0 0 0 1 -1 0 0 1 0 -1 0 0 -1 -1 0 0 1 0 0 0 0 1              (S18)
This network admits an autocatalytic submatrix, and we obtain an autocatalytic core of Type III consisting of the species DCC, I1, DCU and Z, as can be seen in Fig. 9b. The new reagent EDCC serves as a feedstock compound, that allows to dispense new DCC, that can now serve as an autocatalyst. The generality of the mechanism follows from an exchange that could have been performed with different reactants than DCU and Z to yield an equivalent network, as shown in Refs [38? ].

Autocatalysis in RAF sets

A definition of Reflexively Autocatalytic Foodgenerated sets can be found in Ref. [20]: a set of reactions R is RAF if every reaction is catalyzed by at least one molecule involved in a reaction in R, and every reactant in R can be constructed from the food set f by successive applications of reactions from R. The philosophy behind the RAF set is that 'every reaction' in a RAF-set can be accelerated by the catalysts themselves. RAF sets do not perform autoinduction, except when members of the food set f within a RAF-set also serve as allocatalysts.

In the RAF-set formalism, reaction and catalysis are distinct mathematical objects. Graphically, RAF-sets are typically represented as bipartite graphs (10a), with nodes (white squares) corresponding to reactions, which connect to nodes (colored circles) which serve as reactants (entering bold arrow) and products (leaving bold arrow) via directed edges. A dashed arrow connecting to a reaction indicates that a species is a catalyst for a reaction. Within the RAF framework, the following terms are employed as distinct: i) autocatalytic reaction ("is a single chemical reaction for which one of the products also catalyzes the reaction"), ii) autocatalytic cycle ("is a sequence of reactions that, once completed, results in two (or more) copies of the molecule that was started with"such as the toy formose reaction), iii) autocatalytic set (or RAF set).

It is important to note that the RAF-set formalism and the distinctions above depend on the chosen level of coarse-graining of the description. In that description (see Fig. 10a) allocatalysts in the same allocatalytic cycle are represented by a single species. The combination of reactions that form the allocatalytic cycle is represented by a single dashed line, connecting to an uncatalyzed reaction. Here, there exists a level of description where a RAF set appears as autocatalytic reaction involving catalytic cycles. By comparing Figs. 10a, we see that the level of description is critical in assessing whether a network is a RAF-set or not. From the definition of a RAF-set, the detailed version of the network in Fig. 10a would not be a RAF, but the less detailed description next to it would be.

Another example is the toy formose reaction, where two different choices of coarse-graining yield two different description in the RAF framework. Even when neglecting the role of catalytic base and metal ions, the toy formose reaction without coarse-graining is not a RAF-set [29], since its individual steps are not catalyzed:

C 1 + C 2 ----C 3 , C 1 + C 3 ----C 4 , C 4 ----2C 2 . (S19)
Framed solely in terms of C 1 and one autocatalyst (e.g. C 2 , C 4 ) one could propose a description in which we are oblivious of these steps, and write it as an autocatalytic reaction in the sense of RAF sets, such as

2C 1 + C 2 ----2C 2 , 4C 1 + C 4 ----2C 4 .
(S20)

This would work as long as we consider them satisfactory levels of description, depending on the context (e.g. further reactivity of intermediates, separation of timescales): either of these hypothetical cases satisfy the RAF criteria. An interesting contrast occurs between RAF sets and our formalism: RAF-sets require a level of coarse-graining without intermediates to define catalysis in terms of single catalysts, which ensures a compact and unique description. Our formalism requires each catalytic cycle to have at least one intermediate to satisfy nonambiguity. Although less compact, the description is flexible: adding further steps is always possible and does not alter the conclusions, which must remain in accord with chemical definitions.

In describing catalysis in terms of uncatalyzed reactions, it becomes possible to formalize the underlying structure of catalysis in different models. We will now illustrate how this formally establishes that all autocatalysis in the GARD model qualifies as a RAF-set.

GARD model

GARD stands for graded autocatalysis replication domain [18], and is a model for the autocatalytic assembly of amphiphile assemblies (e.g. micelles). It describes micelles with a composition n n n = {n 1 , ..., n s }, that follow an evolution equation

dn i dt = k + i ρ i N -k - i n i   1 + 1 N s j=1 β ij n j   . (S21)
The surfaces are in contact with a reservoir, that contains species Z i at concentration ρ i and that can enter the surface, which has an area proportional to N . The incorporation happens with a base rate of k + i , but can be facilitated by other amphiphiles, for which the catalytic rate enhancement is characterized by β ij . A special ingredient in GARD is the division process, which splits a mature aggregate in two new ones, after achieving a maximal size.

We examine here stoichiometric mechanisms leading to (S21). Let us consider amphiphiles A and B, which can be in the micelle or reservoir, labelled I and II (see Fig. 10). To catalyze incorporation of the other, a complex is formed with a reservoir species, and subsequent dissociation takes place in the micelle

A I + B II ----[AB] I ----A I + B I , (S22) B I + A II ----[BA] I ----A I + B I . (S23)
This network corresponds exactly to the network discussed in Fig. 10a, where the products of two allocatalytic cycles serve as a catalyst for each other, a Type II core in our stoichiometric formalism. More generally, labelling amphiphiles as A k (k ∈ {1, 2, .., s}, the entry β ij encodes the contribution for the allocatalytic cycle

A II i + A I j c c c * ------ -c c c * A I i + A I j (S24)
where c c c * denotes a reaction vector for the catalytic cycle, for any description that verifies nonambiguity. When i = j, this is a Type I autocatalytic cycle (Fig. 10c). When i = j (cross-incorporation), Type II autocatalytic cycles are obtained, which are built up from n sequential allocatalytic incorporation steps and which end in the incorporation of the original amphiphile. Fig. 10 shows examples for n = 2 and n = 3. The importance of the allocatalysis step (S24) is graded by the entry β ij . For a given n-step autocatalytic motif to exist, we require

n k=1 β s k+1 s k > 0, s 1 = s 2 = .. = s n . (S25)
In practice, all β ij > 0, so all motifs exist in principle.

The starting point of GARD is (S21), i.e. a coarsegrained description in which allocatalysis can be described as a single step and a single allocatalyst. From Fig. 10, we see that we can obtain networks in GARD that would be RAF-sets in the RAF-framework by coarse-graining an incorporation cycle to convert it to catalysis in the RAF sense. The illustrated procedure extends to all autocatalysis in GARD, which is fully characterized by the continuation of the structures in Fig. 10 to their n-step Type II analogues. (S25) guarantees that each reaction in GARD is catalyzed. It follows that all autocatalysis in the GARD model can formally be treated as a RAF-set.

Interestingly, the RAF-set formalism treats catalysis as pertaining to chemistry in single phases, with the environment supplying food locally through rapid exchange. In GARD, we instead have phase-transfer catalysis between a bulk medium (II) and an interface (I). A species in the bulk then serves as the food. Once the exact same species enters the interface, however, it may cease to be abundant or have a fixed concentration due to rapid exchange. It may thus no longer have the properties ascribed to the food set in the bulk and should ipso facto be treated as a different species as described in the final section of the main text.

Appendix E: The extinction problem and fixation

In this section, we show that the extinction problem (finding the extinction probability at long times, P ex ) can be solved by a mapping to a branching process. We will first derive how in a system initially at steady-state perturbed with dilute autocatalysts, key statistical properties emerge: autocatalytic species can be treated as independent, and their environment as fixed. Extinction becomes exponentially less likely as the population size continues to grow, which means P ex can be determined in the dilute limit.

We then proceed by applying the framework to a variety of networks.

Context

We consider a reaction network, in a macroscopic system (letting N denote the amount of chemical species, let us say N > 10 23 , or similarly, let the system volume V be large) in a steady-state. For simplicity, let us first consider a CSTR (single phase, ideally mixed), with a residence time τ , corresponding to a uniform degradation rate k d . Now, we perturb the steady-state with a handful (O(1)) of new (not yet present in the system) autocatalysts {X k } = {X 1 , X 2 , ..., X s } that are part of the same autocatalytic core.

We consider that the population can grow due to catalysis by autocatalysts, and the population decays by degradation reactions and effective degradation. For example, outflow out of a CSTR is considered as a first-order degradation process:

X k → ∅ (S1)
The problem we wish to solve is the extinction problem: For an initial population of autocatalysts

{N X k } = {N X 1 , .., N X s } what is the probability P ex ({N X k })
, that, after a long time, the autocatalyst population goes extinct (∀k N X k = 0)?

Large system limit

Let us first note that we (deliberately) consider the initial stochastic kinetics in a large system, with a small number of autocatalysts, such that reactions among autocatalysts of the kind X k + X j --→ "P roduct(s)" (S2) are exceedingly rare and slow (the probability that a given X k molecule encounters another X j in a given time-frame scales with N X j /N ), where N X j is initially of the order 1. It follows survival of autocatalytic cycles requiring such reactions in the forward sens is hampered in large systems.

Reactions that are first order in terms of autocatalysts are not hampered

X k ----"P roduct(s)" (S3) X k + Y j ----"P roduct(s)" (S4)
where Y j is a feedstock compound that was already (abundantly) present in the system at a fixed molar fraction x Y j . The probability for one X k molecule to encounter a Y j does not change with N , as x Y j remains fixed (and macroscopic).

Note furthermore that, when autocatalysts are rare, reactions producing more autocatalysts are 'irreversible'

X k + Y j --→ X l + X m (S5)
in the sense that the reverse reaction is exceedingly more rare than the forward reaction.

Solving the quadratic equation yields two solutions

P ± ex = 1 ± (1 -2p 2 ) 2p 2 . (S18)
For our problem, the 'physical' solution is P + ex while p 2 ≥ 1/2. Beyond that point, P + ex > 1, while we require 0 ≤ P ex ≤ 1, so P - ex = 1 becomes the only physical solution.

P ex = 1 p2 -1, p 2 ≥ 1 2 , 1, p 2 < 1 2 , (S19)
The average number of descendants is 2p 2 , which means that below p 2 = 1/2 (the decay threshold) B 1 is on average replaced by less then one B 1 species and extinction is guaranteed.

A Branching Process for the same type I cycle

To illustrate that there is a variety of choices for the stochastic process under study, we will here consider an alternative choice for the simple Type-I cycle, which is more generally applicable. Noting that a single cycle is successful with probability p 2 , we now consider the number of successful cycles a single B 1 creates, knowing that at some point degradation intervenes with probability 1 -p 2 . A succession of k cycles preceding a failure will thus lead to k successors

B 1 p k -→ kB 1 .
(S20)

The probability to spawn k descendants, p k , is the probability of k successful Bernoulli trials followed by failure, and follows a geometric distribution

p k = p k 2 (1 -p 2 ). (S21)
Injecting this distribution in (S9), we find a geometric series

P ex = ∞ k=0 (1 -p 2 )(p 2 P ex ) k = 1 -p 2 1 -p 2 P ex . ( S22
)
Multiplying both sides by 1 -p 2 P ex then yields the quadratic equation (S17) previously found for the Birth-Death process. This is necessary, since we are calculating the same quantity P ex for the same network. It highlights that we may construct a variety of branching processes to find P ex .

Microscopic expressions for p2

We can construct p 2 from microscopic details. In terms of transition probabilities, we find p 2 by summing over all trajectories in (S12):

p 2 = ∞ k=0 Π + 1 (Π - 2 Π + 1 ) k Π + 2 = Π + 1 Π + 2 1 -Π - 2 Π + 1 (S23)
A more detailed description is possible when our description is Markovian (i.e. reactions are sufficiently elementary). Let w + k denote a forward transition rate, to go from X k to X k+1 . Let w - k denote a backward transition rate from X k+1 to X k and w ∅ k the degradation rate for X k . We may then write

Π + 1 = w + 1 w + 1 + w ∅ 1 , Π - 2 = w - 1 w - 1 + w ∅ 2 + w + 2 , (S24) Π + 2 = w + 2 w - 1 + w ∅ 2 + w + 2 .
(S25)

The irreversible reaction limit

Simple Type-I cycles have been studied in the limit where all reactions proceed irreversibly [17,19,40] (∀k w - k → 0). In this limit backward reactions are ignored (Π - 2 = 0), which for our example leads to

p 2 = Π + 1 Π + 2 = w + 1 w + 1 + w ∅ 1 w + 2 w + 2 + w ∅ 2 . ( S26
)
In this limit, there is only one trajectory that contributes to p c , and each step involves a competition between the forward reaction and degradation only.

In studies using simple Type-I cycles, the fraction

ζ k = w + k w + k + w ∅ k , ( S27)
has been referred to as the specificity of reaction step [19, 39? , 40] k, which for irreversible reactions coincides with the transition probability Π

+ k lim w - k-1 →0 Π + k = lim w - k-1 →0 w + k w - k-1 + w + k + w ∅ 1 = ζ k (S28)
For simple Type-I networks with n reaction steps (n distinct edges), the irreversible limit ((S26)) generalizes to [39? ]

p 2 = n k=1 Π + k = n k=1 w + k w + k + w ∅ k , ( S29
)
which we will show more formally in the next section.

The simple Type I cycle with n steps

To expand our discussion on simple Type I cycles, we will now derive a general expression for p 2 , when there are n steps, of which the first n -1 are treated as reversible. The problem is illustrated in Fig. 11c.

Let us denote P X k →X j the probability to reach X j , starting from X k . For P X 1 →X 2 , there is only one trajectory:

X 1 Π + 1 --→ X 2 , ( S30
)
and hence P X 1 →X 2 = Π + 1 . For P X 2 →X 3 , we need to consider that we can go back and forth between X 2 and X 1 :

X 1 Π + 1 ---- Π - 2 X 2 Π + 2 --→ X 3 . (S31)
We can absorb the contribution of hopping back-and-forth in the factor Γ 2 :

P X 2 →X 3 = Γ 2 Π + 2 (S32) Γ 2 = ∞ k=0 (Π - 2 Π + 1 ) k = 1 1 -Π - 2 Π + 1 (S33)
Now, let us consider this argument when we want to find

P X 3 →X 4 for X 1 Π + 1 ---- Π - 2 X 2 Π + 2 ---- Π - 3 X 3 Π + 3 --→ X 4 . (S34)
Now, we need to consider that we can go back and forth between X 3 and X 2 , and that at X 2 we can go back and forth between X 2 and X 1 (captured by Γ 2 ) before moving to X 3 . We absorb all this in the factor Γ 3 :

P X 3 →X 4 = Γ 3 Π + 3 (S35) Γ 3 = ∞ k=0 (Π - 3 Γ 2 Π + 2 ) k = 1 1-Π - 3 Γ2Π + 2 (S36)
We can repeat this argument, e.g. for the kth step

X 1 Π + 1 ---- Π - 2 X 2 Π + 2 ---- Π - 3 ... Π + k-1 -------- Π - k X k Π + k --→ X k+1 , ( S37
)
which then yields the following recursive expressions for k ≥ 1

P X k →X k+1 = Γ k Π + k (S38) Γ k+1 = ∞ s=0 (Π - k+1 Γ k Π + k ) s = 1 1-Π - k+1 Γ k Π + k . (S39)
where Γ 2 imposes that Γ 1 = 1.

The total probability p 2 to reach X n from X 1 and then perform the final irreversible fragmentation reaction r n , is then

p 2 = n-1 k=1 P X k →X k+1 Π + n Γ n = n k=1 Π + k Γ k , (S40)
When backward transitions become negligible (∀k Π - k → 0), we have ∀k ≥ 1 Γ k = 1, and p 2 then acquires its well-known limit expression described by (S29).

A Branching Process for a Type II cycle

As an example of a simple Type II cycle, we consider the network in Fig. 11e, starting with the species C 1 . Our approach will be reminiscent of our branching process in Sec. E 7, but with a repeated branching step.

With a probability p C , C 1 will successfully perform the allocatalytic cycle r 1 + r 2 + r 3 (with some possible back-and-forths), yielding overall

C 1 p C --→ C 1 + D 1 .
(S41)

The probability P D k that k successful cycles occur before the first failure (i.e. degradation C k → ∅) is

P D k = (1 -p C )p k C . (S42)
Corresponding effectively to the overal reaction

C 1 P D k --→ kD 1 . (S43)
Now, let p D be the probability that D 1 succesfully completes a cycle r 4 + r 2 + r 3 (including back-and-forths):

D 1 p D --→ C 1 + D 1 . (S44)
The probability of k successful cycles before failure becomes

P C k = (1 -p D )p k D . (S45)
Combined, a single C 1 has been then replaced according to We find 3 cores of Type I (2 equivalent, up to the choice of reaction to link 5 and 9), and 4 cores of Type II. Species that would need to be provided externally in order to isolate cores are indicated in blue and connected to reactions with dotted lines. k denotes the probability that, starting at B k , the next transition will be a step forward in the cycle, Π - k a step backward, and Π ∅ k a degradation. Fragmentation (yielding 2B 1 ) and degradation (∅) are absorbing states, attained with probabilities pc and 1 -pc respectively. c) Autocatalytic core for a Type I cycle with n nodes. d) Transition network for c e) schematic for an autocatalytic core for a Type III cycle. f) Transition network for e) g) General schematic for an autocatalytic core for a Type II cycle with a single fragmentation step. h) Transition network for g.

Π 2 Π n-1 Π n B 2 B 1 B n B n-1 C 2 C 3 D 1 C 1 1 2 3 4 C 2 C 1 - - Π 2 C 1 + D 1 C 3 D 1 Θ 1 Θ 1 - f ) e) a) b) c) d) g) h) X 1 X 2 X s X n X 2 X 1 + X s X 1 - - ... - X s - X n Π 2 Π n Π s-1 - Π s Π s ...
A 2 A 1 + B 1 A 1 - - A 3 Π 2 e) f ) c) a) d) b) Π 3 Π 3 A 1 A 1 + B 0 A 0 - - A 2 Π 2 Π 0 Π 0 g) j) k) h) i) A 2 B 1 + C 1 A 1 -

Figure 1 .

 1 Figure 1. Different representations for allocatalysis (a,b,c) and autocatalysis (d,e,f). a) Combining reactions (1')+(2') affords an allocatalytic cycle that converts A to B. b) stoichiometric matrix of a), the dashed square encloses the allocatalytic submatrix ν ν ν for network b). c) Graph representation of the allocatalytic subnetwork. d) Combining (1")+(2") affords an autocatalytic cycle converting A to B. e) stoichiometric matrix of d), the dashed square encloses the autocatalytic submatrix ν ν ν for network e). f) a graph representation of the autocatalytic subnetwork.

  Fig 1a and 1d represent examples of allocatalysis and autocatalysis, respectively, with their respective submatrices ν ν ν and ν ν ν , and hypergraph representations Fig 1c and Fig 1f.

Figure 2

 2 Figure 2. a) Five minimal motifs. Orange squares indicate where further nodes and reactions may be added, provided this preserves the motif type (I,II,III,IV,V) and minimality. b+c) Examples of chemical networks, along with their autocatalytic cores. blue: external species, yellow: autocatalysts. b) Type I: Breslow's 1959 mechanism for the formose reaction[33] c) Type II: Another autocatalytic cycle in the formose reaction. Species denoted as Cx inside the nodes refer to molecules containing x carbon atoms, which are shown below in standard chemical representation.

Figure 3 .

 3 Figure 3. a) Pex as function of ζ (legend: Pex(ζ) for Pex < 1) for b) 5 autocatalytic networks of similar size, starting at the dashed node. N1: Type I cycle. N2: Type II with one fork. N3: Type II, two nonoverlapping allocatalytic cycles, a common motif in GARD with a 1st order RAF representation. N4: Type II: allocatalytic cycles connected by intermediate steps. N5: Type V. Symbols: Pex after 1000 simulated trials, detailed in SI appendix Sec. 8, lines: exact solution, derived in SI Appendix Sec. 6

Figure 4 .√ 13 -3 2 .

 4132 Figure 4. Multicompartment autocatalysis. a) Reaction network with two reactions and five species in a single compartment. The network does not contain any autocatalytic core, thus cannot perform autocatalysis. b) Same reaction network as in (a), but duplicated in two compartments α and β, coupled by the selective exchange of species A and A 2 B. A Type II core, highlighted in orange, emerges. c) Open reactor with two compartments, a semi-permeable membrane, degradation and exchange reactions. Chemostatted species have a lighter background: U and V in α and AB in β. d) Extinction probability Pex for multicompartment autocatalysis in (c), starting from a single Aα, as a function of exchange rate k ex and degradation rate k d , relative to other relevant reaction rates fixed at k. Slanting asymptote: exchange-limited survival kex = 2k d . Vertical asymptote: reaction-limited survival k d /k = √ 13-3 2 . Dashed white line: transition between extinction and potential fixation (Pex < 1). Expressions for Pex and asymptotes are derived in SI Appendix Sec. 7.

Definition 8 .Definition 9 .

 89 A reaction is simple if has a single reactant and a single product. Consider a simple reaction R with reactant x and product y, with respective stoichiometries -a and b. The contraction of R consists of: (i) removing R; (ii) merging x and y into a single species z, and; (iii) multiplying by a (resp. b) the stoichiometric coefficients associated with z for all reactions formerly associated with x (resp. y).

Proposition 3 .Proposition 4 .

 34 In a core, a square autonomous submatrix must have a product outside its set of species.Proof. Be A a square autonomous submatrix of M and write M = A C B D . We need to show that B > 0.As A is autonomous, every A-reaction has a reactant in the set of A-species. By Property 2, reactants are always solitary, thus B ≥ 0. Now suppose B = 0. Then det(M ) = det(A).det(D). Either det(A) = 0, implying det(M ) = 0, contradicting M invertible, or det(A) = 0, then A is productive, contradicting the minimality of M . A core is strongly connected.

C 1 -

 1 Figure 5. a) Hypergraph paths and cycle examples. b) Examples of hyper-ears (left) and associated proto-ears (right)obtained by removing green species. The syntax below graphs (see section B 3) describes the relationships between non-P species and path P. Nodes 'u' are products of R, nodes 'v' are reactants of reactions producing x, '-' is any series of reactions with a single reactant and product, '*' denotes the closure of C by the green path. c) Autocatalytic cores. Edges are oriented consistently along cycles, so that reaction have a single reactant. Orange squares are chains of arbitrary length made of reactions with a single reactant species and product species. Edge-to-node connections are weighted by a stoichiometric coefficient, represented explicitly only for Type I by a fork (stoichiometry of 2). C, C and C are cycles. In Type II, the dotted path may comprise multiple cycles similar to the green box. In main text Fig.2, only the case of a single green box is represented for simplicity. In Type I, detC = 0; in Type II, det(C) = 0 and det(C ) can be any value; in Types III-V, det(C) = det(C ) = det(C ) = 0. d) Generic hyper-ear structure of cores. e) Nested and entangled back-branches. f) Example of allocatalytic cycle.

Figure 6 Figure 7

 67 Figure 6. a) A decorated Toy Formose reaction given by the submatrix ν ν ν * , obtained by removing C 1 . The replication cycle g g g D 3 is illustrated in blue. In this network, only species C 2 , C 3 and C 4 are autocatalysts. b) The minimal formose reaction in its autocatalytic subnetwork, an example of an SFA. Arrows illustrate the replication cycle g g g C 2 .

Figure 8

 8 Figure 8. a) Autocatalysts in the Reverse Krebbs cycle, which yield a Type II core. b) Autocatalysts in the Calvin Cycle.We find 3 cores of Type I (2 equivalent, up to the choice of reaction to link 5 and 9), and 4 cores of Type II. Species that would need to be provided externally in order to isolate cores are indicated in blue and connected to reactions with dotted lines.

Figure 9

 9 Figure 9. a) Chemical network for the cavitand-amplification network. blue: feedstock compounds and waste. beige: autocatalysts. b) Type III autocatalytic core for the cavitandamplification network.

Figure 11 .

 11 Figure 11. a) A Type I subnetwork. B 1 and B 2 reversibly interconvert, B 2 can also irreversibly form two B 1 , marked by the forked edge. b) transition network, Π +k denotes the probability that, starting at B k , the next transition will be a step forward in the cycle, Π - k a step backward, and Π ∅ k a degradation. Fragmentation (yielding 2B 1 ) and degradation (∅) are absorbing states, attained with probabilities pc and 1 -pc respectively. c) Autocatalytic core for a Type I cycle with n nodes. d) Transition network for c e) schematic for an autocatalytic core for a Type III cycle. f) Transition network for e) g) General schematic for an autocatalytic core for a Type II cycle with a single fragmentation step. h) Transition network for g.

Figure 12 .Figure 13

 1213 Figure12. a) three motifs with the same Pex(A k ), when the transition network follows the same symmetry as the network structure. b) A more general case: a multiple of allocatalytic cycles of size n. c) 6-membered Type I cycle. d) a 6-membered Type II cycle. e) a 6-membered Type III cycle. f) Transition network for one of the two allocatalytic cycles. g) a 6-membered Type III cycle. The transition of a precursor to allocatalyst is not mediated by an allocatalytic cycle, and hence this is not a RAF. h) Transition network for g). i) a trio of symmetric analogues. j) a Type V autocatalytic core. k) Transition network for j.

  Step k: Suppose D k exists, k < n -1. Re-arrange M such the top left block A of size k corresponds to species set D k and reaction set R(D k ). As A is autonomous, by Property 3, there exists a species x k+1 outside of D k which is downstream D k , hence D k+1 exists.

	Proposition 5. In a core, every species is involved in a
	cycle.
	Proof. Obvious from Property 4, considering the back
	and forth paths joining any two species.

Proposition 6. Consider a partition of a core into two species sets V and W . V cannot be upstream of W , i.e. reactions with products in W cannot have all their reactants in V .

Proof. Suppose on the contrary that M can be written

M = A B 0 C

where A and B span V , C spans W , and B ≤ 0. Given the latter inequality, Property 1 implies that A is non-empty and autonomous. Consider γ > 0 such that M.γ 0. Be α and β the respective restrictions of γ to the reaction spaces of A and B. Then A.α + B.β 0. As B.β ≤ 0, we have A.α 0, contradicting the minimality of M .

ACKNOWLEDGMENTS

AB acknowledges stimulating discussions with D. van de Weem and Z. Zeravcic. AB and DL acknowledge support from Agence Nationale de la Recherche (ANR-10-IDEX-0001-02) and PSL IRIS OCAV). PN acknowledges S. Krishna and C. Flamm for discussions. AB, DL and PN acknowledge Andrew Griffiths for careful reading of the manuscript. AB acknowledges the donation of Jonathan Rothberg for research on artificial life.

Constant composition, constant transition rates

The effect of rare autocatalysts on the steady-state composition (maintained by influx and degradation in a CSTR) is initially small: every reaction introduces changes in molar fraction of the order 1/N (or in concentration terms, 1/V ∝ 1/N ). For large N , the alterations of the composition will be vanishingly small while the autocatalysts are rare.

Consequently, we can approximate the reactor composition in which an autocatalyst is placed, as the steady-state composition before perturbation. We then assume that the molar fractions of species Y k consumed by autocatalysts are sufficiently abundant i.e. x Y k 1/N , which was also required for (S4). For sufficient N , deviations from this approximation become vanishingly small.

A given X k will therefore have a fixed transition rate w + k = kx Y j to perform (S4). Similarly, for (S3) and CSTR degradation, a fixed transition rate w = k is found, depending only on a rate constant.

Independence

In the rare-autocatalyst regime, all reaction steps we consider for autocatalysts are first-order, and they occur at fixed rates. It follows that autocatalysts do not influence each other, and they can each be treated independently. Thus, we can treat each autocatalyst type separately: P ex ({N X k }) = P ex (N X 1 )P ex (N X 2 )..P ex (N X s ). (S6)

Also each individual autocatalyst can be treated as such:

Where P ex (X k ) denotes the extinction probability of a population initially composed of a single X k species. We thus need a method for finding P ex (X k ). This will be obtained by mapping the problem to a branching process.

A branching process

An attractive method we propose for finding P ex (X k ) (from here on simplified to P ex ), is by considering the distribution of 'descendants' X k a species X k will generate. From this, we construct a Branching Process, represented chemically as a single parent molecule X s yielding k descendants:

with p k a distribution of the number of descendants which depends on the network topology.

Knowing p k suffices to find P ex , since the probability to go extinct is the probability that all descendants independently ( (S7)) go extinct:

We will now highlight some possible choices for Branching Processes and their associated p k .

A Birth-Death Process for the Type I cycle

Consider a simple Type-I cycle such as in Fig 11a . Here, simple refers to there being a direct path of (effective) unimolecular steps between the starting compound (B 1 ) and final compound (B 2 ), followed by a single fragmentation step producing two B

Starting from the marked node (B 1 ), let p 2 be the probability of successfully forming 2B 1 , i.e.

where p 2 contains the contribution of all possible trajectories (here: going back and forth between B 1 and B 2 ) that precede the irreversible fragmentation step, i.e.

Where Π + 1 , Π - 2 and Π + 2 are success probabilities for the single reaction steps. These transitions compete with irreversible degradation processes

Where Π ∅ 1 , Π ∅ 2 are success probabilities for the degradation reaction. Due to total probability conservation, we have

Ultimately, a B 1 species will either be replaced by 2 new ones (2B 1 ), or none (∅):

which is a chemical representation of the simplest type of branching process: a birth-death process. In (S9), the only nonzero terms will come from 0 descendants (p 0 ) and 2 descendants (p 2 ):

where

From (S9), we then find P ex . By rewriting the geometric series due to (S49), we have

which admits the solutions P ex = 1 and

12. Microscopic expressions for p C and p D Let us first consider how p C can be constructed from smaller reaction steps. To do so, we observe that the first step must be

with Π + 1 the probability of success, competing with degradation

for which Π ∅ 1 = 1 -Π + 1 . Arrived at C 2 , going back-and forth reversibly becomes possible for neighboring nodes C 1 , C 3 and D 1 :

Where Π + k , Π - k , Θ - k , Θ + k all denote success probabilities. Starting at C 1 , a successful trajectory necessarily starts with reaction r 1 (Π + 1 ) and ends with r 2 + r 3 (Π + 2 Π + 3 ), all in the forward sense. In between, we can go back-and-

) any number of times. Summing over all possible trajectories, p C then becomes

which sums to

Starting at D 1 , a successful trajectory necessarily starts with r 4 (Θ + 1 ) in the forward sense, to form C 2 . From there on, a successful trajectory follows the previous calculation, i.e. p D = Θ + 1 /Π + 1 p C :

In the limit where all reactions proceed irreversibly forward (Sec. E 9), p C and p D only have a contributing from a single straight trajectory

Type III cycles with one fragmentation step

Let us now consider a Type III network composed of n species {X 1 , .., X n } and reaction steps {r 1 , .., r n }, where the final fragmentation step r n produces

as shown in Fig. 11e We may then introduce the success rates for the allocatalytic cycles for X 1 and X s

Which was exactly our starting point in Sec. E 11. Starting at X s , we may thus directly use P ex , upon replacing p C with p s and p D with p 1 , thus finding

We now turn to the problem of finding expressions for p c,s and p c,1 . By structural analogy to simple Type I cycles (apart from the fragmentation step), we may again write

Symmetric motifs

When the network motif is symmetric in structure, and if the transitions preserve this symmetry, this can be exploited to simplify calculations and gain insight in topological aspects of autocatalysis and robustness. Experimentally, this symmetry rarely applies for the transitions, but it can be made applicable for the purpose of our analysis, e.g. by setting transition probabilities to values that reflect the structural symmetry.

Consider a series of m linked allocatalytic cycles (Fig. 12b), all consisting of n nodes and n edges which are structurally equivalent:

which loops back at the mth cycle:

As before, Π + k denotes a forward transition probability, and we also introduce reverse reactions and degradation in the usual sense:

Now, let us choose transition probabilities such that they are equivalent among the equistructural allocatalytic cycles, i.e. periodic in n: Π + k = Π + k+n and idem for reverse steps (Π - k = Π - k+n ) and by extension, degradation

which is seen readily in Fig. 12a by rotating the networks and interchanging the labels. Applying this symmetry to Fig. 12a for the six-membered cycle (n = 2, m = 3), we can thus write P ex (X 1 ) = P ex (X 3 ) = P ex (X 5 ).

(S71) P ex (X 2 ) = P ex (X 4 ) = P ex (X 6 ).

We characterize the success of each equivalent allocatalytic cycle (n steps) by the same probability p 2 . We can then express the extinction probability in terms of the reaction products of one allocatalytic cycle:

where we have used the symmetry in (S72), which yields the exact second order equation we obtained for a simple Type I cycle of n steps. We can thus resolve the general problem using our previously derived solutions.

Appendix F: Expressions for Fig 3, a survey of Pex for various structures

In Fig. 3 in the main text, the behavior of P ex is compared for a number of networks (N 1 to N 5 ), in the limit where all reaction steps proceed irreversibly (Sec. E 9), and where all reaction steps do so with a common success probability ζ (also known in the literature as specificity). Of course, this is an abstraction that is hard to realize experimentally, and its purpose is the following: by controlling for kinetics, we can systematically investigate and compare how survival is affected by network topology. In this section, we will derive the functional dependence of P ex on ζ for the structures discussed in the main text.

N1: a 6-membered Type I cycle

In Sec. E 10, we derived the general solution for nmembered Type I cycles in terms of the probability p 2 to reach X n and perform r n . Starting from X 1 , we can now recover the solution for n = 6

For n = 6, p 2 becomes

In the particular cases of irreversible reactions and for a uniform specificity (∀k

, we obtain p 2 = ζ 6 . Upon injection in the solution P ex = (1 -p 2 )/p 2 (for p 2 ≥ 1/2) ((S19)), we then find

N2: a 6-membered asymmetric Type III cycle

In Sec. E 13, the general solution for n-membered Type III cycles with one fragmentation reaction was derived. For a 6-membered cycle with the fragmentation reaction

which corresponds to network N 2 in the main text.

Having X 4 as the starting species, we can express P ex (X 4 ) in terms of (S66)

In the irreversible reaction limit with fixed specificity (∀k

N3: 6-membered Type II network with RAF representation

The network N 3 consists of two nonoverlapping allocatalytic cycles, which produce each other's allocatalyst:

Choosing our transition probabilities Π + k = Θ + k , we can exploit the network symmetry as outlined in Sec. E 14 and P ex (A 1 ) = P ex (B 1 ). Let p c be the probability that either A 1 or B 1 successfully finishes an allocatalytic cycle. We may then write

The success rate p c can be expressed in the familiar way

In the irreversible limit (∀k ≥ 1, Γ k = 1) with fixed specificity ζ, p c = ζ 3 , and thus we find for p c ≥ 1/2

N4: 6-membered symmetric Type III network

The network N 4 consists of two nonoverlapping allocatalytic cycles, which produce a precursor (A 0 , B 0 ) for each other's allocatalyst:

Choosing our transition probabilities to follow the symmetry of the network, we can write Π + k = Θ + k , and P ex (A k ) = P ex (B k ). Finally, we note that B 0 can either i) degrade with probability Θ ∅ 0 = 1 -Θ + 0 , or ii) form B 1 with probability Θ + 0 = Π + 0 , such that

) Denoting p c the probability that A 1 performs a successful allocatalytic cycle (yielding A 1 + B 0 ), we can write the extinction probability as

which upon injecting (S14) yields the following expression for P ex (A 1 )

Solving the quadratic equation (S16), we find

which yields P ex = 1 and

We can write p c in terms of back-and-forths starting at A 1 , terminating with an irreversible fragmentation

) so that in the irreversible limit with fixed specificity (∀k ≥ 1 Π k+ = ζ, Π k-= 0), we obtain

An alternative way of seeing this is that, by symmetry, P ex is the same as that for a 3-membered Type III cycle with one fragmentation step (forming X 2 ), for which we can directly use the solution derived in Sec. E 13.

A trio of symmetric analogues

As derived in Sec. E 14, the interlinked allocatalytic cycles of size n behave, due to symmetry, as an n-membered simple Type-I cycle. Reproducing the solution for the 2-membered cycles (see Sec. E 6) in the irreversible limit with ∀k ≥ 1 Π + k = ζ, Π - k = 0, we thus find

6. N5: a type V core

In N 5 we have the reactions

If our transitions follow the symmetry of the network, we have

Denoting p c the success probability of the allocatalytic cycle, we can write P ex in terms of (S23)

which is the solution found for the 2-membered cycle. Noting that

the irreversible limit with fixed specificity (∀k ≥ 1 Π + k = ζ, Γ k = 1) yields, as in Sec. E 6,

There are two degradation reactions (r 6 and r 7 ), both taking place in the compartment β:

Here r k refers to the reaction label as found in Fig. 4d, α, β to the compartments and A, AB, A 2 B to species as defined in Fig. 4d. We are guaranteed that, starting from A α or AB α , eventually A β will be formed with probability 1. We can thus simplify the problem by only considering reactions r 3 to r 5 and species A β , A 2 B β and A 2 B α , i.e. we can then write for the transition probabilities

We readily find that this effective description obeys the solution for a simple Type-I network with 3 members, with p 2 a success probability for a cycle

Where we have used Γ k for calculating back-and-forth trajectories, as derived in sec. E 10:

with Γ 1 = 1. Let us now give a microscopic kinetic interpretation to the competing processes, by considering sufficiently elementary transitions on the level of a single species for which we can introduce rate constants

present with a fixed molar fraction x AB β ), which plays the role of a feedstock species in compartment β in r 3 proceeding forward.

• w - 5 = k - 4 exchange, from compartment α to β, when r 4 proceeds backward.

• w + 5 = k + 5 release of A 2 B α through locally irreversible reaction r 5 .

We then obtain

Noting that Γ 1 = 1, the product Γ 2 Γ 3 simplifies to

This allows to fully express P ex in terms of 8 microscopic coefficients. Reactions r 1 and r 2 would give 4 more rate constants and two more molar fractions (for chemostatted species). However, their values do not alter P ex .

For the purpose of illustration, we will consider the competition between exchange, degradation, and other transitions. To do so, we choose one rate for exchange k + 4 = k - 4 = k ex and one rate for degradation k d 6 = k d 7 = k d . Furthermore, we let the sequestration-release steps be equally probable in compartment β, and match release in α: w + 3 = w - 4 = w + 5 = k. The transition success probabilities can then be expressed in terms of two ratios

which upon substitution yields

This permits to construct the phase diagram for P ex in Fig. 4D in terms of the variables ∆ and Ξ.

Phase boundaries and limits

Let us first find an expression for the boundary between autocatalysis and deterministic extinction, which occurs when P ex = 1 and (S7) coincide, i.e. when p c = 1/2, which upon substitution of (S12) becomes:

In terms of Ξ and ∆, we have

Which rearranges to a linear dependence in Ξ ∆ 2 Ξ + ∆ 2 + 3∆Ξ + 2∆ -Ξ = 0, (S19) from which we obtain for the phase boundary

(S20)

In the regime where reactions outpace degradation (∆ 1), extinction will be due to rate-limiting exchange (Ξ 1). Taking (S20), dividing by ∆ and letting ∆ → 0, the ratio Ξ/∆ tends to

as also seen in the phase diagram. When exchange is very rapid (Ξ → ∞), it ceases to be rate-limiting, and degradation will only compete with other reactions. This occurs when we let the denominator of (S20) become 0, i.e.

1 -3∆ -∆ 2 = 0, (S22) which has solutions

≈ 0.30 is the only physical solution, as can also be seen in the phase diagram.

Appendix H: Stochastic simulations

We can numerically sample the extinction probability P ex ({N X } 0 ), with {N X } 0 the initial autocatalyst population, by performing stochastic simulations of the kinetics using Gillespie's algorithm which start from the initial composition {N X } 0 , using transition rates as discussed previously. The compositions obtained from the algorithm are used to sample extinction or survival, as detailed below.

When the total autocatalyst population in a run reaches zero, this is sampled as an extinction event. When the population is not extinct at the end of a prescribed number of reactions n r , the total number N ac of autocatalysts is compared with a threshold population N tresh . If N ac > N tresh , the run is sampled as survival (The error in this approximation scales as P N tresh ex ). Otherwise, the composition after n r steps is used as an initial composition to repeat this protocol, again for n r steps (and so forth if again the threshold is not reached). Once either the threshold is reached or extinction occurs, it is sampled accordingly and the run is terminated.

For networks N 1 to N 5 , the symmetry of the problem is used to sample reactions: a random number χ between 0 and 1 is drawn, and if χ < ζ, a forward reaction step is performed for an autocatalyst picked at random. If χ > ζ, a degradation step is performed for an autocatalyst picked at random.

In Fig. 13a a short time interval is shown for 10 simulation runs performed for the multicompartment network using, for completeness, all 5 autocatalysts AB α , A α , A β , A 2 B β , A 2 B α (Since AB α , A α return to A β with probability 1, the exact solution for P ex can also be found with the three-species network as shown in Sec. G). We start with N A β = 1.0 and all other autocatalysts at 0. Rate parameters for r 3 to r 7 follow the conventions outlined above, and are set at k ex = 1.0,k d = 0.3, k + 3 = 100, k = 10.0. r 1 is treated as irreversible with an effective rate constant k + 1 = 10.0 and exchange of A is performed with a rate k ex . From 100000 simulation runs (N tresh = 40) we obtain the numerical estimate of P ex = 0.698 ± 0.004. This is consistent with our exact solution (S7), which yields P ex = 0.6999.

Among the 10 short runs for the multicompartment network in Fig. 13a, 4 runs reach extinction within the simulation time, whereas the autocatalysis in one run (green) has reached an exponential growth regime that is close to deterministic with a total population of 40 autocatalysts.