Onset of natural selection in auto-catalytic heteropolymers
Reduction of information entropy along with ever-increasing complexity are among the key signatures of living matter. Understanding the onset of such behavior in early prebiotic world is essential for solving the problem of origins of life. To elucidate this transition, we study a theoretical model of information-storing heteropolymers capable of template-assisted ligation and subjected to cyclic non-equilibrium driving forces. We discover that this simple physical system undergoes a spontaneous reduction of the information entropy due to the competition of chains for constituent monomers. This natural-selection-like process ultimately results in the survival of a limited subset of polymer sequences. Importantly, the number of surviving sequences remains exponentially large, thus opening up the possibility of further increase in complexity due to Darwinian evolution. We also propose potential experimental implementations of our model using either biopolymers or artificial nano-structures.
The second law of thermodynamics states that the entropy of a closed system increases with time. Life represents a remarkable example of the opposite trend taking place in an open, non-equilibrium system . Indeed, both information and thermodynamic entropies decrease in the course of Darwinian evolution reflecting ever-increasing complexity of living organisms and their communities . Interestingly, both the second law and the concept of entropy were introduced by Clausius in 1850s, roughly at the same time as Darwin developed and published his seminal work. A century later the connection between life and entropy was highlighted in the classical work of Schrödinger titled ”What is life?” . According to him, living systems are characterized by their ability to ”feed on” and store the negative entropy (which he referred to as ”negentropy”) . In the same work, he effectively predicted the existence of information-storing molecules such as DNA. Soon after, Brillouin established  the connection between the thermodynamic negentropy and its information cousin defined by Shannon .
The emergence of life from non-living matter is one of the greatest mysteries of fundamental science In addition, the search for artificial self-replicating nano- and micro-scale systems is an exciting field with potential engineering applications [6, 7, 8, 9]. The central challenge in both of these fields is to come up with a simple physically-realizable self-replicating system obeying the laws of thermodynamics, yet ultimately capable of Darwinian evolution when exposed to non-equilibrium driving forces.
Chemical networks of molecules engaged in the mutual catalysis have long been considered a plausible form of prebiotic world [10, 11, 12, 13]. Furthermore, a set of mutually-catalyzing RNA-based enzymes (ribozymes) is one of the best known examples of experimentally realized autonomous self-replication. This is viewed as a major evidence supporting the RNA-world hypothesis (see e.g. Refs. -). The ribozyme activity requires relatively long polymers made of hundreds of nucleotides with carefully designed sequences whose spontaneous emergence by pure chance is nearly impossible. Thus, to make the first steps towards explanation of the origin of life, one needs to come up with a much simpler system capable of spontaneous reduction of the information entropy, that would ultimately set the stage for Darwinian evolution e.g. towards functional ribozymes and/or autocatalytic metabolic cycles.
A promising candidate for such mechanism is provided by template-assisted ligation. In this process pairs of polymers are brought together via hybridization with a complementary template chain and eventually ligated to form a longer chain (see Fig. 1a-b). Unlike the non-templated reversible step-growth polymerization used in Ref. , this mechanism naturally involves the transmission of an information from a template to the newly-ligated chain, thus opening an exciting possibility of long-term memory and evolvability. An early conceptual model involving template-assisted polymerization was proposed by P. W. Anderson and colleagues [20, 21]. It has also been a subject of several more recent experimental and theoretical studies [24, 25, 22]. In particular, the model by Hordijk et al.  makes a connection between the classical Kauffman model of autocatalytic sets  and polymer systems capable of template-assisted ligation (see also  for further development of that approach). Recently, we theoretically established  that a cyclically-driven system of this type is capable of producing long, mutually-catalyzing chains starting from a primordial soup dominated by monomers. In the current study we focus on statistics of sequences of these chains and discover that the dynamics of the system naturally results in a dramatic reduction of the information entropy in the sequence space.
Here we further develop the model introduced in Ref. . It describes the emergence of heteropolymers out of the “primordial soup” of monomers by virtue of template-assisted ligation. Our system is driven out of equilibrium by cyclic changes in physical conditions such as temperature, salt concentration, pH, etc. (see Fig. 1).
We consider a general case of information-coding heteropolymers composed of types of monomers capable of making mutually complementary pairs. Polymerization occurs during the ”night” phase of each cycle when existing heteropolymers may serve as templates for ligation of pairs of chains to form longer ones. When the end groups of two substrate chains are positioned next to each other by the virtue of hybridization with the template, a new covalent bond connecting these end groups is formed at a certain rate (see Fig. 1b). During the ”day” phase of each cycle all hybridized pairs dissociate, and individual chains are fully dispersed (see Fig. 1a).
One of the key results of our previous work  is the existence of the optimal hybridization overlap length for template-substrate binding. In this work, for the sake of simplicity we assume that a single pair of complementary monomers is sufficient to bind a substrate to the template. This can be interpreted as if each of monomers in the present model is in fact a ”word” composed of smaller elementary letters, such as e.g. RNA or DNA bases. Within this interpretation, the number of such ”composite monomers” can be exponentially large: , where is the number of elementary letters ( in the case of RNA). As in  we ignore the process of spontaneous, non-templated ligation [25, 19].
In our model monomer types are labeled in such a way that type is complementary to type . One of the key concepts in our analysis is that of a ”2-mer” referring to a monomer immediately followed by the monomer and found anywhere within any heteropolymer. Note that, similar to DNA/RNA complementary strands, polymers in our system are assumed to be directional and anti-parallel when hybridized to each other. Therefore, a 2-mer formed from monomers and is complementary to the 2-mer . It can serve as a template catalyzing the ligation of two substrate chains with monomers and located at their appropriate ends. (see Fig. 1b).
Let denote the overall concentration of 2-mers of type , i.e. the total number of consecutive monomers of types and found anywhere within any chain, divided by the volume of the system. We will refer to the matrix formed by all as the 2-mer matrix. Let denote the concentration of all chains ending with a monomer of type , while - the concentration of all chains starting with a monomer of type . When two ends and of such chains meet due to hybridization with a complementary template , they are ligated at a certain rate to form a new 2-mer . We describe this process by a three-body mass-action term . Here is the ligation rate averaged over the duration of the day-night cycle with the understanding that it happens only during the night phase. 2-mers in our system are assumed to spontaneously break up at a rate . Thus we extend the original model by introducing an explicit sequence dependence of ligation () and breakage () rates. Master equations, describing the slow dynamics in our system occurring over multiple day/night cycles, are:
This mass-action description implies that our system stays well below the saturation regime during the night phase. In other words, we assume that template-substrate hybridization probability is determined by the association rate and not by the competition of multiple different substrates for the same binding site on a chain. This is realized when the duration of the night phase of the cycle is shorter than the typical association time for hybridization. Importantly, this regime also ensures that there is no template poisoning, i.e. the probability of two complementary 2-mers binding each other (and thus loosing their catalytic activity) remains low.
One could write a similar set of kinetic equations describing the dynamics of concentrations of “left” and “right” ends of chains, and . Instead, we use the conservation of overall concentrations of monomers of each type to obtain the explicit algebraic expressions for and in terms of the 2-mer matrix:
Here is the overall concentration of monomers of type in the pool, both free and bound. At the start only free monomers are present (), and thus the initial conditions are given by .
Our model allows for an alternative interpretation that does not involve breakage of intra-polymer bonds. Specifically, as shown in the Supplementary Materials (SM), Eqs. (1-2) also describe a system subject to uniform dilution at rate and the influx of fresh monomers . In the light of this interpretation, below we focus on the case of all equal to each other. Without loss of generality, they can all be set to unity: . This defines the fundamental timescale in our system as either the average lifetime of a single bond or the inverse of the dilution rate.
Spontaneous entropy reduction
In our previous study  we worked within random sequence approximation. If all monomers have identical total concentrations , this approximation corresponds to all matrix elements being equal to each other. For general initial conditions these elements would be proportional to . The key hypothesis proposed but not tested in Ref.  is that the system dynamics would eventually favor the survival of a subset of the ”fittest” sequences at the expense of the others, thus breaking the random sequence approximation. Here we test this hypothesis by simulating the dynamics of the model given by Eqs. (1-2) with . We start with a system characterized by a weak variation in individual ligation rates and concentrations . We choose them from a log-normal distribution with their logarithms having standard deviation 0.1 and means and respectively. Our choice of parameters is motivated by the need to understand the the limit of infinitesimally weak variation of rates and concentrations. For this combination of parameters, the Eqs. 1 are initially linearly unstable with respect to formation of all 2-mers. However, no 2-mer would be formed until either it or its complementary partner is present in the system at least in some infinitesimal ”seed” concentration. Once such seed is introduced, the corresponding pair of mutually complementary 2-mers and would be exponentially amplified. In our simulations we used the same small seed concentration of for all 2-mers
The key parameter we use to quantify the emergent complexity in our system is the information entropy of 2-mers based on their relative concentrations and defined in the standard Boltzmann-Shannon manner:
Fig. 2 shows the time dependence of this entropy in 5 different realizations of and . The entropy starts at its maximal value , and after a brief dip followed by a rebound, it steadily declines as a function of time. Such behavior is a remarkable manifestation of the non-equilibrium nature of our system, as the entropy changes in the direction opposite to that dictated by the second law of thermodynamics. To reveal the source of this entropy dynamics, in Figs. 2b and 2c we show the 2-mer matrix at two time points during our simulations. At all of the 2-mers have grown up from their seed to substantial concentrations. Remarkably, the subsequent dynamics leads to a complete extinction of the majority of 2-mers ultimately giving rise to the 2-mer matrix at shown in Fig. 2c. The time dependence of the logarithm of the number of surviving 2-mers is shown as red lines in Fig. 2a. The ultimate number of survivors is just below (out of ) represented by the lower horizontal dotted line at in Fig. 2a.
Competition between 2-mers and the number of survivors
The observed behavior can be understood from the analysis of Eqs. (1-2). For a fixed set of concentrations and , Eqs. 1 form a set of linear kinetic equations with respect to 2-mer concentrations . Furthermore, this set of equations breaks into independent blocks of equations describing the dynamics of mutually complementary 2-mers and . For a small subset of self-complementary 2-mers , occupying a diagonal of the 2-mer matrix, such a block is represented by a single equation. In all other cases it involves a pair of equations for and coupled via a 2x2 matrix:
Because the trace of the matrix is always negative, at least one of the eigenvalues has to have a negative real part, while the real part of the other one could be positive, negative, or zero depending on the value of the matrix determinant . A negative value of the determinant corresponds to a positive eigenvalue and hence to the exponential growth of two complementary 2-mer concentrations observed at the initial stage. As growing 2-mers gradually deplete and , increases and may eventually turn positive. In this case both eigenvalues become negative. This triggers the exponential decay of concentrations and ultimate extinction of the corresponding pair of 2-mers. A small subset of 2-mers survive and reach the steady state. For these survivors, the determinant has to become exactly zero: . These conditions for surviving 2-mers can be rewritten as
while for all extinct 2-mers .
Now we can put the upper bound on the number of surviving 2-mers in the steady state of the system. This is accomplished by comparing the total number of constraints given by Eqs. 5 to the number of independent variables. Since matrix is fixed, the only variable parameters in Eqs. 5 are the left and right end concentrations and . While naively, the number of such variables is , Eqs. 5 always contain them in combinations . Therefore, for the purpose of our counting argument, only these products should be considered as independent variables. The number of constrains (Eqs. 5) that are simultaneously satisfied cannot be greater than . Each of these equations corresponds to either a pair of mutually complementary 2-mers or to a single self-complementary 2-mer. Denoting the total number of surviving 2-mers as , and the number of self-complementary surviving 2-mers as , the number of equations for surviving 2-mers is given by . Thus the upper bound on the number of surviving 2-mers is given by
Note that for large the number of surviving 2-mers is dramatically lower than - the total number of possible ones. This explains the entropy reduction observed numerically (see Fig. 1). The parameters of the system were chosen in such a way that initially all 2-mers grow exponentially. Since the rate of this early exponential growth depends on and , it differs from one 2-mer to another. This results in a transient behavior where the inhomogeneity of 2-mer concentrations is amplified giving rise to an early decrease in entropy (see the dip around in Fig. 2a). As concentrations and start to get gradually depleted, the growth saturates, giving time for slower-growing 2-mers to catch up with the faster-growing ones around (see Fig. 2b). As a consequence, the entropy recovers close to its maximal value. After that, a new process starts in which 2-mers actively compete with each other for remaining left and right ends. When the determinant for a particular pair of 2-mers and changes its sign to positive, that pair of 2-mers starts to exponentially decay and eventually goes extinct. This process continues until the number of remaining 2-mers (for which ), falls below the upper bound given by the Eq. 6. These surviving 2-mers are visible in the heatmap in Fig. 2C.
A useful visualization of the emergent state of the system is a so-called de Bruijn graph shown in Fig. 3a. It represents each of monomer types as a vertex, and each of surviving 2-mers as a directed edge connecting vertices and . The weight of every edge is proportional to the steady state 2-mer concentration . A de Bruijn graph is a common representation of heteropolymer ensembles such as e.g. DNA sequences of all the chromosomes in an organism. It is straightforward to construct it from a known pool of sequences. However, the inverse problem of reconstruction of the statistics of a sequence pool from the de Bruijn graph is highly non trivial. Each of the polymers in the pool can be represented as a walk on this graph. The simplest case is when the consecutive steps of this walk are uncorrelated with each other. This means that the walk is a random Markov process with the probability of a step given by , while the probability of termination of a polymer at the vertex given by . Note that the entropy defined above and plotted in Fig. 2a is exactly the information entropy of a pool of polymer sequences generated by such Markov process . Longer-range correlations are not captured by the present model, but can in principle could emerge due to the effects outlined in the Discussion section, leading to a further reduction of the information entropy in the system.
The de Bruijn graph can be complemented by another, more compact, graphical representation specific to our system. Since mutually complementary 2-mers appear in pairs and (with the exception of self-complementary 2-mers ), each such pair can be depicted as a single undirected edge connecting vertices to . In this representation each edge represents two 2-mers, while each vertex stands for either or monomer, depending on whether it is the first or the second letter within the 2-mer. For our system this undirected graph (see 2b) has a number of remarkable properties derived in SM. First, it is a so-called ”pseudoforest” : each of its individual connected components contains no more than one cycle. This allows us to refine and give a topological interpretation to Eq. 6: , where is the number of trees (components without cycles) in the pseudoforest. Second, only the odd-length cycles (1,3,5, etc.) are allowed in this graph.
Fig 3c shows the distribution of the number of surviving 2-mers (or equivalently of directed edges in the de Bruijn graph shown in Fig. 3a) in 250 realizations of the system with different values of and . Fig 3d shows the distribution of the number of edges in undirected graphs such as one shown in Fig. 3b for these realizations. As discussed above, the deviation of this last quantity down from is equal to the number of trees in the pseudoforest. As shown in Fig. 3d, it follows approximately exponential distribution with the average around 1.6 (the red line in Fig. 3c). At the same time, the distribution of the number of surviving 2-mers () has a peak around . The quantity is always positive and approximately follows a Poisson distribution with the average of (red line in Fig. 3c).
Variability of the surviving set of 2-mers
The set of surviving 2-mers and their concentrations depend on a number of parameters: ligation rates , total monomer concentrations , and, possibly, seed concentrations of individual 2-mers.
We analyzed the sensitivity of the steady state of our system with respect to all of these parameters one-by-one. First we fixed both , and and analyzed the final 2-mer concentrations for a large number of random realizations of small (but positive) seed concentrations. We found the final state to be completely reproducible as long as all seed concentrations are non-zero. Note that due to autocatalytic nature of 2-mer dynamics, Eq. 1, a pair of complementary 2-mers with zero seed concentrations would never emerge on their own.
Next, we fixed the ligation rates to their values used to construct the heatmaps shown in Fig. 2 and networks in Fig. 3. We then simulated 100 realizations of the system with pulled from a log-normal distribution with and . Fig. 4a shows the heatmap of the fraction of realizations of in which each individual 2-mer survives in the steady state. In Fig. 4b we present the same results in the form of the histogram (blue bars). The majority of 2-mers (324 out of 400 visible as the leftmost bar in Fig. 4b) got extinct in all of 100 realizations. While the number of surviving 2-mers in each realization never exceeded (see Fig. 2c), the overall number of 2-mers that survived in at least one realization of was substantially larger: 76. Out of them, 20 ”universal survivors” were present in all 100 realizations. Furthermore, as can be seen by comparing heatmaps in Figs 2d and Fig. 4a, these 20 2-mers typically have high steady state concentrations . To further investigate the correlation between 2-mer concentration and its survivability, we analyzed a larger set of 250 realizations of the system with the same fixed ligation matrix but variable monomer concentrations . The resulting distribution of 2-mer concentrations shown in Fig. 4c is clearly bimodal. The bimodality is also apparent from an example of de Bruijn graph shown in Fig. 3a, where approximately half of all links (thick lines) correspond to more abundant 2-mers, while the other half (thin lines) - to 2-mers present at low concentrations. The high-concentration peak of the distribution in Fig. 4c is dominated by the contribution from 20 universal survivors shown as red line. Note that despite an increased number of realizations, the set and number of universal survivors stayed the same as in Fig. 4a-b.
We further investigated the variability of the set of surviving 2-mers as a function of width of the log-normal distribution of monomer concentrations. As shown in Fig. 4d, the number of universal survivors (red line) systematically decreases with , ultimately reaching 0 for . Consistent with this trend, the bimodality of the concentration distribution in Fig. 4c also disappears for larger valu. Note that all numbers of 2-mers shown in Fig. 4d were normalized by , i.e. the upper bound of the number of surviving 2-mers in each realization. The average number of survivors in a single realization (green line) does not have a notable dependence on . The blue line in Fig. 4d shows the number of 2-mers (normalized by ) that survived at least once among 100 realizations of . Note this curve grows significantly with reaching the value as high as 5. This corresponds to half of all possible 2-mers having a chance to survive at in least one of these realizations.
The major conclusion following from this study is that our model system of mutually catalyzing heteropolymers has a natural tendency towards spontaneous reduction of the information entropy. This represents an effective reversal of the second law of thermodynamics in this class of systems. While ”violations” of the second law are indeed expected in externally driven non-equilibrium systems, the observed ”reversal” has much greater implications. Both living organisms and other self-organized systems such as human culture, economics, and technology, are characterized by an ever increasing complexity, indicating the ongoing reduction in the information entropy.
The thermodynamic entropy of a system of heteropolymers is composed of two distinct parts : (i) the translational and configurational entropy of polymer chains; and (ii) the information entropy associated with sequence statistics. Our current model hints at a hierarchical scenario of entropy reduction in populations of heteropolymers. First, the translational entropy is reduced due to template-based polymerization as studied in our previous work  within Random Sequence Approximation (RSA). Then RSA breaks down at the level of 2-mers due to their competition with each other for a limited resource of monomers. Such symmetry breaking in the sequence space results in a dramatic reduction in the information entropy (ii). At this point, sequences of the entire pool of chains are generated as Markovian random walks on the de Bruijn graph (Fig. 3a). One can imagine a further reduction of the information entropy due to the emergence of correlations between consecutive steps of this walk.
There are multiple physical scenarios outside the scope of our current model that would lead to such longer-range correlations in the sequence space. They include e.g. a dependence of association rates in Eq. 1 on lengths of the three chains involved in the process of template-assisted ligation. Another intriguing scenario is a spontaneous emergence of chains with weak catalytic activity above and beyond their role as templates for ligation. For instance, some sequences may facilitate breakage or conversely promote ligation reactions, either among some specific sequences or universally. Such sequences would provide a missing link between the prebiotic soup considered here and the emergence of the first ribozymes in the RNA world scenario.
A common pattern in functional RNA-based structures such as ribozymes and ribosomes is the presence of hairpins and loops. Note, however, that mass-action term in the Eq. 1 assumes that the template and the two substrates belong to three different chains thus ignoring these higher order structures. A proper model description taking them into account is a natural next step in development of our approach. An important future question to address is whether our system would naturally evolve towards or away from loops and hairpins. One might expect that a hairpin composed from several pairs of mutually complementary 2-mers would have a self-healing property: any broken bond would be quickly repaired since the template and one of the substrates belong to the same chain and hence are more likely to be hybridized.
An interesting feature of the sequence statistics emerging in our model is that the entropy is not reduced to its absolute minimum that would correspond to a unique ”master sequence”. In the de Bruijn representation such master sequence would look e.g. like a single cycle (or several unconnected cycles) in which for every monomer there is unique right neighbor following it on every chain. In this case the walk on the de Bruijn graph would be completely deterministic. In contrast, in our model each monomer typically has two possible right neighbors in de Bruijn graph. In the limit of relatively small variations used to construct the Fig. 3a, stronger links corresponding to ”universal survivors”, produce -independent backbone of the graph akin to the master sequence. Conversely, weaker links allow for infrequent hopping between different parts of the graph leading to deviations from this master sequence. The opposite limit of large variations in monomer concentrations is especially promising from the point of view of further evolution. In that limit, there is no dominant master sequence. This dramatically expands the explored region in the sequence space: now, for every step of the Markovian walk there are on average two comparable probabilities for selecting the next monomer. As a result, the number of possible sequences of length in our model, , remains exponentially large, yet dramatically reduced compared to its random sequence limit, . Furthermore, in the limit of large , different realizations of give rise to significantly different sets of surviving 2-mers. Note that, unlike the ligation rates , monomer concentrations (or equivalently their influxes ) could vary significantly from one spatial location to another (See  for a study of spatial inhomogeneity in the prebiotic context). This allows for an effective exploration of various regions (of size each) of the global sequence space, rather than converging to the same subset of sequences over and over.
Our model describes a simple yet general mechanism for spontaneous entropy reduction in systems capable of template-assisted ligation. There are multiple possible experimental realizations of such systems based on either traditional DNA/RNA biochemistry or artificial micro/nano structures. The most direct implementation of our model would be a system composed of words made of a string of nucleotides bound together by strong (e.g. DNA-type) bonds. They have to be designed to form mutually complementary pairs that are orthogonal to each other, i.e. words from different pairs have no long overlaps. These words would play the role of composite monomers in our model that could be subsequently connected to each other with weaker, breakable (e.g. RNA-type) bonds. As discussed earlier, our model is directly applicable to the scenario in which all bonds are unbreakable, while the whole system is uniformly diluted and fresh (unbound) monomers are supplied at a constant rate. This greatly expands possibilities for its experimental implementation. Similarly, our dynamics can be implemented e.g. using the DNA origami nanoblocks introduced in Ref. . It should be emphasized, that in order to achieve the behavior described by our model, the experiments need to be conducted well below the saturation regime, i. e. the night phase should be shorter than a typical association time for hybridization.
Yet another interpretation of our model does not involve any long polymers at all. In this case 2-mers are represented by physical dimers made of only two monomeric units incapable of forming longer chains. Our model predicts that, even in this simple system, the compositional entropy would drop because of the extinction of most of the dimers leaving only survivors. On the one hand, this further extends possibilities for experimental implementation. For instance, one could construct the DNA-based system described above but limited to no more than 2-word chains. This would greatly reduce the complexity of the screening process.
On the other hand, this dimer interpretation has an intriguing connection to the Kauffman model of autocatalytic chemical reaction networks . In the Kauffman case, some of the molecules types in the pool are capable of catalyzing the synthesis of others from ”raw materials” (abundant small metabolites) ultimately resulting in the emergence of metabolic autocatalytic cycles in large systems. A recent model of such chemical reactions  shows that the system self-organizes to a state finely tuned to the external driving force. This can be interpreted as maximization of the rate of negentropy adsorption from the environment. In our case dimers correspond to mutually catalytic entities while monomers represent raw materials. While in the current implementation of our model, the catalytic cycles could only involve two mutually complementary dimers, it is straightforward to generalize the model to allow the mutual catalysis of any pair of dimers. We expect our findings about the reduction of entropy to be fully transferable to that case. Thus, our model has a potential of bridging the gap between two traditionally competitive visions of the Origin of Life: ”information first” and ”metabolism first”.
This research used resources of the Center for Functional Nanomaterials, which is a U.S. DOE Office of Science User Facility, at Brookhaven National Laboratory under Contract No. DE-SC0012704.
-  J. Lovelock, GAIA â A New Look at Life on Earth (Oxford University Press, New York, 1979).
-  Strictly speaking, due to the exchange of energy and material with the environment, the appropriate measure of the reduction of thermodynamic entropy by living systems is the local accumulation of Gibbs free energy. This was pointed out by Schrödinger in Ref. .
-  E. Schrödinger, What is Life â the Physical Aspect of the Living Cell, (Cambridge University Press, Cambridge, 1944).
-  L. Brillouin, Science and Information theory, (Dover, Mineola, New York, 1956).
-  C.E. Shannon, A Mathematical Theory of Communication, Bell System Technical Journal,27, 379 (1948).
-  T. Wang, R. Sha, R. Dreyfus, M. E. Leunissen, C. Maass, D. J. Pine, P. M. Chaikin, and N. C. Seeman, Nature 478, 225 (2011).
-  Z. Zeravcic, and M. P. Brenner, Proc. Natl. Acad. Sci. USA 111, 1748 (2014).
-  Z. Zeravcic, M. P. Brenner, Spontaneous emergence of catalytic cycles with colloidal spheres. Proceedings of the National Academy of Sciences 114, 4342-4347 (2017).
-  J. M. Dempster, R. Zhang, M. de la Cruz, Self-replication with magnetic dipolar colloids. Physical Review. E 92, 42305 (2015).
-  M. Eigen and P. Schuster, Naturwissenschaften 64, 541 (1977).
-  F. J. Dyson, Journal of Molecular Evolution 18, 344 (1982).
-  S.A. Kauffman, Autocatalytic sets of proteins, Journal of Theoretical Biology 119, 1 (1986).
-  S. Jain and S. Krishna, Autocatalytic sets and the growth of complexity in an evolutionary model, Phys. Rev. Lett. 81, 5684 (1998).
-  M. P. Robertson and G. F. Joyce, The origins of the RNA world, Cold Spring Harbor Perspectives in Biology 4, a003608 (2012).
-  J. A. Doudna, J. W. Szostak, RNA-catalysed synthesis of complementary-strand RNA, Nature 339, 519 (1989).
-  T. A. Lincoln and G. F. Joyce, Self-sustained replication of an RNA enzyme, Science 323, 1229 (2009).
-  W. Gilbert, Origin of life: The RNA world, Nature 319, 618 (1986).
-  L. E. Orgel, Prebiotic chemistry and the origin of the RNA world, Critical Reviews in Biochemistry and Molecular Biology 39, 99 (2004).
-  C. B. Mast, S. Schink, U. Gerland, and D. Braun, Escalation of polymerization in a thermal gradient, Proc. Natl. Acad. Sci. USA 110, 8030 (2013).
-  P.W. Anderson, Suggested model for prebiotic evolution: The use of chaos, Proc. Natl. Acad. Sci. USA 80, 3386 (1983).
-  P. W. Anderson, and D. L. Stein, in Self-Organizing Systems edited by F. E. Yates, A. Garfinkel, D. O. Walter, and G. B. Yates, (Springer US, 1987) pp. 445-457.
-  W. Hordijk, S. A. Kauffman, M. Steel, Required Levels of Catalysis for Emergence of Autocatalytic Sets in Models of Chemical Reaction Systems. International Journal of Molecular Sciences 12, 3085-3101 (2011).
-  H. Fellermann, S. Tanaka, S. Rasmussen, Sequence selection by dynamical symmetry breaking in an autocatalytic binary polymer model. arxiv.org, arXiv:1708.04779v1 (2017).
-  R. Rohatgi, D.P. Bartel, J. W. Szostak, Kinetic and mechanistic analysis of nonenzymatic, template-directed oligoribonucleotide ligation, J. Am. Chem. Soc. 118, 3332 (1996).
-  J. Derr, M. L. Manapat, S. Rajamani, K. Leu, R. Xulvi-Brunet, I. Joseph, M. A. Nowak, and I. A. Chen, Prebiotically plausible mechanisms increase compositional diversity of nucleic acid sequences, Nucleic Acids Research, 40, 4711 (2012).
-  A.V. Tkachenko, S. Maslov, Spontaneous emergence of autocatalytic information-coding polymers. J Chem Phys, 143, 045102 (2015).
-  G. B. Dantzig, Linear Programming and Extensions, (Princeton University Press, Princeton NJ, 1963)
-  D. Andrieux, P. Gaspard, Nonequilibrium generation of information in copolymerization processes. Proc. Natl. Acad. Sci. USA 105, 9516-9521 (2008).
-  S. Walker, M. A. Grover, N. V. Hud, Universal Sequence Replication, Reversible Polymerization and Early Functional Biopolymers: A Model for the Initiation of Prebiotic Sequence Evolution. PloS one 7, (2012).
-  J. M. Horowitz, J. L. England, Spontaneous fine-tuning to environment in many-species chemical reaction networks. Proc. Natl. Acad. Sci. USA 114, 201700617 (2017).
Equivalence between bond breakage and dilution
Here we demonstrate that the set of equations (1-2 ) also describes the system in which intra-polymer bonds are unbreakable, while concentrations of all molecules are diluted at rate and fresh monomers of type are supplied at rate . Indeed, the dilution adds terms to the r.h.s. of equations 1. The dynamics of individual monomer concentrations (both bound and unbound), , are given by equations . After a brief transient regime, all monomer concentration reach a steady state value so that Eqs. 2 become automatically satisfied. In this way the problem with dilution is exactly mapped onto the problem where bond breakages are allowed.
Topological analysis of the undirected graph
Here we discuss the undirected graph representation of the pool of heteropolymers. Since our system always contains pairs of mutually complementary 2-mers and (with exception of self-complementary 2-mers ), one can represent this pair with a single undirected edge connecting to . These edges form an undirected graph shown in Fig. 3b. Note that due to these rules an edge connecting vertices and in this graph represents a pair of 2-mers and shown as two edges in Fig. 3a or two matrix elements in Fig. 2c. For simplicity in Fig. 3b we did not assign weights to these symmetric edges. Each edge of this graph corresponds to exactly one equation in the set of Eqs. 5. Hence, the number of undirected edges is equal to and according to Eq. 6 it cannot exceed . On the other hand, based on network topology, this number can be expressed as . Here is the number of connected components of the graph, while - the number of independent cycles defined as the minimal number of edges one needs to cut to remove all cycles. The inequality (6) means that the number of independent cycles cannot be larger than the number of components. Furthermore, this inequality must be also satisfied for each of the individual connected components of the graph because the number of equations (edges) cannot exceed the number of independent variables (the number of vertices in this components). In other words, each of the components may contain no more than one cycle. Graphs with this property are known as ”pseudoforests” . For unicyclic components (i.e. those that include exactly one cycle), the numbers of edges and vertices are equal to each other, while for each of the tree (cycle-free) components, the difference between these two numbers is equal to . This gives a topological interpretation to the number of surviving 2-mers in terms of the number of tree-like components of the undirected graph, :
which automatically leads to inequality (6).
One can also demonstrate that only the cycles of odd lengths (1,3,5, etc.) are allowed in our system. Indeed, for a hypothetical even-length cycle , one can construct a combination of Eqs. 5 of the following form:
Here All the variables and at the left-hand-side of this equation cancel, making it an invariant that depends only on ligation rates. Therefore this equation cannot be satisfied for a generic matrix , which rules out the existence of even cycles in most of the cases.