Thermodynamic Basis for the Emergence of Genomes during Prebiotic Evolution

The RNA world hypothesis views modern organisms as descendants of RNA molecules. The earliest RNA molecules must have been random sequences, from which the first genomes that coded for polymerase ribozymes emerged. The quasispecies theory by Eigen predicts the existence of an error threshold limiting genomic stability during such transitions, but does not address the spontaneity of changes. Following a recent theoretical approach, we applied the quasispecies theory combined with kinetic/thermodynamic descriptions of RNA replication to analyze the collective behavior of RNA replicators based on known experimental kinetics data. We find that, with increasing fidelity (relative rate of base-extension for Watson-Crick versus mismatched base pairs), replications without enzymes, with ribozymes, and with protein-based polymerases are above, near, and below a critical point, respectively. The prebiotic evolution therefore must have crossed this critical region. Over large regions of the phase diagram, fitness increases with increasing fidelity, biasing random drifts in sequence space toward ‘crystallization.’ This region encloses the experimental nonenzymatic fidelity value, favoring evolutions toward polymerase sequences with ever higher fidelity, despite error rates above the error catastrophe threshold. Our work shows that experimentally characterized kinetics and thermodynamics of RNA replication allow us to determine the physicochemical conditions required for the spontaneous crystallization of biological information. Our findings also suggest that among many potential oligomers capable of templated replication, RNAs may have evolved to form prebiotic genomes due to the value of their nonenzymatic fidelity.


Introduction
All biological organisms are evolutionarily related. The salient characteristics of life (reproduction and selection) must have therefore emerged either gradually or abruptly from inanimate chemical processes some time in the early history of the Earth. Our ever-increasing knowledge on the biochemical and genetic basis of modern life forms should guide the quest to understand this transition, in addition to the chemistry of potential building blocks [1,2] and geochemical considerations [3,4]. The lack of fossil evidence forces us to rely on model building, which can often be tested experimentally in the laboratory [5]. One of the simplest and most promising is the RNA world hypothesis [1,6,7], which proposes RNA molecules as precursors to modern life forms consisting of DNAs as carriers of genomes and proteins as molecular machines. Continued progress in experimental studies has yielded a diverse range of evidences supporting this hypothesis. In particular, plausible synthetic routes to nucleotides [2] and oligomers [8] have been demonstrated. RNA ribozymes capable of catalyzing RNA replications have been designed and synthesized via in vitro selection [9,10]. Extensive studies of RNA folding landscapes further demonstrate the capability of RNAs to function both as carriers of genotypes and phenotypes [11,12].
Conceptual difficulties to this scenario include the need for the existence of sufficiently concentrated and pure building blocks (chirally selected nucleotides for RNAs) and the necessity to explain subsequent evolutions of multi-chemical autocatalytic systems [13]: the incorporations of proteins and nonreplicative metabolic networks. In this context, Nowak and Ohtsuki recently considered a model describing a pre-evolutionary stage with nonreplicative chemical selection [14]. The undeniable strength of the RNA world hypothesis, nevertheless, is that it has the potential to provide an empirically well-tested pathway for the transition from chemistry to biology, irrespective of its factual historical relevance. The relative simplicity of the model should also allow quantitative descriptions that can complement empirical approaches.
Our focus in this paper, in particular, is the transition from the first RNA molecules formed, which must have been pools of nearrandom RNA sequences, to the first genomes coding for RNA ribozymes. Crucial in understanding such an emergence of the first RNA genomes is the error threshold predicted by the quasispecies theory [15][16][17]. At this threshold, the structure of a population of RNA sequences shifts from being dominated by a stable genome ('master sequence') to becoming random pools, or vice versa. This transition can also be described and understood in the context of more general population dynamics models [18,19], for which many exact results have now been obtained based on statistical physics approaches [17,[20][21][22][23][24]. The error catastrophe transition is in the forward direction, and has thus been likened to 'melting' by Eigen [15]. The transition has recently been observed in behaviors of modern RNA viruses exposed to mutagens [25,26]: a moderate artificial increase in mutation rates of viruses can lead to a complete extinction of virus populations. The error threshold is roughly proportional to the inverse of genome length, which also raised the question of how genomes long enough to encode error correction could have evolved under high error rates (Eigen's paradox) [15,27,28]. Notably, Saakian et al. [29] have recently applied analytical treatments of quasispecies theory to consider this question. Higher organisms keep error rates down to levels that are orders of magnitude lower than achievable by polymerases only, using sophisticated error correction mechanisms including mismatch repair complexes. Tannenbaum et al. [30,31] have studied the quasispecies models of organisms posessing mismatch repair genes, finding transitions analogous to the classic error catastrophe transition in repair-deficient mutator frequencies.
The prebiotic evolution in the RNA world is in the opposite direction of the error catastrophe transition, and may thus be referred to as 'crystallization.' In an equilibrium fluid, whether one observes melting or crystallization is determined by the changes in temperature and pressure. Can we find analogous conditions for the emergence of the first genomes? Addressing this question requires connections to thermodynamics of RNA synthesis. Recent developments in the theory of nucleotide strand replication [32][33][34] provide a promising new direction to bridge the gap between the basic chemical thermodynamics of RNA synthesis and molecular evolution. The mean error rate of replication increases as the reaction condition approaches equilibrium, contributing to entropy production [32]. With a combination of this singlemolecule thermodynamics and quasispecies theory, a surprisingly complete analogy to equilibrium fluids was proposed [34], where volume, pressure, and temperature are replaced by replication velocity, thermodynamic force, and inverse fidelity, respectively, with counterparts of condensation, sublimation, critical point, and triple point. Based on the analysis of a model replication kinetics equivalent to the Jukes-Cantor model of DNA evolution [35], it was suggested that the prebiotic evolution of RNA strands may have been biased by a thermodynamic driving force toward increasingly higher fidelity of polymerase ribozymes below a certain threshold [34].
To what extent these theoretical predictions are applicable to the actual prebiotic evolution that occurred in the past must ultimately be judged based on quantitative empirical data from existing and new experiments. Here, we extend our previous work [34] and assess the applicability of this thermodynamic theory of molecular evolution to prebiotic evolution, using experimental data for polymerization kinetics currently available in the literature. Our results based on these empirical data provide a strong support for the main conclusion of the theory, that there is a thermodynamic driving force biasing random sequence evolutions in the absence of genomes toward higher fidelity in a certain regime of parameter spaces. With considerations of the timedependent evolutionary behavior of RNA populations, we furthermore show that it is possible to estimate the time scales that would have been required for a random sequence pool to crystallize a newly discovered master sequence under a given thermodynamic condition. These results also shed new light on Eigen's paradox. Most importantly, our approach enlarges the scope of both the quasispecies theory-based discussions of the stability of genomes and biochemical approaches to RNA replication by introducing the concept of thermodynamic driving forces and constraints in molecular evolution.

RNA replication kinetics and thermodynamics
The thermodynamic theory of molecular evolution [34] combines the kinetics and thermodynamics of RNA replication on a single-molecule level with population-level features. We first consider the molecular level description of RNA synthesis (or elongation): an elementary step of insertion by addition of a nucleotide ( Figure 1A) consumes a nucleoside triphosphate (NTP) and produces a pyrophosphate (PPi). Its driving force F is given by where F 0 is defined such that F~0 at equilibrium (see Methods). We may estimate the equilibrium constant from DG 0~5 :3 kcal=mol of DNA phosphodiester bond formation and DG 0~{ 10:9 kcal=mol of NTP?NMPzPPi (NMP: nucleoside monophosphate), yielding DG 0~{ 5:6 kcal=mol [36,37]. This value likely overestimates the magnitude of F 0 because it ignores the unfavorable entropy change of binding a free NTP monomer, leading to F 0 * v 9.
One may seek the origin of the observed high-fidelity of polymerization reactions [38] in the relative thermal instability of incorrectly formed Watson-Crick base pairs. However, the stability differences between correctly and incorrectly inserted nucleotide pairs are small: an experimental estimate based on melting temperature measurements for the difference in free energy between incorrect and correct pairs yielded DDG~0:33 kcal=mol [39], which we adopted in this work. This value is the average of the relative stabilities of nucleotides G, C, and T (0:25, 0:33, and 0:41 kcal=mol, respectively [39]) with respect to A in a DNA 9mer duplex terminus against the template base T. The precise value depends on the identity of the base pairs at the terminus and at the neighboring position immediately upstream: for DNAs, duplex stabilities including effects of mismatches can be reliably estimated based on nearest neighbor interactions [40]. Longerranged interactions presumably play more important roles for RNAs, which form secondary structures and higher-order folds [41], affecting DDG values. Frier et al. [42] provided values of free energy contributions to the duplex stability from all 16 possible terminal RNA base pairs and mismatches next to 4 distinct base pairs upstream (Table 4 in Ref. [42]). From these data, we calculated DDG~0:82+0:24 kcal=mol, comparable to k B T.
In quantitative descriptions encompassing both the high kinetic selectivity and this marginal stability difference, it is

Author Summary
A leading hypothesis for the origin of life describes a prebiotic world where RNA molecules started carrying genetic information for catalyzing their own replication. This origin of biological information is akin to the crystallization of ice from water, where 'order' emerges from 'disorder.' What does the science of such phase transformations tell us about the emergence of genomes? In this paper, we show that such thermodynamic considerations of RNA synthesis, when combined with kinetics and population dynamics, lead to the conclusion that the 'crystallization' of genomes from its basic elements would have been spontaneous for RNAs, but not necessarily for other potential building blocks of genomes in the prebiotic soup.
important to fully take into account the reversibility of the reactions [32]. We adopt the simplest description of the kinetics of polymerization, specified by 16 forward and reverse rates, k z (aDb) and k { (aDb), respectively, each corresponding to the insertion and its reverse of a nucleotide (a~A,U,G,C) against a template base (b~A,U,G,C; Figure 1A). In reality, these rates do depend on the identity of base pairs immediately upstream [40,41], which may lead to stalling after incorrect incorporations [28]. More importantly, however, these rates also depend on ½NTP and ½PPi. We estimated the forward rates from the available experimental data of primer extension under the farfrom-equilibrium limiting condition [9,28,[43][44][45][46][47][48][49]. The backward rates can then be related to the forward rates via equilibrium stability.
In general, the overall elongation reaction of a single nucleotide goes through a transition state, whose activation energy is differentially affected by the action of polymerases. If one ignores the reverse reaction under the condition of ½PPi?0, the Michaelis-Menten kinetics applies for the primer extension. In the limit of small ½NTP, we then have k z ?k pol =K d , the latter representing the apparent second-order rate constant with the substrate dissociation constant K d and the turnover rate of product formation k pol [50]. Measurements of polymerase-catalyzed reactions show the selectivity reflected in differences in K d for correct and incorrect base pairs to be orders of magnitude larger than equilibrium stability differences [39]. Examples currently found in the literature are shown in Tables 1 and 2, including those for activated nonenzymatic polymerization (DNA replication  Table 1. Reference base incorporation rates k z (aDb) of NTPs (rows) against template bases (columns).
A. Nonenzymatic [28] A T G C without enzymes) determined recently by Chen et al. [28]. Table 1, in particular, shows the dramatic increase in the degree of relative stabilization of the transition states for correct base pairs in modern polymerases. The evolution of polymerases has entailed two aspects: the facilitation of the overall elongation rate and the amplification of the preferential attachment of correct versus incorrect nucleotides. As we show below, this latter aspect of selectivity evolution leads to a phase transition-like behavior, profoundly affecting population dynamics of evolving macromolecules.
To characterize this dual aspect of enzyme-catalyzed polymerization reactions, we adopt a 'reduced' description involving two key characteristics of forward rates: the mean base incorporation rate k and the relative inverse fidelity c (the ratio of incorrect to correct insertion rates). Precise definitions of these quantities in terms of kinetic rates emerge from the mean field theory (see Table 2. Reference base incorporation rates for DNA polymerases. A. Sulfolobus solfataricus P2 DNAP IV (Dpo4) [44] A T G C Methods): where b Ã is the Watson-Crick complementary base of b and the angled brackets denote a harmonic mean over distribution u(b) of template bases: : ð3Þ Figure 2 shows the distribution of these quantities among nine polymerase systems whose polymerization kinetics have been determined experimentally (Tables 1 and 2), in which we observe qualitative trends of the evolutionary changes reflected on the values of k and c: the k values of modern polymerases are *10 7 times larger than the activated nonenzymatic rate, while the nonenzymatic fidelity (c*10 {3 ) implies that the Watson-Crick structure in the absence of enzymes already supports a fairly high level of fidelity. The arrows illustrate the direction of evolutionary changes that must have occurred from the nonenzymatic to protein-based polymerases via the polymerase ribozymes in the RNA world. The nonenzymatic data are for the templated oligomerization of activated nucleotide analogs, the nucleoside 5 0 -phosphorimidazolide, where PPi is replaced by the imidazole group [28]. Zielinski et al. have compared the kinetics of RNA versus DNA elongation of the activated system [51]. They concluded that RNA elongation is more efficient because its A-form helical structure positions the 3 0 -OH group towards the incoming monomer, whereas contributions of wobble-pairing appeared to facilitate mismatches. This study suggests that the nonenzymatic kinetic rates for RNAs may have higher k and c values than for DNAs. We nevertheless expect their order of magnitudes to be similar.

Mean field theory
The kinetic rates and thermodynamic conditions (the value of F ) allow us to extract, using simulations in general (see Methods), the main stationary properties of RNA elongation: the mean elongation velocity n (the average number of nucleotide pairs added per unit time) and error rate m (the average fraction of mismatched nucleotide pairs). They differ from their respective microscopic counterparts, k and c, because of varying contributions of the reverse rates as functions of F . Importantly, exact analytic expressions for the stationary properties can be obtained if the kinetic rates have sufficient symmetry: the set of k z (aDb) for all a is independent of the identity of b ('symmetric template models;' see Methods and Figure 3).
In Ref. [34], an important special case of symmetric template models equivalent to the Jukes-Cantor model of DNA evolution [35], was considered. The Jukes-Cantor model is a two-parameter model, while general symmetric template models have four parameters. However, to quantitatively assess the applicability of the theory based on empirical data of RNA replication kinetics, it is necessary to allow all 16 values of k z (aDb) to be independent empirical parameters. Here, we used a version of the mean field theory that generalizes the analytic results with the following expressions for the elongation velocity n and error rate m (see Methods): X a q(aDb)~1, ð5bÞ where q(aDb) denotes the probability to find a base-paired to b, and Eq. (5b), the normalization condition for q(aDb), determines n(b), the mean velocity of nucleotide addition against template base b. Equation (5a) is a generalization of the equilibrium Boltzmann distribution, to which it reduces to when n(b)~0, and is exact for symmetric template models ( Figure 3). Because the complete reproduction of an RNA strand requires a pair of replications, we also considered the net error rate m 2 of two consecutive replications: For connections to thermodynamics, one must calculate the entropy production (in units of k B ) per monomer addition [32]: where the first term in the square brackets represents the contribution of monomer consumptions to dissipation and the second term corresponds to the disorder creation by copying errors. The average in Eq. (7) is an arithmetic mean since F is not a rate and should match the external thermodynamic force given by Eq.
(1) in stationary states. The quantity {g(aDb) is the free energy change in units of k B T (or the negative of entropy production) for the addition of a nucleotide a against b, with which the backward rates k { are expressed in terms of forward rates k z via The dependence of g(aDb) on concentrations of monomers is nontrivial because four NTPs compete for a single site. Relative stabilities of correct versus mismatched base pairs (DDG), in contrast, are expected to be largely insensitive to concentrations. The following form of g(aDb) reflects this expectation [34]:  (Table 2F). All symbols were calculated from numerical simulations. C-D: Mean velocity (C) and error rate (D) for the pol b kinetics, both with full experimental kinetics ( It can be shown that m ranges from a minimum m min~3 c far from equilibrium (F ??, n?k), leading to Eq. (2b), to a maximum m max~3 =(3zg) at equilibrium (F~n~0) (see Methods). This variation of m with varying thermodynamic force F can be interpreted as follows: near equilibrium, both the correct (faster) and incorrect (slower) incorporation steps are balanced by their reverse steps, leading to comparable net incorporation statistics. Far from equilibrium, the reverse rates become negligible and the faster correct incorporation dominates.
In Figure 3A, we show that the mean field theory is exact for arbitrary symmetric template models. Figure 3B supports the siteindependence of q(aDb) [Eq. (35) in Methods] for more general 16-parameter cases. Comparisons of the mean field theory predictions for elongation properties of pol b kinetics (Table 2) with simulations ( Figure 3C-D) show that the theory generally gives reliable results. The Jukes-Cantor reduction of empirical rates [Eq. (4)] based on Eqs. (2) is also seen to give a good approximation over all parameter ranges, showing that the analytical theory developed in Ref. [34] provides accurate descriptions of realistic kinetics. Nevertheless, for the best numerical accuracy of predictions based on experimental kinetics, we based our main results in the following sections on stochastic simulations. Importantly, however, the mean field theory in the current application yields the definitions given by Eq. (2), in addition to the analytical limits of velocity and error rate (see Methods), which we verified exactly from simulations.

Single-molecule properties
We applied this single-molecule description of RNA replication to three experimental systems: nonenzymatic reactions [28], 'Round-18' (R18) polymerase ribozyme [9], and poliovirus polymerase (3D pol ) [43] (Figure 3), each representing the beginning, intermediate, and late stages of evolution. As has been previously observed in Ref. [34] for the Jukes-Cantor model, the qualitative trend shown in Figure 4 parallels that of fluids undergoing vapor-liquid transitions with decreasing temperature when pressure, volume, and temperature are replaced by thermodynamic force F , velocity n, and inverse fidelity c, respectively. The correspondence of c to temperature in fluids, in particular, is natural because it is a microscopic measure of randomness destroying genomic information. Figure 4A,B shows that for high inverse fidelity c (nonenzymatic), the elongation velocity n and error rate m monotonically increase and decrease, respectively, with increasing F . A critical point is crossed (ribozyme) as c decreases, and n and m become nonmonotonic (3D pol ) with discontinuous jumps for decreasing F ('evaporation'). The error rate m 2 exhibits the same qualitative behavior ( Figure 4B). These results verify the biological applicability of the theoretical predictions made previously in Ref. [34], based on known experimental kinetic data of systems representing key milestones of evolutionary processes ( Figure 2).
The key question then is: how would these changes in the elongation behavior of RNA replication actually have occurred during the prebiotic evolution? To address this question, we modeled the increases in fidelity from the nonenzymatic starting point by uniformly rescaling the incorrect incorporation rates [k z (aDb), a=b Ã ] of the set of nonenzymatic kinetics (Table 1) to produce different values of c. Simulations identified the critical point suggested in Figure 4A,B at c c~7 :6|10 {4 and verified the limits of error rates at and far from equilibrium predicted by the mean field theory exactly ( Figure 4C,D).

Phase behavior
We then scanned the variation of these phase behavior for different values of c and F to generate the phase diagrams shown in Figure 5, which confirms that the qualitative features of the Jukes-Cantor model phase diagram [34] are preserved for empirical RNA replication. However, as opposed to the results in Ref. [34] that represent generic predictions, Figure 5 is based on empirical nonenzymatic kinetics and its uniform rescaling, with no other adjustable parameters. The discontinuous jumps shown in Figure 4C,D correspond to the limit of stability ('spinodal'; thick black lines) of the 'liquid' or L phase (high n-low m 2 state) against the 'gas' or G phase (low n-high m 2 state). In equilibrium, the location of a phase transition in the phase diagram is determined by the equality of free energies of the two phases [52,53]. Here, we adopted the assumption that if multiple stationary states exist for a given F , the state with higher n (and lower m 2 ) is chosen. This assumption is based on the relationship connecting the entropy production rate _ S S to F and the velocity n l of RNA replicator l present in the system, where N is the total number of replicators and n g is the mean velocity (or 'fitness'). The analogy to equilibrium phase behavior also excludes the first-order character of liquid-solid transitions, which for the current case is continuous. Equation (10) is a special case of a general relationship between nonequilibrium fluxes and conjugate forces [52]. In this formulation, a state with high n g contributes more to entropy production. This assumption is consistent with the standard interpretation of the replication rate as a measure of fitness [15,54]. The multiplicity of stationary states at a single-molecule level is supported by the recent demonstration of a real-time sequencing-by-polymerization technique [55], where it was reported that polymerases interconverted between two distinct velocities during DNA elongation for a given reaction condition ( Figures 3C and S3 of Ref. [55]). A complete kinetic characterization of the w-29 polymerase used in this experiment would allow us to make a more quantitative assessment of this interpretation.

Population dynamics
In considering the thermodynamic interpretation of the population dynamics of RNA sequences, we adopt the following physical model (Figure 1): during evolutionary drifts of a random population in sequence space, a particular sequence that folds and catalyzes the replication of RNAs with the same sequence (and no others) is 'discovered.' (In reality, a ribozyme would more likely have had catalytic activities for arbitrary sequences. The selectivity toward its own sequence, instead, would have arisen from the need for spatial diffusion in order to act on other sequences.) This sequence therefore has a higher k value [Eq. (2a)] compared to others, leading to the single-peak Eigen landscape [Eq. (18) below]. Our goal in this and the following subsections is to describe the growth and stability of this master sequence. In Ref. [34], the basic quasispecies theory under the single-peak landscape was combined with the theory of a single-molecule elongation. We expanded this treatment by considering different scenarios of how F and c may have been distributed in RNA populations ( Figure 1B).
For the inverse fidelity c, one may first assume that it is nearly uniform (or regard it as an average over replicators) in a population, as has been assumed implicitly in Ref. [34]. We also assumed that only the RNA strands with a certain polarity (analogous to the positive or negative-sense polarities of viral genomes [56]) have catalytic activities, such that a pair of replication events is necessary to reproduce a polymerase ribozyme. This feature makes the current treatment more realistic for RNA prebiotic evolution compared to those in Ref. [34]. The following derivation of the thermodynamic quasispecies theory in this subsection otherwise adopts the approach therein [34].
In a population of self-replicating RNAs with genotypes labeled by index i, the genotype i catalyzes replications with rates where k i is the equivalent of Eq. (2a) for the genotype i specified by a fitness landscape. The relative rate k(aDb) specifies the rate of addition of nucleotide a against base b, all normalized such that S P a k(aDb)T b~1 . In this model, therefore, all genotypes have the same set of relative enzymatic rates for nucleotide pairs (and the same value of c), while differing in their absolute magnitude of catalysis, k i . This assumption of uniform inverse fidelity is reasonable for populations with genotypes distributed within a small neighborhood of a master sequence (or a small random subspace in the absence of a master) in the sequence space. The elongation velocity n i is given by where (the relative velocity) n is now determined from Eqs. (5c) with k z (aDb) replaced by k z (aDb), and the replication rate of genotype i is because a pair of replication events requires the addition of 2M nucleotides, where M is the length of genome. The mutation rate Q ij from genotype j to i is given by where d ij is the Hamming distance (the number of nucleotides that are different) between i and j. Denoting the number of individuals (of the polarity that has catalytic activity) with genotype i as n i , the evolving population in the Eigen model [15,16] without constraints on the population size obeys the dynamical equation, _ n n i~X j Q ij r j n j : At any time t, the total number of all individuals (population size) N is given by N~P i n i , which from Eq. (15) changes via _ N N~rN, where r~P i r i p i is the population growth rate (mean fitness) with the frequency of genotype i, p i~ni =N. Therefore, for a given population characterized by the set fn i g, the corresponding entropy production rate is given by Eq. (10) with Similarly, under an idealized condition where replication occurs together with degradation [15,17,57], a population can evolve under a constant F with a fixed mean population size. In this case, a replication event occurs with the same rate as the random degradation of a replicator. The evolution equation becomes such that N is constant. For the fitness landscape, we adopted the single-peak Eigen landscape: where k 0 is a constant with the unit of a rate and Aw1 is the relative fitness of the master sequence. Under these simplifying approximations, the standard quasispecies theory becomes applicable directly, with connections to thermodynamics made by m 2~m2 (c,F ) and n~n(c,F ). These fundamental relationships linking elongation properties to thermodynamic and kinetic parameters can be written in implicit but closed analytical forms [34] for the Jukes-Cantor model. In the numerical approach adopted here for arbitrary rates, simulations are first performed for a given set of rate constants and c values to obtain averages of n, m 2 , and F values as functions of g, as illustrated in Figure 1 of Ref. [34]. The implicit parameter g is then eliminated to obtain n and m 2 as functions of F ( Figure 3C-D). For the region in which multiple branches of n exist for a given F , the branch with the largest n (L phase) is chosen.
In the infinite population limit, the quasispecies is either dominated by the 'master species' with the mean fitness r~(k 0 n=2M)A(1{m 2 ) M^( k 0 n=2M)Ae {m 2 M , where (1{m 2 ) M is the probability of replicating M sites over two consecutive cycles without error [see Eq. (6)], or by the 'mutant species' with fitness r~k 0 n=2M. From Eqs. (13) and (16), the mean fitness is therefore given by where m Ã 2 (c,F )~l denotes the threshold error rate for which n g becomes the same for the master mutant species. Equation (19) implies that a constantm 2 contour (red lines in Figure 5) is the 'melting line' separating a 'crystalline' (C) phase from the L phase [34]. As shown in Figure 5, the L-C transition line meets the L-G line at the triple point, below which the C and G phases meet directly ('sublimation'). Our results show that this fairly complete analogy to the equilibrium phase behavior of fluids discovered first in Ref. [34] is indeed equally applicable in more realistic considerations of RNA prebiotic evolution.
The L-C transition line lies at the heart of the crystallization of genomes that may occur during evolutionary walks [13] in sequence space. The presence of the L phase distinct from the G phase below the critical point has an important consequence to such sequence explorations: despite the absence of a stable genome, analogous to liquid phases with short-range orders, RNAs in the L phase with m 2 *10 {2 ( Figure 5D) would still exhibit sequence correlations for a significantly large number of generations. We may use the Jukes-Cantor relationship between the error rate and the cumulative mean Hamming distance d from an ancestral sequence after l generations [35], A typical sequence in the G phase with m 2^0 :5 ( Figure 5D), for instance, would evolve to reach d~M=2 in just l^1:6 generations, on average, whereas in the L phase with m 2^0 :02, it would do so in l^41 generations. Therefore, when a system 'evaporates' into the G phase, an ancestral sequence gets lost in a couple of generations. In contrast, conditions in the L phase, with error rates comparable to those in the C phase nearby in the phase diagram, would greatly facilitate crystallizations of viable genomes. In interpreting the physical distinction between L and G phases, it is useful again to compare them with their analogs in equilibrium fluids, the liquid and gas phases in a container. The pressure of a fluid in equilibrium is controlled by the external force per unit area of the container, which matches the average of microscopic forces per unit area exerted by molecules on the wall interior. At high temperatures (the average kinetic energy of molecules), a given external pressure can be balanced by the mean force of a state (gas), where density is low and molecules rarely interact. The equilibrium density is then roughly proportional to pressure and inversely proportional to temperature. At low enough temperatures, a given external pressure can also be matched by a different phase (liquid) with a much higher density held together by intermolecular attractions. Both gas and liquid phases are characterized by the lack of long-range order. The sharp boundary between them appears when temperature goes below the critical value because the effect of molecular interaction renders a certain range of pressure values unstable.
Analogously, an RNA molecule replicating in a chemical reservoir is driven by the external thermodynamic force given by Eq. (1), which matches the average entropy production per monomer addition. For large c values, the replication is nearly random and the second term of Eq. (7), the sequence disorder contribution to the entropy production, is constant (^ln 4), making the dependence of internal F on g monotonic [34]. With a sufficiently small c, in contrast, the sequence disorder nearly vanishes, reducing the entropy production. This change is compensated by the dominance of faster correct incorporation steps, with the corresponding increase in velocity and decrease in error rates. A given value of external F can be matched either by a state with low velocity and high errors (G phase), or by one with high velocity and low errors (L phase), each distinguished by the relative importance of the two terms in the square brackets in Eq. (7). The sharp boundary between them appears because, for intermediate values of n, stationary states become unstable against fluctuations. The neighborhood of regimes where the C phase is stable is dominated by the L phase ( Figure 5) in which the error rate is comparable to those in the C phase, if c is subcritical.
For the population as a whole, random drifts in c due to sequence explorations are not isotropic but, rather, are biased toward the direction of increasing n g . In our previous work [34], a threshold was identified within the phase diagram separating regimes where the direction of this bias shifts. We sought the analog of this threshold in Figure 5 corresponding to the evolution of RNAs, where the region in which l~Ln g =Lcv0 in the L phase (white dotted lines in Figure 5A, B) includes the nonenzymatic fidelity value and links it to the C phase. Inside the C phase, l is always negative (Figure 6A). Once a population has c values to the left of the white dashed line in Figure 5A, random drifts in sequence space would be biased toward increasingly higher fidelity, leading to crystallization and stable genomes.

Stochastic evolutionary dynamics
We next relaxed the assumption that c is uniform within a population (c?c i is the value for genotype i). Equation (13) is then replaced by for which numerical simulations have to be used. An efficient method to extract collective population dynamics of competing molecules is again provided by the Gillespie algorithm [58], which was first applied to the quasispecies dynamics by Nowak and Schuster [59]. The set of possible reactions corresponding to Eq. (17) a population can undergo are written as where R i is a replicator of genotype i. The mutation matrix is given by where m i~m2 (c i ,F) is the error rate of reactions catalyzed by the genotype i. We tested this simulation algorithm using the special case of the exponential growth of a population with no degradation [Eq.

Crystallization kinetics
We used the constant-N stochastic evolutionary dynamics simulations to examine the temporal evolution of quasispecies. The inverse fidelity c i was assumed to depend on genotype i via the same form of single peak landscape as for fitness: Figure 6. Dependence of mean fitness on fidelity. A: Mean fitness as a function of c at constant F . The slope l~Ln g =Lc is negative below a threshold c for each F (white dashed lines in Figure 5A,B and green dashed lines in Figure 5C,D, respectively). The discontinuous jump for F~0:4 and the cusps at smaller c values correspond to G-L and L-C transitions, respectively. B: Mean fitness averaged over starvation processes (F 0~5 ) for different initial thermodynamic force F i (see Figure 10). The slope l s is negative below a threshold c for each F i (green dotted lines in Figure 5A where we took c 1~1 :0|10 {4 and c 2~4 :23|10 {3 (the nonenzymatic value) in Figure 9. The time scale of simulations is set by using k 0 =½NTP~6:65|10 {9 mM {1 s {1 from the nonenzymatic replication (Table 1) in Eqs. (18) and (22). We assumed ½NTP~1 mM as a representative chemical environment, such that k 0~0 :210yr {1 . Figure 9 shows two typical trajectories starting from an initial pool of random sequences of length M~10 2 , containing a single replicator designated as the master sequence. This 'seeding' of the population by a master sequence mimics the situation where a genotype with a significantly higher fitness is discovered during random drifts. The resulting evolution in Figure 9 is analogous to 'crystal growths,' in which the frequency of the master sequence steadily grows to reach a value consistent with the stability of the C phase ( Figure 5): a master sequence with A~e and c i~1 0 {4 spreads and dominates the population under thermodynamic force F~1:0. The corresponding growth under F~0:5, which corresponds to the vicinity of the L-C boundary in Figure 5, is much weaker and slower, suggesting that the phase diagram remains valid for inhomogeneous c i . The estimated time scales in Figure 9 (based on the activated nonenzymatic rates and ½NTP*1 mM) further suggest that the crystallization of a genome can occur within *10 4 yrs under suitable conditions. However, as in equilibrium fluids, it will never occur if thermodynamics precludes a stable C phase.

Starvation process
We also considered an alternative setup where a population growth occurs in a closed system, which leads to an evolutionary change we refer to as the 'starvation process.' Similar situations were also considered in Ref. [60]. During an idealized starvation process, a single genotype is placed inside a medium containing a given amount of NTPs and PPi's with the corresponding initial thermodynamic force F i . The population growth leads to the gradual depletion of NTPs and accumulation of PPi's, lowering F .  The error rate therefore increases over time. The growth of the population would come to an end when the condition finally reaches equilibrium (F~0). The resulting collection of RNAs in reality may then disperse into fresh media, restarting new rounds of starvation processes.
We introduce the fractional population size j with respect to the asymptotic population reached in the limit of equilibrium (see Methods). A single process starts with F~F i , where n g is maximum and j~0, and may undergo up to two transitions (C-L and L-G) if cvc c to reach equilibrium, where F ?0, j?1, and n g ?0 ( Figure 10). In Figure 6, the mean fitness as a function of c averaged over a starvation process (B) is compared with that without the averaging (A). For c close to the nonenzymatic value (vertical dotted line), l s~L Sn g T=Lc is negative in the L phase for F i above a threshold. The dependence of the l s v0 region on F i (green dashed line in Figure 5A,B) closely resembles that of the lv0 region on F (white dashed lines in Figure 5A,B). We therefore conclude that an environment that supports repeated starvation processes with an initial F above this boundary for a given c promotes evolution that lowers c. It is worthwhile to note that this conclusion was reached without invoking any significant simplifying assumptions other than the experimentally characterized kinetics of nonenzymatic replication (Table 1), thermodynamic considerations, and the quasispecies theory, except the uncertainty in values of F 0 . We verified that the conclusion remains valid for all possible F 0 (Figure 11).

Evolution of longer genomes
The conclusion that there was an underlying driving force biasing fidelity increases in the absence of genomes is particularly powerful because it is independent of the physical mechanisms implementing it. A likely mechanism for such changes is the evolution of error correction with the necessary increases in genome length M. The Eigen's paradox arises because such an increase would lower m Ã 2 (the melting line recedes toward smaller c in Figure 5A,B). The melted population, however, would be driven to recrystallize a new, longer genome because lv0 in the L phase. Growths in genome lengths most likely occurred with insertions, which is beyond the scope of our treatment that only considered base substitution errors. Saakian has studied the evolutionary model of parallel mutation-selection scheme with insertion and deletion [61]. Similar approaches combined with our findings may offer more detailed insights on how genome growths may have been facilitated by thermodynamic driving forces. In addition, we have restricted our study here to a single chemical system (RNAs). It would be of interest to apply similar approaches to more complex systems containing multiple ingredients, including peptides.
Together, our findings in Figure 5 suggest that the initial nonenzymatic fidelity of RNA lies within the threshold favoring fidelity increases. Rather than being coincidental, this feature may explain nature's choice of NTPs as the media for encoding biological information. Many possible alternative oligomers capable of templated replications have been proposed as precursors to RNAs [7]. Their corresponding monomers, however, would have had widely different fidelity values, and one system (NTPs) that happened to lie within the lv0 boundary presumably evolved the RNA quasispecies cloud towards smaller c, eventually crystallizing the first genomes.

Thermodynamic force
An elementary elongation reaction can be written as where a~A,U,G,C and R n is the RNA primer of length n. The entropy of the system plus reservoir is S~S(n,N a ,N p ), where N a and N p are the total numbers of monomers a and PPi, respectively. The entropy production rate is [34,52] where f~TLS=Ln is the force acting on the growing primer (in length units of e.g., base pair rise), m p~{ TLS=LN p , m a~{ TLS=LN a , n~_ n n~_ N N p~P a n a , n a~{ _ N N a is the consumption rate of nucleotide a. The force f is analogous to the external tension balancing the entropic force of rubber elasticity [52]. In conditions where the external force is not controlled, it may be replaced by frictional drag on polymerases, which would depend on elongation velocity. The constant F a is given by In Eq. (28)

Equilibrium condition
A useful physical insight to g defined in Eq. (9) can be gained by considering the condition of equilibrium, which can be derived from Eq. (7) as or with Eq. (8), k z (aDb)~k { (aDb)q eq (aDb), the detailed balance. In equilibrium, on the other hand, we can calculate q eq by considering a two-level system with a ground state and m{1fold degenerate excited states with energy gap DDG (m is the total number of NTP types; m~4 in the main text): Comparing this with Eqs. (9) and (29), we obtain which we also verify in the next subsection directly from the mean field theory. We may interpret m{1 and g in Eq. (31) as the entropic and energetic factors for mismatches: they are costly individually (by a factor of g) but there are many of them (mw1). Equation (29) shows that the free energy parameter g is defined with a constant term such that it absorbs the partition function that normalizes q eq (aDb) in equilibrium. The mean field theory generalizes Eq. (29) into Eqs. (5a) and (8) for nonequilibrium conditions.

Mean field theory of RNA replication
Previously, we derived a mean field theory for the templated replication and showed with numerical tests that it was exact for the two-parameter Jukes-Cantor model rates [34]. Here, we reproduce the analytical derivation and expand it to show that the mean field theory becomes exact for symmetric template models (i.e., four-parameter models with the set of k z (aDb) for all a independent of the identity of b), of which the Jukes-Cantor model is a special case. We also test this conclusion by comparing the analytical results with numerical simulations (Figure 3). The particular version of the mean field theory we adopted for general rates in this paper [Eq. (5)] can then be considered as a generalization of these expressions.
Considering the elongation of RNA depicted in Figure 1A, the master equation for the probability P n (a n Db n ) to have a chain with length n and sequence a n :fa 1 , Á Á Á ,a n g under a template b (assumed infinitely long; may refer to the whole template sequence or a single base) can be written as d dt P n (a n ,tjb)~k z (a n jb n )P n{1 (a n{1 ,tjb) z X a 0 k { (a 0 jb nz1 )P nz1 (a 0 a n ,tjb){k { (a n jb n )P n (a n ,tjb) where a,b~1,2, Á Á Á ,m (m~4 and a,b~fA,U,G,Cg in the main text). We introduce the reduced distribution [32] P (l) n (a n a n{1 Á Á Á a n{lz1 ,tjb)~X fa i g,iƒn{l P n (a n ,tjb) :p n (tjb) q(a n Á Á Á a n{lz1 jb), where p n (tDb) is the probability to find the chain with length n at time t under the given template b, and q is the conditional probability of having the indicated sequence for the given chain of length n. In writing q as time independent, we have implicitly assumed the stationary limit where q represents the asymptotic monomer distributions near the terminus of the growing chain. Our numerical simulations confirm the existence of such distributions for l~1. Conversely, the chain length distribution p n (tDb) supports a peak, SnT~P n np n (tDb), moving with a constant velocity, n~P n ndp n (tDb)=dt, which depends on the entire template sequence b in general. A special case of particular interest is the symmetric template models, for which the set of rates would be specified completely by at most m parameters instead of m 2 . The simplest example is the two-parameter model [34] given in Eq. (4). For such models, we may write The physical interpretation behind Eq. (34) is that the monomer addition and deletion at the n'th site on the template are solely determined by the set of rates k z (a n Db n ), which are independent of the identity of b n . The probability for chain growth, therefore, should be independent of template sequences. We also assume that which is expected because the rates k + (aDb) are all local in their dependence on nucleotides. Numerical tests suggest that Eq. (35) is generally valid for arbitrary rate constants ( Figure 3B). Using Eqs. (34), (35) and summing both sides of Eq. (32) over a i (i~1, Á Á Á ,n{1), we have _ p p n q(a n jb n )~k z (a n jb n )p n{1 z J { p nz1 {J z p n {k { (a n jb n ) p n ½ q(a n jb n ), where J z~X a k z (aDb), ð37aÞ are the total fluxes for chain growth and shrinkage. Any expression involving summations over a of k z (aDb) is independent of b for symmetric template models. Equation (36) is valid for any fa n ,b n g, which we replace by fa,bg. We note that P n p n~1 , P n _ p p n~0 , and P n n _ p p n~n . If we sum both sides of Eq. (36) over a using P a q(aDb)~1, multiply by n, and sum over n, we get If we sum both sides of Eq. (36) with respect to n, or with Eq. (38), Equations (37b), (38), (40), and P a q(aDb)~1 which was used in the derivation, form a set of self-consistent equations for q(aDb) [mz2 equations, m copies of Eq. (40), the normalization, and Eq. (37b); for mz2 unknown, q(aDb), n, and J { ]. Because n is independent of b, this set of equations is most conveniently solved by imposing the normalization condition to Eq. (40) to determine n: We note that Eq. (41) leads to a unique n independent of b because of the symmetry of rates with respect to b. The mean error rate is given by again independent of b. In Figure 3A, we test Eqs. (40) and (41) with an example set of symmetric template model rates. The comparison with numerical simulations shows that the expressions are exact for symmetric template models. We also tested the factorization assumption, Eq. (35), for the general case (pol b kinetics, Table 2) in Figure 3B, which suggests that it is generally valid.
Without the symmetry of kinetic rates with respect to template base identity, the main difficulty in the analytical treatment is that Eq. (34) is no longer valid, and the velocity n depends on the entire template sequence. There are a number of ways to generalize the exact expressions, Eqs. (40) and (41), into cases where the kinetic rates do depend on the identity of template bases. One way, demonstrated in Ref. [34], is to introduce averages over template bases to Eq. (36) to symmetrize q(aDb) and n over the distribution of b. An average over b of the right hand side of Eq. (41) determines the mean velocity n independent of b (Eq. (12) of Ref. [34]). Here, we used a different approach, generalizing Eq. (40) into Eq. (5a), which introduces a template-dependent velocity n(b). We found this mean field theory to give better agreements with simulations for asymmetric rates especially when combined with averages over b defined as the harmonic mean, i.e., Eq. (3). Because harmonic means are not additive, the summation must precede the average in Eq. (5d) within the current approach. Equation (5) becomes exact for symmetric template models.
In applying the mean field theory expressions, all forward rates are first scaled with k [k z (aDb)?k z (aDb)~k z (aDb)=k], and Eq. (41) is solved for each b to find n(b). Although it is a quartic equation with respect to n(b), there was always only one solution for 0vn(b)v1. The mean velocity n, error rate m, and thermodynamic force F then follows from Eqs. (5c), (5d), and (7). The dependence of n and m on F are obtained by treating g as a parameter to be eliminated (Figure 3) [62].

Limiting behavior of elongation properties
Equilibrium is reached when n~0 and F~0, which occurs when n(b)~0 for all b: from Eqs. (5a), (8), and (9), we verify Eqs. (29) and (31), and the maximum error rate is which would be equal to 3/4 if g~1 and m~4. Far from equilibrium, where F ??, g~? and k { (aDb)~0, Eqs. (5a) and (41) give n(b)? P a k z (aDb):n max (b) and n?n max :k, where k is given by Eq. (2a). The error rate in this limit approaches its lower bound, which defines the inverse fidelity parameter c via Eq. (2b). The analogous limits of m 2 can also be calculated from Eq. (6): and m 2,max~X a=c X b e ½geq(ajb)zgeq(bjc) c (m{1) e geq e geq g ze geq e geq g z(m{2)e geq e geq g 2 where in the second equality in Eq. (46), the average over c was omitted because g(aDb) is symmetric [Eq. (9)] and the three terms in the square brackets correspond to a~b=c, a=b~c, and a=b=c cases, respectively. For m~4, m 2,max~3 =4~m max if g~1, and m 2,max~0 :732 with g~1:73 (dotted red line in Figure 5D).

Stochastic simulation of single-molecule elongation
Numerical simulations of single-molecule RNA replication kinetics were performed with the Gillespie algorithm [58] applied to Eq. (32) [32]. A sufficiently long template sequence fb n g was pre-generated for the simulations using random sequences with equal distributions u(b)~1=4. The initial condition was chosen as n~0 (with rate constants assigned arbitrarily for nv2), and only conditions that lead to positive velocities were considered. For a given set of k z (aDb) [or k(aDb)~k z (aDb)=k], a value of g is chosen, Eq. (8) is used to generate k { (aDb), and simulations are run to obtain n~n=t, where n and t are the length of chain grown and time elapsed, respectively. The error rates m and m 2 , and the thermodynamic force F are obtained by first calculating q(aDb) over the chain and using Eqs. (5d) and (7). This procedure is repeated for different values of g to yield n and m 2 as functions of F . Typically, simulations were run up to t~10 8 k {1 and properties were averaged over the entire chain grown. For poliovirus 3D pol , simuations were also performed using templates generated by repeating the poliovirus sequence [63] (diamonds in Figure 4B).

Stochastic simulation of evolutionary dynamics
In the stochastic form of the Eigen model given by Eqs. (23), the total rate of transformation at a given time t is where P i Q ij~1 was used. In a simulation, a random number 0vz 1 v1 with a uniform distribution is drawn, and determined the time tzDt of the next replication/degradation event. A second random number 0vz 2 v1 was drawn next, which chooses one (k, k~1, Á Á Á ,2N g ) of 2N g reactions (N g replications and N g degradations, where N g is the total number of distinct genotypes present within the population) from Eqs. (23) following Ref. [58]: where a k~r k n k (1ƒkƒN g , replication), rn k (N g z1ƒkƒ2N g , degradation):

& ð50Þ
For a degradation event, the number of replicators is updated as n k{Ng ?n k{Ng {1. For a replication, a progeny genotype j is produced from i~k by attempting to mutate each nucleotide into 3 different bases with probability m i =3, followed by the update n j ?n j z1. Since the total number of possible genotypes 4 M is exponentially large for even moderate values of M, exact enumerations of n i for all possible genotypes was avoided. Instead, the simulation proceeded by first creating from the initial distribution a list of genotypes for which n i w0, and adding newly encountered genotypes to the list as mutations occurred.

Test case for quasispecies dynamics
For testing the Gillespie simulation of RNA population dynamics, we used the single-peak Eigen landscape (18) without back-mutation, for which the quasispecies dynamics (15) can be easily integrated. Although more advanced methods pioneered by Saakian and coworkers [17,[20][21][22]24] allow exact analyses of the Eigen model, the following simple treatment suffices for our purpose of testing numerical simulations because for moderately large M, the effect of back-mutations become negligible. Writing a vector n~½n 1 n 2 T , where n 1 and n 2 are the total numbers of individuals with the master sequence and mutants, respectively, and ignoring back-mutations, Eq. (15) can be written as where R~q r 1 0 (1{q)r 1 r 2 and q~(1{m 2 ) M . Diagonalizing R and integrating, we get n(t)~e Rt : n(0)~V : e l 1 t 0 0 e l 2 t " # : where l 1~q r 1 and l 2~r2 are the two eigenvalues of R and V is the eigenvector matrix. We then obtain the time-dependent master sequence frequency p 1~n1 =(n 1 zn 2 ), p 1 (t)~q A{1 A{1{(1{q)Ae (r 2 {q r 1 )t : ð54Þ

Starvation process
During an idealized starvation process, a single genotype is placed inside a medium containing N (0) m NTPs and N (0) p PPi's with the corresponding initial thermodynamic force F i . As replication progresses, F decreases via F~F 0 z ln N m N p~F 0 z ln assuming rapid mixing, and N grows via _ N N~rN. Equation (55) can be solved for N to give Introducing the fractional population size j with respect to the asymptotic population N ? reached in the limit of equilibrium (F ?0), j: or F~ln Equation (58) with n~n(F ) and m~m(F ) give the dependence of mean fitness on j during a starvation process with the initial condition F~F i . This procedure assumes that the early stages of growth with finite N for which n g deviates from Eq. (19) make negligible contributions. The mean fitness averaged over the process is Sn g T~1 N ?
This integral was performed using trapezoidal rules to obtain Figure 10.