Thermodynamic and Kinetic Sequence Selection in Enzyme-Free Polymer Self-Assembly inside a Non-equilibrium RNA Reactor

The RNA world is one of the principal hypotheses to explain the emergence of living systems on the prebiotic Earth. It posits that RNA oligonucleotides acted as both carriers of information as well as catalytic molecules, promoting their own replication. However, it does not explain the origin of the catalytic RNA molecules. How could the transition from a pre-RNA to an RNA world occur? A starting point to answer this question is to analyze the dynamics in sequence space on the lowest level, where mononucleotide and short oligonucleotides come together and collectively evolve into larger molecules. To this end, we study the sequence-dependent self-assembly of polymers from a random initial pool of short building blocks via templated ligation. Templated ligation requires two strands that are hybridized adjacently on a third strand. The thermodynamic stability of such a configuration crucially depends on the sequence context and, therefore, significantly influences the ligation probability. However, the sequence context also has a kinetic effect, since non-complementary nucleotide pairs in the vicinity of the ligation site stall the ligation reaction. These sequence-dependent thermodynamic and kinetic effects are explicitly included in our stochastic model. Using this model, we investigate the system-level dynamics inside a non-equilibrium ‘RNA reactor’ enabling a fast chemical activation of the termini of interacting oligomers. Moreover, the RNA reactor subjects the oligomer pool to periodic temperature changes inducing the reshuffling of the system. The binding stability of strands typically grows with the number of complementary nucleotides forming the hybridization site. While shorter strands unbind spontaneously during the cold phase, larger complexes only disassemble during the temperature peaks. Inside the RNA reactor, strand growth is balanced by cleavage via hydrolysis, such that the oligomer pool eventually reaches a non-equilibrium stationary state characterized by its length and sequence distribution. How do motif-dependent energy and stalling parameters affect the sequence composition of the pool of long strands? As a critical factor for self-enhancing sequence selection, we identify kinetic stalling due to non-complementary base pairs at the ligation site. Kinetic stalling enables cascades of self-amplification that result in a strong reduction of occupied states in sequence space. Moreover, we discuss the significance of the symmetry breaking for the transition from a pre-RNA to an RNA world.


Introduction
Extant biological systems use different molecules for the storage of genetic information than for the catalysis of biomolecular reactions. Inevitably, the question arises of which came first, the informational polymer carrying the instructions for the enzymes, or the enzymes assembling the polymers? As RNA cannot only store genetic information but also fold into catalytically active structures [1][2][3][4][5], it is central to one of the most prominent hypotheses for the emergence of living systems [6][7][8][9][10][11]. However, this RNA world hypothesis does not explain the origin of the catalytic RNA molecules, called ribozymes [12]. While recent experimental work revealed potential prebiotic pathways to synthesize nucleotides [13][14][15], the mechanisms assembling these building blocks into functional molecules are only beginning to be explored [16][17][18][19][20][21]. The smallest ribozymes known today are 30 to 100 nucleotides long [22][23][24]. More complex ribozymes that could, e.g., assist replication are likely to have a minimum length of more than 150 nucleotides [25][26][27]. For polymers of a length between 30 to 150, a total of 10 18 to 10 90 distinct sequences are possible. However, the subset of catalytically active sequences is generally believed to be tiny. Hence, the spontaneous emergence (and maintenance) of a functional ribozyme in a random pool of oligonucleotides seems highly unlikely. Therefore, it remains unclear how the transition from a pre-RNA world lacking ribozymes to an RNA world where essential processes are driven by ribozymes could occur.
To understand this transition, one must study the dynamics in sequence space, which emerges when mononucleotides and short oligonucleotides inside a reaction volume collectively self-assemble into longer strands [28,29]. The self-assembly is governed by the process of templated ligation [30][31][32][33][34]. In this process, two strands that are hybridized adjacently on a third template strand are covalently linked. In contrast to random ligation concatenating two arbitrary strands, the process of templated ligation is sequence-selective for two reasons, a thermodynamic and a kinetic one. First, complementary nucleotide pairs at the hybridization sites increase the stability of the complex of strands, and thereby also increase the probability for a new covalent bond to be formed. However, complementary sequences of the same length comprising different sequence motifs can have different binding stabilities. The stability of a complex is determined by the hydrogen bonding between complementary base pairs and the stacking interactions between neighboring base pairs [35]. Changing the order of the base pairs or flipping one of the pairs generally alters the complex's stability. Hence, certain sequence motifs can be favored over others thermodynamically [35][36][37]. Second, the kinetics of the ligation step are motif-selective: Non-complementary nucleotide pairs in the vicinity of the ligation site stall the formation of a new covalent bond. As a result, the formation of new strands from shorter fragments that do not match the template strand is also suppressed kinetically [32,[38][39][40].
Experimentally probing the enzyme-free self-assembly of long strands from a random pool of mononucleotides and short fragments is challenging. Typically, the experiments require long times, while the reaction yields remain low and undesired side products obscure the results. Moreover, tracking the evolution of the whole sequence pool simultaneously remains an unsolved technical challenge [20,21,41,42]. Due to these constraints, non-enzymatic self-assembly experiments either employed initial oligonucleotides with precisely designed sequences limiting the product space [43][44][45][46][47] or focused on primer extension scenarios. In the latter scenario, a defined primer that is statically bound to a defined longer template strand gets extended by mononucleotides and short oligomers [38,[48][49][50][51][52][53][54]. Two explorative experimental studies investigated the emergence of progressively longer strands from DNA-oligomers [18,55], both using DNA ligases to accelerate the assembly dynamics and to obtain better yields, and employing temperature cycling for strand separation. In Ref. [55], all possible 12-mers that can be formed from a binary alphabet of A and T are present initially. The assembly dynamics give rise to structured sequence pools characterized by a reduced sequence entropy compared to a random pool. The emerging longer strands are either characterized by a large A or T content since mixed strands are more prone to self-inhibition due to hairpin formation. In Ref. [18], three pairs of carefully designed complementary sequences composed of 20 nucleotides were used as basic building blocks. The authors demonstrated that certain subsets of sequence motifs composed of two basic building blocks form cooperative networks. Since the initial building blocks are already quite long in both studies, the binding energies of bound strands are large, such that small differences in the stacking energies associated with adjacent nucleotide pairs become irrelevant. However, subtle differences in the stacking energies might trigger sequence selection already on the level of the shortest oligomers, i.e., dimers and trimers, for which dissociation occurs spontaneously, rather than being induced externally. Since a sequence bias emerging early on might feedback onto itself, it could have a substantial impact on the pool of longer strands at later stages. In summary, an experimental study exploring growth dynamics into longer polymers starting from a pool of small building blocks is still missing.
Investigating the collective growth from small building blocks theoretically or by means of computer simulations in a model including the essential features of self-assembly, i.e., sequence-dependent (de)hybridization and ligation dynamics, is also challenging: First, the number of possible complex configurations grows exponentially fast as strands become longer. Second, there is an intrinsic separation of time scales between the fast dissociation of short and the slow dissociation of long hybridization sites and the slow ligation step. To date, no theoretical study on self-assembly via templated ligation accounted for the motifdependent thermodynamic and kinetic aspects of hybridization and bond formation (see Section 4.2). Therefore, the following questions remained open: (1) What are the emerging dynamics in sequence space as strands grow longer? (2) Which are critical factors that enable self-enhancing sequence selection? (3) How do motif-dependent thermodynamic and kinetic parameters affect the selection process?
In Ref. [56], we developed a simulation method that partially handled the complexity of the self-assembly process. In this first study, we treated the sequence dependence of the (de)hybridization dynamics in a mean-field picture, in which the dissociation rate only depends on the length of the hybridization site. Our study identified several growth regimes arising from the competition of timescales for dissociation and extension. Moreover, we showed that, depending on external control parameters, the strand length distribution in the stationary state can exhibit a non-monotonous shape characterized by a distinct strand length. For the present study, we extended the simulation method to explicitly treat sequences, including sequence-dependent thermodynamics and kinetics. The 'RNA reactor' simulations that we report here assume a closed reaction volume, initialized with mononucleotides and a few dinucleotides, with an unbiased nucleotide distribution (symmetric initial condition in sequence space). Within the RNA reactor, oligomers grow via templated ligation and degrade via hydrolysis. Eventually, the sequence pool converges to a non-equilibrium stationary state characterized by its length and sequence distribution.
To address the above questions, we consider different model variants. We start with a simple reference scenario, where kinetic stalling is absent and the stacking energies for all complementary neighboring nucleotide pairs are identical. This scenario distinguishes solely between complementary and non-complementary pairings. We then introduce thermodynamic and kinetic sequence selection, both separately and in combination, and compare the resulting four different scenarios. Our main finding is that, under the conditions assumed here, thermodynamic discrimination within hybridized strands is not sufficient by itself to promote self-enhanced sequence selection that drives the sequence pool significantly away from the random state. However, distinct patterns in sequence space arise if non-complementary strand termini at the ligation site slow down the ligation step significantly (kinetic stalling). In this case, a small thermodynamic bias for certain sequence-motifs triggers a self-enhancing dynamics, such that the thermodynamically favored sequence motif dominates the stationary state.

Strands and Complexes
We consider a binary alphabet composed of two complementary nucleotides, denoted as X and Y for generality. A molecule containing L nucleotides linked covalently is called a strand of length L (see Figure 1a). A single nucleotide is a strand of L = 1. Strands are directed and point from the 5 to the 3 end, which we also refer to as the − and the + ends. We allow strands to hybridize to each other, but do not account for the possibility of self-folding. An entity formed by several hybridized strands is referred to as a complex. All staggered conformations that can arise from a set of single strands are allowed inside the RNA reactor, regardless of the number of strands and mismatches, i.e., non-complementary nucleotide pairs (see Figure 1b,c and Figure A1 in Appendix B). However, branched hybridization structures and other nonlinear complexes involving loops are excluded.
We call a complex that contains two or three strands a duplex or triplex, respectively. The overlapping horizontal region between two strands is referred to as a hybridization site. Moreover, the vertical interface between two strands hybridized adjacently on a third strand is called a ligation site. Figure 1. Schematic illustration of the dynamics inside the RNA reactor. The elementary processes are hybridization, dehybridization, ligation on the template, and hydrolysis with corresponding elementary rates k on , k off and k lig , and k cut . The elementary rates k on , k off , k lig are functions of the sequence context: (a) Strands have a binary sequence and are directed. L denotes their length. (b) When two molecules collide, they can form χ different hybridization complexes. (c) Hybridization sites within complexes (horizontal interfaces) can contain mismatches. Two strands (− and + strand) located adjacently on another strand may get joined covalently via templated ligation. The speed of the ligation reaction depends on the complementarity κ of the nucleotide pairs at the ±1 and ±2 position (red box). Non-complementary pairings lead to kinetic stalling. (d) The stability of a hybridization site is governed by the hybridization energy ∆G hyb . ∆G hyb is obtained by summing over stacking energies γ associated with nearest-neighbor blocks (purple box) and considering terminal nucleotide pairs. ∆G hyb and γ depend on the structural and sequential context. Mismatches weaken the binding. (e,f) Covalent bonds within single strands or single-stranded segments may get cleaved via hydrolysis at a constant rate. The resulting unactivated strand termini are assumed to be rapidly reactivated.

Elementary Reactions
Strands and complexes form new complexes via hybridization, dehybridization, templated ligation and hydrolysis (see Figure 1). All reactions are assumed to be elementary and occur with sequence-and structure-dependent rates k on , k off , k lig , and k cut . Assuming constant environmental conditions, k on and k off are related to the hybridization energy ∆G hyb associated with a hybridization site via the thermodynamic consistency requirement [57] where β = (k B T) −1 , k B is Boltzmann's constant, and T denotes the (absolute) temperature, V and N A are the reaction volume and Avogadro constant and c • = 1 mol/L is the reference concentration. We will express all concentrations as a multiple of the reference concentration. Moreover, in the following, we use the dimensionless hybridization energy Γ is obtained by summing over dimensionless motif-dependent stacking energies of nearestneighbor blocks [35][36][37] (see Section 2.4). Γ thus reflects the number of complementary and non-complementary nucleotide pairs and their arrangement. Generally, mismatches increase the hybridization energy, therefore reducing the stability of a complex. The formation of new covalent bonds requires energy, which needs to be provided by the environment in the form of an activation chemistry [42,58,59]. We assume that the RNA reactor is constantly fueled with the activation chemistry, enabling a fast chemical (re)activation of the termini of all strands present in the system. With that, two strands that are located next to each other on a third strand can always ligate. The fast activation step is not modeled explicitly. The rate k lig at which two neighboring strands ligate depends on the paired nucleotides in the vicinity of the ligation site. Mismatches lead to kinetic stalling, i.e., a reduction of the ligation speed [32,[38][39][40].
We model the kinetic stalling using the kinetic stalling factors Φ ± ≤ 1. The stalling factors Φ ± are functions of the complementarities κ ±i ∈ {1, 0} of the paired nucleotides in the vicinity of the ligation site. The value 1 indicates a complementary pair, whereas the value 0 indicates a non-complementary pair of nucleotides. Φ − takes the complementarities κ −1 , κ −2 of the two nucleotides in the − direction of the ligation site into account, while Φ + is an equivalent expression for the two nucleotides in the + direction (see Figure 1c,d and Section 2.5 for more details). The two stalling factors are then multiplied with the basal ligation rate λ. With that, the ligation rate k lig becomes Since random ligation of two strands in the absence of a template is weak compared to templated ligation [31][32][33][34], we neglect it in our model. While covalent bonds within double-stranded parts of complexes are assumed to be stable against hydrolysis, bonds within single-stranded sections get cleaved [60][61][62][63] (see Figure 1e,f). The corresponding rate is assumed to be sequence-independent, k cut = const. (4) with that, the overall degradation rate for a single strand of length L is (L − 1)k cut , for example. In real systems, k cut varies by several orders of magnitude as a function of environmental parameters and crucially depends on the polymer's backbone chemistry [60,[62][63][64][65]. Note that templated ligation and cleavage are irreversible, since the respective reverse reactions (random ligation and "templated cleavage") are absent in our model.

Kinetics of Hybridization and Dehybridization
Since Equation (1) only constrains the ratio of k on and k off , an additional kinetic parameter is required to fix the kinetics of the model. However, the chosen parametrization has only a minor effect on the global kinetics, given that ligation and hydrolysis are rare compared to hybridization and dehybridization (see Section 2.8). Our approach uses a constant rate of collision between two complexes k coll = (VN A c • t 0 ) −1 , where t 0 is the collision time scale. In the following, we express all times in units of collision time scale t 0 .
In general, two colliding complexes can form multiple hybridization configurations via χ distinct hybridization channels (see Figure 1b). We assume no bias for any channel, such that the probability of choosing one particular channel is Hence, the rate for a hybridization via a given channel is whereas the dehybridization rate becomes If χ = 0, no hybridization can occur. This is the case if one of the colliding complexes is a duplex without any overhang. A parametrization attributing the hybridization energy Γ to the dehybridization rate is common in theoretical approaches and is consistent with experiments [66][67][68]. The kinetic model resulting from the specific choice of k on and k off , i.e., Equations (6) and (7) was developed, described in detail, and rationalized in Ref. [56], where we studied self-assembly in a sequence-independent model with hybridization energies simply being proportional to the overlap lengths. The kinetic assumptions Equations (6) and (7) reduce the computational complexity considerably, while still sampling complexes in a thermodynamically consistent way (see Appendix C).

Hybridization Energy
Detailed models for the free energy of given RNA and DNA secondary structures [35][36][37] build on the so-called stacking interactions of neighboring nucleotide pairs at the hybridization sites [35,69]. Every nearest-neighbor interaction, i.e., every block of two adjacent nucleotide pairs, is associated with a motif-dependent stacking energy. These stacking energies additively contribute to the total free energy. Additional contributions to the total free energy take into account nonlinearities of secondary structures such as loops, branching points, and particular end configurations.
Our coarse-grained model that excludes nonlinear complex structures conserves the essential feature of the detailed nearest-neighbor models. The central element of our energy model is the stacking interaction of two neighboring nucleotides pairs P i and P i+1 with where dots symbolize hydrogen bonds between complementary nucleotides. To every block of adjacent nucleotide pairs [P i P i+1 ], we assign a dimensionless stacking energy γ([P i P i+1 ]). (Note that the last two pairs are non-complementary. Therefore, the nucleotides are not connected via a dot.) The hybridization energy is then given by the sum over all stacking energies and contributions − and + accounting for the terminal nucleotide pairs at the − and the + end of double-stranded segment (see Figure 1d), i.e., The contributions ∓ for the ∓ end also depend on the structural and sequence context. If the ∓ terminal nucleotide pair forms a dangling end, i.e., is preceded or followed by an unpaired nucleotide, we have ∓ = 0. If the terminal nucleotide pair is part of a ligation site, there also is a contribution ∓ = 0. If otherwise, it corresponds to blunt end of a complex, and we have ∓ = 0 (see Appendix A for details). For simplicity, we assume symmetric stacking energies, i.e., γ([P i P i+1 ]) = γ([P i+1 P i ]). Moreover, complementary nearest-neighbor blocks are either alternating if or homogeneous if Here, the − symbol stands for a covalent bond. We denote stacking energies assigned to alternating and homogeneous blocks by γ alt and γ hom . Motivated by the observation that γ alt = γ hom in DNA and RNA systems (see Table 1) [36,37], we treat the energy difference as a variable parameter. Without loss of generality, we assume ∆γ ≤ 0 for our model. Moreover, we assume constant stacking energies γ 1nc and γ 2nc for nearest neighbor blocks containing one or two non-complementary nucleotide pairs. Since blocks containing mismatches weaken the binding, their stacking contributions are positive. The contribution for a block with two mismatches is larger than for a block with only one mismatch. In summary, the block-wise contributions obey the hierarchy where γ com is the average energy value of complementary blocks (see Table 2), i.e.,  (12)). Note that the sign of ∆γ depends on whether A and U or T or G and C are considered for the binary system.

Kinetic Stalling
Our kinetic stalling model describes the experimentally observed sequence dependence [32,38,39] in a simplified way, using only two parameters, σ 1 , σ 2 . Mismatches directly at the ligation site affect the ligation speed more substantially than distant ones. If the nucleotide pair at the ±1 position is non-complementary (κ ±1 = 0), a mismatch at the ±2 position (κ ±2 = 0) amplifies the stalling effect. Otherwise, a mismatch at the ±2 position has no effect, i.e., where σ 1 ≤ σ 2 . If the hybridization site in the + or − direction contains only one nucleotide pair (see Figure 1c), we use Equation (15) with κ ±2 = 1. The strength of the stalling effect depends on the underlying activation chemistry as well as the type of nucleotides being used [32,38,39]; therefore, we treat σ 1 and σ 2 as variable parameters (see Table 2).

Effective Cyclic Environment
According to the energy model defined in Equation (9), hybridization energies for long, primarily complementary hybridization sites become arbitrarily negative. Hence, the corresponding dehybridization rates converge to zero exponentially. As a result, strands can be bound in duplexes without single-stranded overhangs over long times. This effect is called template inhibition and leads to freezing of the dynamics [41,61,70]. To overcome this problem, we assume cyclic variations of the physico-chemical conditions (temperature, pH , or salt concentrations) inside the RNA reactor such that all hybridized strands separate within the period time τ [19,71]. Aforesaid oscillatory conditions arise for example due to convection flows induced by temperature gradients or micro scale water cycles at a heated gas-liquid interface. Both scenarios arise naturally in rock fissure in the vicinity of hydrothermal vents [72][73][74][75][76]. They are modeled effectively by introducing a lower bound for the dehybridization rate [56,77], i.e., modifying Equation (7) such that where τ = k −1 low . With that, the (dis)assembly dynamics of long complementary complexes do not obey the thermodynamic consistency requirement Equation (1) anymore. Nonetheless, the kinetics are still plausible: The constant collision rate is a reasonable approximation for a collision process with a diffusion coefficient decaying with length compensated by a cross section growing with length [56].

Validity of Our Model and Application to Primer Extension
In the Results section, we focus on self-assembly scenarios where all strands (apart from monomers) are equally important because there are no distinct template, primer, and substrate strands as in typical primer-extension situations. However, in Appendix F, we show that our modeling of the kinetic stalling and the (de)hybridization kinetics in combination leads to copying dynamics in primer-extension situations consistent with the experimental literature.

Parametrization of Rates
We can parametrize every rate constant k * introduced so far by a dimensionless length l * such that k * = e γ com l * .
This presentation will prove convenient in the later analysis of the results as it connects time scales to length scales. For example, l low = 7 tells us that entirely complementary hybridization sites composed of more than seven nucleotides dissociate as quickly as altogether complementary hybridization sites comprising exactly seven nucleotides. Parameters used in the following are summarized in Table 2. Moreover, l lig = 10 signifies that the timescale of a dehybridization for a hybridization site counting more than ten nucleotides would be slower than the bare ligation timescale if the lower bound with l low would not have been introduced.

Implementation
To simulate the model dynamics in C++, we use an extension of the framework developed in [56], based on an optimized Gillespie algorithm [78][79][80]. The simulation only keeps those species in memory that have a non-zero copy number. If a species appears (vanishes), the corresponding species object is created (deleted) dynamically.

Boundary Conditions and Observables
We aim to investigate the model dynamics starting on the lowest level, i.e., where mononucleotide and a few short oligonucleotides collectively evolve into larger entities. Will the sequences of longer strands be random, or will they show patterns? We chose the arguably simplest setting for our study, which is a closed reaction volume that does not exchange complexes with the environment. In such a setting, we expect the dynamics to settle to a stationary state eventually. We initialize the reaction volume symmetrically with 5000 nucleotides distributed over 4920 mononucleotides and 40 dimers (see Figure 2a). Moreover, we adjust the reaction volume such that the total nucleotide concentration is given by c tot = 0.01 c • . The ratio of the initial monomer to dimer concentration is c init 1 : c init 2 = 123 : 1. The reactor is initialized symmetrically with mononucleotides and a few dimers such that the amounts of X and Y nucleotides are equal and that all four dimer sequences have the same concentrations (see Section 3.1). Within the RNA reactor, oligomers grow via templated ligation and degrade via hydrolysis. Eventually, the sequence pool converges to a non-equilibrium stationary state characterized by its length and sequence distribution (see Figure 1); (b) To characterize the dynamics in sequence space, we introduce the zebraness ζ on the level of single strands and the system-level zebraness Z. The zebraness ζ of a single strand is the fraction of zebra motifs, i.e., alternating binary motifs contained in the strand. In contrast, the system-level zebraness Z measures how zebra-like, i.e., alternating or homogeneous, the pool is as a whole. Z corresponds to the total number of zebra motifs spread over all strands normalized with respect to the overall number of binary motifs within all strands present in the reactor.
Our focus is on the evolution of the length distribution and the dynamics in sequence space. The length distribution c L expresses the concentration of strands of length L, irrespective of whether they are part of a complex or not. We denote the mean length by L. To describe the dynamics in sequence space, we aim for a simple observable with an intuitive and straightforward interpretation. Therefore, we introduce the zebraness as a characterization of a strand's sequence. The zebraness ζ(S) of a strand S of length L S is the number of alternating "zebra" submotifs X − Y or Y − X within its sequence divided by the number of binary motifs L S − 1 (see Figure 2b). With that, a random sequence S r is expected to have ζ(S r ) = 0.5 on average. Moreover, the system-level zebraness Z characterizes how zebra-like the ensemble of strands is. It is given by where the summation is performed over all individual strands with L > 1. A system containing homogeneous strands only would have Z = 0, whereas, for a system exclusively composed of strands with alternating sequences, we would have Z = 1. All plots show ensemble averages which are taken over 20 independent realizations of the dynamics.

Overview of Key Findings
Before we present the detailed analysis of the four different scenarios outlined in the introduction, we briefly summarize our key findings. First, we study the simplest variant of our model where both kinetic stalling and energetic bias are absent, i.e., ∆γ = 0 and σ 1 = σ 2 = 1. This scenario only distinguishes between nearest-neighbor blocks containing zero, one, or two mismatches. Alternating and homogenous blocks have identical energetic properties. While non-complementary pairings decrease the complex's stability, erroneous pairings at the ligation site do not reduce the bare ligation rate. The first model variant does not give rise to motif selection; the composition of the sequence pool remains entirely random, i.e., Z = 0.5.
In the second scenario, we introduce an energetic bias ∆γ < 0 for alternating blocks while still assuming a non-discriminative ligation. The energetic bias favors the hybridization of strands with zebra-like sequences. This time, a weak zebra pattern with Z > 0.5 is induced transiently during the initial growth phase. However, the pattern vanishes almost completely as the system approaches the steady-state (see Figure 3).
The dynamics in sequence space change drastically if kinetic stalling with σ 1 , σ 2 < 1 is applied. If non-complementary nucleotide pairs at the ligation site slow down the formation of a covalent bond, distinct patterns in sequence space can emerge. In the third scenario, we investigate the correlation between the strength of the kinetic stalling effect and the reduction of possible states in sequence space assuming identical energetic properties for alternating and homogeneous blocks, i.e., ∆γ = 0. Within this setting, we observe a spontaneous symmetry breaking in sequence space. Independent realizations of the dynamics evolve to stationary states, dominated by either zebra motifs with Z < 0.5 or homogeneous motifs with Z > 0.5 (see Figure 4). Moreover, we see that a dominant pattern emerges such that Z → 0 or 1 if the stalling effect is strong enough.
In the fourth scenario, we show that a slight energetic bias ∆γ < 0 can become selfamplifying if kinetic stalling is present (see Figure 5). Depending on the strength of the kinetic stalling, the system converges to either a partial or pure zebra state characterized by either Z > 0.5 or Z → 1.

Reference Model without Energetic Bias and Kinetic Stalling
This section aims to answer whether the energetic discrimination of matches and mismatches alone is sufficient to give rise to spontaneous symmetry breaking in sequence space such that Z = 0.5. To this end, we study the simplest variant of our model with neither energetic bias nor kinetic stalling, i.e., ∆γ = 0 and σ 1 = σ 2 = 1.
Initially, the growth dynamics of the mean length L is slow until t ≈ 8.8 × 10 9 (see the dark blue curve in Figure 3a). At this point, the mean length L starts to increase rapidly. We refer to this time point as the onset of growth and denote it by t. After the steep increase, L reaches a plateau value. The inset shows the steady-state length distribution displaying a double-exponential shape.
The ensemble average of the zebraness Z initially fluctuates and then converges to Z = 0.5 (see Figure 3b). Looking at single trajectories (see Figure 3c) reveals a behavior similar to the ensemble average. The initial values of Z ≶ 0.5 on the single trajectory-level are due to small numbers of strands with L > 1. A value of Z = 0.5 hints towards an entirely random sequence pool but does not exclude motif correlations on larger scales. However, analyzing distributions of longer motifs reveals that the final sequence composition is indeed random (see Appendix D).  (19). For ∆γ < 0, L reaches a maximum before decaying gradually to the stationary value. The inset shows the steady-state length distributions. (b) If there is no energetic bias, i.e., ∆γ = 0, no distinct patterns emerge in sequence space, and hence Z = 0.5 (see also Appendix D). If an energetic bias ∆γ < 0 is applied, Z grows initially and then decays when L ≈ 7. The final value is slightly above the random state Z = 0.5 and below the simple thermodynamic estimate Z * (see Equation (22)). (c) Single realizations of the dynamics for ∆γ = 0 behave similar to ensemble average. Strong fluctuations for small times stem from low numbers of strands with L > 1. (d) The fraction of mismatches m first decreases and then increases as the mean length becomes longer. (e) The fraction of concealed mismatches, i.e., mismatches not affected by energetic discrimination grows simultaneous with the mean length. (f) Over time, concealed erroneous ligations become frequent and destroy the initial sequence bias.
The evolution of the mean length shows some interesting features. After a lag phase, its increase becomes exponential at t = t (see Figure S1 in the Supplemental Material). Formally, we define the onset of growth t by intersecting the tangents to the L-curve at t = 0 and the point where the increase is strongest (dashed line in Figure 3a; for details, see Appendix E and Figure S2). We observe that t coincides with the moment in time at which higher-order ligations, i.e., ligations involving at most one monomer become more abundant than ligations joining two monomers to a dimer (see Figure A3 in Appendix E).
Moreover, we can predict the onset with a relative error smaller than 15% by the following formula derived in Appendix E (dotted line in Figure 3).
Here, c init 2 is the initial dimer concentration and k {1,1|2} and k {1,2|2} are the effective rate constants at which new dimers or trimers are formed from monomers or monomers and dimers. k {1,1|2} is given by where K {1,1|2} is the effective dissociation constant averaged over all triplexes involving two monomers and one dimer. An analogous expression exists for k {1,2|2} . Repeating the computer experiment with k cut and k lig different from the standard values given in Table 2 confirmed the validity of Equation (19) (see Figures S3 and S4). Moreover, altering k cut and k lig while keeping the other parameters fixed confirmed that the mean length L in the stationary state does not explicitly depend on these two variables but only on their ratio as expected from dimensional analysis. The dependence of the mean length on the ratio k lig k cut can be derived analytically for a random ligation model [81].

Energetic Bias in the Absence of Kinetic Stalling
In the previous section, we saw that the energetic discrimination between complementary and non-complementary nucleotide pairs alone is insufficient for the spontaneous emergence of patterns in sequence space. Therefore, we now ask whether an energetic bias ∆γ < 0 favoring the binding of zebra motifs can induce zebra patterns that become self-amplifying, while kinetic stalling is still absent (σ 1 = σ 2 = 1).
The energetic bias, ∆γ < 0, causes a transient overshoot in the mean length beyond the steady-state value but otherwise does not strongly affect the dynamics of the mean length ( Figure 3a). The onset of growth t can still be predicted by a formula analogous to Equation (19) (see Appendix E and Figures S5-S7 for plots of single realizations). The steep increase after the lag phase is followed by a gradual descent to the steady state value, slightly below the maximum. Moreover, the steady-state length distribution remains very similar to the scenario without energetic bias.
In sequence space, a simple thermodynamic estimate Z * for the final zebraness can be made based on a two-state system with an energy difference ∆γ Since this estimate neglects any correlation and feedback effects, one could naively expect that the zebranass resulting from the simulated model dynamics reaches a value larger than Z * . Indeed, the observable initially grows beyond the estimate Z * . However, the growth stops when the mean length reaches a value of L ≈ 7. At that point, the zebraness starts to decay and converges to a stationary value simultaneously with L (see Figure 3b). In the stationary state, the zebraness Z is only slightly above the result for random sequences, i.e., Z = 0.5, and below the simple thermodynamic estimate Z * .

Loss of Energetic Discrimination Prevents Sequence Selection
Why are the initially emerging zebra patterns triggered by the energetic bias ∆γ < 0 in Figure 3b neither amplified nor maintained? In the following, we analyze the growth processes in detail to give an intuitive explanation.
Strand growth requires the formation of complexes comprising at least three strands. The more negative the hybridization energies of the hybridization sites in these complexes are, the more stable the configurations become and the higher the probability for a ligation gets. Non-complementary nucleotide pairs increase the hybridization energy and, therefore, weaken the binding. To analyze the effect of these mismatches, we define the overall fraction of mismatches m as where N pairs and N non are the absolute numbers of nucleotide pairs and non-complementary nucleotide pairs in all present complexes. Initially, the fraction of mismatches m decreases since complexes mostly containing complementary nucleotides that emerge during the early growth persist longer and hence contribute more substantially to the average. However, the fraction of mismatches starts to increase when the mean length becomes larger (see Figure 3d). This increase of the mismatch fraction arises from the loss of thermodynamic discrimination induced by the cut-off k low in the dehybridization rate, i.e., the effective temperature cycles (see Section 2.6). Although the hybridization energy may become arbitrarily negative for large hybridization sites, the dehybridization rate can not become smaller than the lower bound k low . The length scale associated with k low is l low = 7 (see Section 2.8). This implies that an entirely complementary hybridization site comprising more than seven pairs has the same stability as a mismatch-free hybridization site composed of exactly seven pairs. Moreover, mismatches in extended hybridization sites might have no effect on the rate for unbinding because the hybridization site still contains a high number of matches. If many matches are present, the hybridization energies are strongly negative such that the lower threshold still determines the rate for dehybridization. This effect enables concealed mismatches. Concealed mismatches are mismatches that do not increase the dehybridization rate k off of a hybridization site. Replacing a concealed mismatch with a complementary pair would not decrease k off further since it is already given by the cut-off, i.e., k off = k low . The longer the strands become during the first growth phase, the more concealed mismatches emerge. With the absolute numbers of mismatches and concealed mismatches in all present complexes N non and N con , we now introduce the fraction of concealed mismatches n con as The evolution of n con shown in Figure 3e reveals that most of the occurring mismatches are concealed, once L has become approximately twice as large as l low . Concealed mismatches also occur at the strand termini at ligation sites and may lead to the formation of new binary motifs which are not complementary to the templating motif at the ligation site. We call such a ligation involving at least one concealed mismatch a concealed erroneous ligation. Dividing the number of concealed erroneous ligations N err per time by the overall number of ligations N lig per time gives the fraction of concealed erroneous ligations n err con , i.e., Every erroneous concealed ligation mitigates the present bias in sequence space and leads to randomness. Erroneous concealed ligation is the reason why the initial sequence patterns decay almost to the random level Z = 0.5. However, not all hybridization sites, particularly the shorter ones, have a dehybridization rate determined by the lower bound. As the initial bias for binary zebra motifs on the system level decreases and sequences become more random, non-concealed mismatches in shorter hybridization sites become more likely. Hence, short hybridization sites not yet affected by the lower bound for the unbinding rate become less stable on average and contribute less to the growth process. This explains why the small maxima in the mean length and the other observables shown in Figure 3a,d-f disappear as the bias for binary zebra motifs fades away.

Kinetic Stalling in the Absence of Energetic Bias
We have seen above that concealed erroneous ligations suppress sequence selection during the self-assembly process. However, in the presence of kinetic stalling, concealed erroneous ligations should be reduced. In this section, we thus investigate the (more realistic) model variant, where kinetic stalling is included. We vary the strength of the kinetic stalling factor σ 1 from 0 to 0.1, while fixing σ 2 = 0.1.
Initially, the dynamics of the mean length L are qualitatively similar to the systems without kinetic stalling studied before (see Figures 3a and 4a). However, the onset of growth t appears later. The time point of the onset t can be predicted by a formula analogous to Equation (19), which considers the kinetic stalling effect, with an error <15%. For σ 1 = 0.05, the values for t from the prediction and the formal definition are highlighted by the dotted and dashed lines in Figure 4a (for details, see Appendix E and Figures S8-S12). On larger timescales, the model including kinetic stalling deviates from the earlier model. After the steep increase, L does not directly settle to a steady-state. Instead, it grows gradually and converges to a constant value eventually. Visualizations with a linear xand a logarithmic or linear y-axis reveal that the initial increase after the lag phase is approximately exponential, while the increase during the second growth phase is approximately linear (see Figures S13-S15). For σ 1 = 0, 0.05 or 0.067, similar stationary mean lengths are reached. However, the relaxation time increases with σ 1 , such that it takes more than ten times longer for a system with σ 1 = 0.067 (see inset of Figure 4a) to converge to the stationary state than for a system with infinite stalling. For σ 1 ≤ 0.067, the steady-state value of the mean length is more than twice as large as for the σ 1 = σ 2 = 1 scenario. Moreover, the length distributions in the stationary state look qualitatively similar to the ones seen earlier (see Figure S14). For σ 1 = 0.1, the increase of the mean length during the second growth phase is small during the time window of observation. From Figure 4a, we can not deduce whether the mean length already approached a stationary value, or whether it will keep on growing. If a stationary value was reached, it would be significantly smaller than for σ 1 = 0.67, 0.05, 0. For σ 1 = 0.1 (as well as for σ 1 = 1), the simulation times are large and prevented us from analyzing at longer time scales. However, plotting the curve for σ 1 = 0.1 in a coordinate system with a linear x-axis might suggest that the system has indeed already converged to a stationary state (see Figure S13). We will discuss the behavior for σ 1 = 0.1 in more detail in Section 3.9 and provide further evidence why the behavior, in this case, might be qualitatively different from the behavior for σ 1 = 0.67, 0.05, 0. Is the novel behavior of the mean length shown in Figure 4a related to a novel motifselective dynamics in sequence space? We now investigate the evolution of the strands' sequences. Since the initial pool of sequences is symmetric and since neither zebra nor homogeneous binary motifs are preferred energetically, we do not expect a preference for a single realization to go to either a zebra (Z > 0.5) or a non-zebra (Z < 0.5) state. Hence, the system-level zebraness Z is not appropriate to describe an ensemble of realizations. As a meaningful observable to quantify the sequences bias on the ensemble level, we, therefore, introduce the sequence order parameter P as During the first growth phase of L, a bias P > 0.5 is established for all values of σ 1 . The dominance of the bias correlates with the strength of the kinetic stalling. During the slow second growth phase, P gradually increases and reaches a stationary value simultaneously with L (see Figures 4b and S11). For σ 1 = 0, 0.05, or 0.067, we observe a value of P ≈ 1 in the stationary state. Hence, on the realization level, the final pool contains either (almost) entirely alternating or homogeneous sequences. We classify such states in sequence space as pure. For σ 1 = 0.1, the final sequence composition within the observation time window is also is non-random. However, the patterns are not pure, i.e., 0.5 < P < 1. We refer to these states as partially mixed states. Whether the system has already converged to a stationary state with P < 1 or will further evolve to a pure state as for σ 1 = 0.67, 0.05, 0 remains unclear at this point (see above). Figure 4c displays the evolution of the zebraness of all realizations forming the steady-state for σ 1 = 0.05. On average, one half of the realizations evolves towards the Z = 1 state, while the other half evolves towards the Z = 0 state. Hence, the symmetry of the initial state is broken spontaneously: either zebra or homogeneous motifs are selected. (See Figures S8-S12 for equivalent plots of single realizations for other σ 1 values.)

Hydrolysis and Stalling Boost Sequence Selection
The previous section revealed a coupling between sequence selection and two distinct growth phases in the kinetic stalling scenario. We now interpret and explain this coupling. Here, we consider the cases σ 1 = 0 and σ 1 = 0.05, where all trajectories eventually converge to a pure state. The case where σ 1 = 0.1 is discussed in Section 3.9.
On the level of individual trajectories, fluctuations lead to a small bias in the motif composition even before the mean length starts to grow rapidly (see Figure 4c). This early asymmetry in the distribution of alternating and homogeneous motifs governs the fate of the realization as seen from Figure 4c.
During the first growth phase (see Figure 4a), monomers and short strands selfassemble into longer strands. The first growth phase ends when most of the initial monomers are consumed. At the end of this initial growth phase, a significant bias towards alternating or homogeneous motifs is present. However, a considerable fraction of binary motifs does not yet reflect system-level bias, i.e., differs from the dominant binary motifs (which are either X − Y and Y − X or X − X and Y − Y). Therefore, mismatches in bound strands are still frequent (see Figure 4d). Moreover, since the average strand length is already significantly above l low , most mismatches are concealed (see Figure 4e).
When monomers and short strands do not dominate the pool anymore, hydrolysis becomes important. Every time a strand breaks, an existing binary motif vanishes. During the second growth phase, fragments of broken strands are reassembled to longer strands via ligation on a template strand. (For an analysis of sequence patterns of strands of specific lengths, see Figure S17). If the kinetic stalling is strong, the ligation of two strands is (almost) impossible if a mismatch occurs at the ligation site and the fraction of concealed erroneous ligations is (close to) zero (see Figure 4f). Hence, every ligation forms a new binary motif that (almost) always complements the templating motif at the ligation site. If the templating motif is zebra-like (homogeneous), the new motif is zebra-like (homogeneous) too. Over time, all binary motifs created during the initial growth phase, particularly those that do not reflect the system bias, get destroyed at a uniform rate. At the same time, binary motifs emerging during the second growth phase likely reflect the system bias. Consequently, the bias becomes self-enhancing (see Figure 4b). Newly created binary motifs enhance the system bias even more, while all motifs that do not reflect the asymmetry in motif space become extinct eventually. As a result, the motif composition becomes more and more ordered and mismatches become rarer (see Figure 4d). Remaining mismatches are now even more unlikely to affect the dehybridization rate since the rate is determined by the lower bound k low for long and primarily complementary hybridization sites. Hence, the fraction of concealed mismatches increases slightly (see Figure 4e). Moreover, if kinetic stalling is finite (σ 1 > 0), the fraction of concealed erroneous ligations also slightly increases (see Figure 4f).
Since the energetic properties are symmetric, every realization randomly approaches either the zebra-like or homogeneous state. As the pool becomes more complementary in its motif composition, triplex configurations achieve higher stability on average. This increase of stability enhances the probability of templated ligations. Therefore, the symmetry breaking in sequence space is concomitant with enhanced growth, which becomes apparent in the further increase of the mean length. For σ 1 = 0 and σ 1 = 0.05, the second growth phase ends when the system has reached an (almost) entirely alternating or homogeneous state (see Figure 4a,b). At that point, the fraction of mismatches is close to zero. The few mismatches that still occur are mostly due to strayed mononucleotides and short oligomers sitting on longer strands. Since these mononucleotides and short oligomers are far from being affected by the lower bound and unbind quickly, the fraction of concealed mismatches takes small values again.

Energetic Bias in the Presence of Kinetic Stalling
The previous section revealed that spontaneous motif selection occurs as a result of kinetic stalling. Without energetic bias, the sequence pool converges to a stationary state which is either dominated by homogenous or alternating sequences. In addition, the energetic symmetry can be broken explicitly if a small energetic bias favoring zebra motifs is applied. To understand the emergent phenomena in this setting, we study two systems, one with strong and one with weak kinetic stalling for various energetic biases ∆γ < 0.
First, we consider the case where σ 1 = 0.05. Most parts of the description in the previous section (∆γ = 0) also apply here (see Figure 5a,b). Again, we can predict the onset of growth t (for details, see Appendix E and Figures S18-S20). The steady-state value of L depends weakly on the energetic bias. However, the final state is reached earlier if the bias is stronger. In sequence space, all trajectories end in a pure zebra state Z = 1. Hence, symmetry breaking in sequence space is now induced energetically as expected because of the explicit symmetry breaking in the energy landscape.
Second, we investigate a scenario with σ 1 = σ 2 = 0.1. The mean length grows strongly in the beginning as before (see Figure 5c and, for more details, see Appendix E and Figures S21-S23). The fast growth phase is followed by either a marginal increase (∆γ = −0.1, −0.2) or decrease (∆γ = −0.3) of L to a stationary value correlating with the strength of the energetic bias. A strong zebra pattern Z > 0.5 is induced in sequence space during the initial increase of the mean length (see Figure 5d). While L grows (decays) during the second phase, Z also grows (decays). Eventually, the sequence pool converges to a partially mixed stationary state with a significant majority of zebra motifs such that 0.5 < Z < 1. The excess of zebra motif again correlates with the strength of the energetic bias. Moreover, from Figures S21-S23, it becomes clear that all single trajectories behave similarly to the ensemble mean, i.e., show steady-state values of the zebraness above 0.5.

Weak versus Strong Kinetic Stalling
In Section 3.6, we speculated whether the system with σ 1 = 0.1 and without energetic bias (blue curve in Figure 4) reaches a stationary state characterized by P < 1 in contrast to the scenarios with σ 1 = 0.067, 0.5, 0, where P = 1 and referred to the corresponding plot with a linear x-axis (see Figure S13). However, this plot did not allow for a clear conclusion either. It could be that the alleged partially mixed stationary state is only transient and that a pure state is reached on much larger time scales. Though the findings from Section 3.8 suggest that the stationary state of the system with σ 1 = 0.1 and without energetic bias is indeed qualitatively different from the scenarios with σ 1 ≤ 0.067 and without energetic and characterized by P < 1. The curves for σ 1 = 0.1 and various energetic biases clearly converge to stationary states with P < 1 (see colored curves in Figure 5c,d). If the stationary state is partially mixed in the presence of an energetic bias, it is not too far-fetched to assume that it is also partially mixed if the energetic bias is absent.
Naturally, the question arises whether a critical value for σ 1 exists above which the system always reaches a pure state with P = 1 for a given value of the energetic bias ∆γ. This first question directly leads to a second question, namely, what would be the nature of the corresponding non-equilibrium phase transition? We leave the answer to this question open for future research. However, finding an answer might be challenging since the relaxation time to the stationary state will probably diverge. At this point, we content ourselves with hypothesizing that two different regimes might exist without drawing an exact border: For strong kinetic stalling, the system converges to a pure state, while, for weak kinetic stalling, it converges to a partially mixed state.
In the light of the above hypothesis, we now analyze the dynamics of mismatches and concealed erroneous ligations in the energetically unbiased system case for σ 1 = 0.1 (see the blue curve in Figure 4d-f as for σ 1 = 0 and σ 1 = 0.05. The self-amplification of the dominant binary motif comes to a halt during the second growth phase (see Figure 4a Consequently, the fraction of concealed mismatches does not decrease again as for σ 1 = 0 and σ 1 = 0.05. Since most of the strands are long enough, most of the mismatches are concealed. Moreover, concealed erroneous ligations are not fully suppressed and still occur in the stationary state. Hence, the weak kinetic stalling scenario includes features of both dynamic regimes described in Sections 3.5 and 3.7.

Summary
Our study investigated the self-assembly of prebiotic polymers via templated ligation inside a non-equilibrium RNA reactor. We identified kinetic stalling as a critical factor for self-enhanced sequence selection. The final sequence space shows no sequence patterns if the underlying stacking energies are uniform without kinetic stalling. In contrast, spontaneous symmetry breaking occurs for strong kinetic stalling. The final pool contains either entirely homogenous or zebra-like sequences. In scenarios without kinetic stalling, any energetically induced sequence bias vanishes almost completely as strands grow. In contrast, in the presence of kinetic stalling, subtle differences in the stacking energies trigger cascades of self-amplification, leading to highly ordered sequence pools. Our results hint towards the existence of two different stalling regimes. We hypothesize that, for strong kinetic stalling, the system converges to a pure state with P = 1, while, for weak kinetic stalling, it converges to a partially mixed state characterized by P < 1.
Initially, the mean length shows burst-like growth dynamics after a short lag phase. The onset of the rapid growth coincides with the time point where higher-order ligations become more abundant than ligations joining two monomers and can be predicted analytically.

Prior Work, Our Model, and Future Extensions
In an earlier model for prebiotic self-assembly, strands only grow via random ligation [82]. There, self-folding and complex formation introduced a protection mechanism against hydrolysis for double-stranded segments. Moreover, the ensemble of strands was assumed to reach a binding equilibrium immediately after a random ligation occurred. In this model, protection against hydrolysis could extend the system's sequence memory. However, the effect was only transient, and all selected patterns vanished eventually. A more recent study combines random ligation and protection against hydrolysis with growth via templated polymerization by mononucleotides [83]. The authors demonstrate that polynucleotides exceeding lengths of 100 can emerge under plausible conditions. Moreover, the authors show that a considerable fraction of the emerging strands forms ribozyme-and tRNA-like secondary structures using folding software. However, the study does not investigate correlations in sequence space on the system level.
Previous theoretical studies considering growth via templated ligation generally explored effective models that reduce the state space to (sub-)sequences without considering complex formation explicitly [33,77,[84][85][86][87][88][89][90][91][92]. Such approaches do not treat (de-)hybridization and ligation as elementary steps. Instead, the reactions are coarse-grained into one extension process. The specification of the corresponding rate neglects the intricacies of the assembly mechanisms and requires a priori assumptions regarding the relevant configurations [33,[84][85][86][87]90]. Moreover, many models ignore that the hybridization energy is a function of the number and nature of the paired nucleotides and use constant (de-)hybridization rates [33,[84][85][86][87][89][90][91]. Other studies treat the sequence dependence of (de-)hybridization employing mean-field approximations where sequence correlations are dismissed [77]. Such simplifications result in systems effectively containing only one type of self-complementary nucleotide [56] and any form of sequence selection is necessar-ily absent. In contrast, our stochastic approach explicitly takes the sequence-dependent thermodynamic and kinetic aspects of templated ligation into account.
Although being already quite complex, our model also made simplifying assumptions. Future studies need to relax some of our assumptions. In particular, one has to consider nonlinear complexes containing loops and multiple branches. Such configurations can give rise to self-folding, self-templating, template inhibition, and gelation [55,82,93]. All these features potentially influence the sequence dynamics. However, we expect that these effects only become important in the long time limit once the strands have reached sufficient size for the formation of secondary structures. Consequently, the discussion of the emergence of structured sequences on shorter timescales in the limit of strong kinetic stalling is expected to hold, even if secondary structures are taken into account. In addition, one has to extend the alphabet size from two to four. The question here is whether pure states containing only a minimal number of sub motifs exist. In the future, the model could also include additional reactions such as non-templated polymerization and ligation, and recombination [82,89,[94][95][96][97][98] or length selective environments [99]. The first two reactions probably play a role in the formation of the first short oligomers, whereas a flow-through system preferentially accumulating long strands can escalate polymerization [100]. Moreover, our study assumed a well-mixed system. Introducing a spatial component together with size-dependent diffusion constants as in Ref. [101], one could study under which conditions local clusters of sequences with specific patterns can emerge and coexist.

Significance for the Emergence of Life
What is the origin of the first ribozymes heralding the transition from the pre-RNA to the RNA world? In Darwinian evolution, the assembly of low-level building blocks into higher-level entities triggered significant developments [134]. In the light of this evolutionary principle, a multi-step process towards greater complexity, eventually resulting in functionality, also seems natural in prebiotic evolution. Here, we studied one of the first steps following the emergence of early nucleotides. This step forms oligonucleotides displaying distinct sequence patterns that could serve as building blocks for the next higher level of self-organization towards functional ribozymes.
In our study, we considered model variants with and without kinetic stalling. Since kinetic stalling is probably inevitable in non-enzymatic templated ligation [32,[38][39][40], the nostalling variant may appear unwarranted. However, this model variant is essential to separate the effects and identify kinetic stalling as a crucial mechanism enabling selfenhancing sequence selection (see Section 3.5). Moreover, the strength of the stalling effect depends on the underlying activation and nucleotide chemistry [32,38,58] and both weak and strong kinetic stalling scenarios, potentially leading to qualitatively different outcomes, are plausible (see Sections 3.6, 3.8 and 3.9). Furthermore, in a pre-RNA world scenario, a primitive ribozyme catalyzing ligations might have a poor ability to discriminate mismatching ends kinetically. In this case, the ribozyme would operate in a regime where thermodynamics mainly control the discrimination between complementary and noncomplementary nucleotides. This regime would be close to the no-stalling model variant.
Moreover, our study revealed that minor differences in the motif-dependent stacking energies significantly affect the dynamics in sequence space. The experiments of Refs. [18,55] probably did not capture this effect, due to the design of the respective experimental systems. These studies used DNA 12-or 20-mers as basic building blocks and a ligase to catalyze bond formation. The ligase requires significant overlaps of both strands with the template to work efficiently. Moreover, the experiments were performed under temperature cycling. The applied temperature cycling was too fast for the long strands to reach a binding equilibrium. Therefore, the hybridization timescale of mostly complementary hybridization sites is always set by the duration of the cold phases. Hence, subtle variations in the stacking energies are not visible. In contrast, the effective cycling in our model is slow enough for thermodynamics to govern the (de-)hybridization of short strands leading to an amplification of the stacking bias (see Section 3.8). The cut-off of the dehybridization rate only affects longer strands emerging from the pool that is already biased.
In non-enzymatic scenarios involving kinetic stalling, the strands formed from the initial pool are already the result of a primary selection process. Selection is not imposed externally but stems from a self-organizing replication network [135]. We showed that the ability to form self-organizing and self-amplifying replication networks is inherent in template-directed growth and does not require higher-level mechanisms such as sequencespecific template inhibition as a result of self-folding, reported previously [55]. In the 1980s, Kauffmann promoted the concept of autocatalytic sets-sets of molecules that mutually catalyze their formations [136]-as a chemical intermediate on the way to biological life [137]. Since then, autocatalytic sets have been the subject of many theoretical and experimental studies [29,30,[138][139][140][141][142][143]. Once our system has reached a stationary state in the strong kinetic stalling regime, it shows the key features of such an autocatalytic set: strands with a specific pattern promote the formation of new strands of arbitrary length, showing the same pattern. Importantly, this concrete realization of an autocatalytic set emerges naturally from an unstructured initial pool without requiring any external (pre-)engineering. This insight could bridge the gap between strand formation and self-sustaining sequence replication.

Data Availability Statement:
The datasets that support the findings of this study, as well as the computer code used to generate these data, are available from the corresponding author upon reasonable request.

Acknowledgments:
We thank Dieter Braun, Annalena Salditt, and the members of the Gerland group for valuable feedback on this manuscript.

Conflicts of Interest:
The authors declare no competing interest.

Appendix A. Terminal Base Pair Energies
The energies j (with the index j denoting to either − or + end) assigned to terminal nucleotide pairs are functions of the complex structure and the sequence context as well. If a terminal nucleotide pair coincides with the end of a complex, a so-called blunt end, there is no additional contribution, i.e., j = 0. In contrast, a terminal nucleotide pair followed or preceded by an unpaired nucleotide, a so-called dangling end, adds a non-zero contribution to hybridization energy, i.e., j = 0 [36,37].
Again, a mismatch results in an energetic penalty j > 0, weakening the binding, while a complementary nucleotide pair leads to a reward j < 0 increasing stability. As for the stacking energies γ([P i P i+1 ]), we distinguish between complementary alternating terminal configurations such as and complementary homogeneous terminal configurations as, for example, We denote the associated energies by alt and hom . As before, we treat the energy difference as a variable parameter. Moreover, for a dangling end configuration involving one mismatch, we assume a constant contribution 1nc . The terminal nucleotide pair contributions obey the inequality with com = ( alt + hom )/2. Parameter values used in Section 3 are summarized in Table A1. For simplicity, we set the energy difference between alternating and homogeneous dangling end contributions to half of the value of the difference between full alternating and homogeneous blocks, i.e., Next, we consider the contribution to the hybridization energy resulting from a terminal nucleotide pair, part of a ligation site. While, in principle, one can choose the contributions due to dangling and blunt ends freely (within a reasonable range), contributions coming from ligation sites are constrained. The reason is that one ligation site is involved in two hybridization sites. To formulate this constraint, we introduce the total hybridization energy β∆G tot (C) of a complex C. β∆G tot (C) is obtained by first summing over the stacking energies of all continuous nearest neighbor blocks within all hybridization sites of the complex, i.e., all blocks of the form where the nucleotide pairs are linked by two covalent bonds (represented by the − symbols). Next, we add the sum over all dangling end contributions. In the last step, we add an energy contribution for every ligation site, i.e., for every noncontinuous nearest neighbor block of the form where one covalent bond is missing. Since the covalent bond does not affect the stacking interaction, the energy contribution to the total hybridization energy β∆G tot (C) of noncontinuous blocks is the same as for the corresponding blocks [37]. For γ as a function of the (non)continuous nearest neighbor blocks, we have, for example, with that, β∆G tot (C) reads Note that we would obtain an identical total hybridization energy if we would replace all non-continuous nearest neighbor blocks in the complex C with continuous blocks, i.e., if we would join the + and − strands at all ligation sites. We now consider two complexes C 1 and C 2 reacting to a new complex C 3 containing one or two new ligation sites. The energy difference between the states before and after the reaction has to correspond to the hybridization energy Γ new associated with the newly formed hybridization site, i.e., Γ new could also be interpreted as the energy that is needed to reverse the reaction. The constraint formulated via Equation (A10) defines the energetic contribution of the newly formed ligation site(s) to Γ new associated with the newly formed hybridization site. In practice, the constraint Equation (A10) reduces to the difference between the stacking contribution(s) γ Newly formed ligation sites also affect the hybridization energies of hybridization sites that were already existing in C 1 or C 2 before the reaction occurred. Hence, these hybridization energies have to be recalculated comparing the total hybridization energies of the new complex C 3 and the compounds that result from virtually dissolving the hybridization site. This procedure is explained in more detail in Appendix B, where we explicitly derive hybridization energies in several exemplary complex configurations.
Note that, in general, This inequation becomes an equation only if the complex C does not contain any ligation sites. Moreover, we can interpret β∆G tot (C) as the total energy stored within the complex. It is equal to the energy difference between the fully disassembled state, where we only have single strands.

Appendix B. Examples for Hybridization Energies
To illustrate the calculation of the hybridization energy, we consider four complexes, sketched in Figure A1.
In example Figure A1a, the duplex C 1 , and the single strand C 2 react to the triplex C 3 . Before the reaction, the total hybridization energies of the complexes are After the reaction, the total hybridization energy is β∆G tot (C 3 ) = hom + γ alt + alt + alt + γ alt + 2γ hom (A14) As stated above, the hybridization energy Γ new associated with the new hybridization site is the difference between the total hybridization energies before and after the reaction, i.e., Since the reaction did not lead to a new ligation site, Γ new corresponds to the sum over all nearest-neighbor blocks plus the danging end contributions. In example Figure A1b, two single strands C 1 and C 2 form a new duplex C 3 . The total hybridization energy before the reaction is zero. Hence, the hybridization energy Γ new of the new hybridization site equals the total hybridization energy β∆G tot (C 3 ) after the reaction, i.e., Again, since no ligation sites are involved, Γ new is equivalent to the sum over stacking and end contributions. In example Figure A1c, the emerging complex C 3 features two new ligation sites. The total hybridization energies of the initial complexes C 1 and C 2 are The total energy of complex C 3 is As before, the difference between the total hybridization energies before and after the reaction determines the hybridization site energy Γ new of the new hybridization site, In the last example, Figure A1d, a mononucleotide C 1 hybridizes to the duplex C 2 resulting in the triplex C 3 . The hybridization energies before and after the reaction are β∆G tot (C 2 ) = alt + γ alt + 2 hom , (A22) Consequently, the hybridization energy Γ new assigned to the mononucleotide is In the examples Figure A1c,d, we also have to recalculate the hybridization energy of the hybridization sites that were present before the reactions since they now also involve an energetic contribution from a ligation site. To this end, we virtually dissolve the hybridization site(s) that were already existing and then virtually reassemble these hybridization sites again. Reassembling the hybridization site(s), we apply the same procedure of calculating the hybridization energy as described before. The updated hybridization energy (energies) obtained that way now include(s) the correct contribution of the newly formed ligation site(s). Moreover, in the examples Figure A1a,c,d, we have to reconsider the channel factors associated with already existing hybridization sites. To update these channel factors, we again virtually disassemble the complex into its compounds (see Figure A1e). We then virtually reassemble the parts and thereby count the number of different configurations that would be possible. With the updated channel factors and hybridization energies, we recompute the dehybridization rates according to Equation (7).

Appendix C. Thermodynamics of Hybridization
With the elementary rates defined in Equations (6) and (7), the total free energy ∆G tot (C) of a complex C is found to be for constant environmental conditions [56]. The first term on the right-hand side is the total hybridization energy defined in Equation (A9), and the second term is a symmetry penalty that occurs if the complex is rotationally symmetric (ρ = 1) and is zero (ρ = 0), otherwise. The free energy β∆G tot (C) is linked to the dissociation constant K D occurring in a mass-action approach where all concentrations are expressed in units of the standard concentration c • via [57] β∆G tot (C) = ln(K D ).
Thermodynamically, the symmetry penalty can be understood as a decrease in the (standard internal) entropy by a factor of ln(2) due to the rotational symmetry. Kinetically, it is rationalized by looking at the interaction probability of two complexes belonging to the same species versus the interaction probability of two complexes representing different species. For equal concentrations, complexes representing different species interact twice as often as complexes belonging to the same species. While the complex resulting from a collision between distinguishable complexes is never symmetric, the interaction of identical molecules always leads to a complex with a rotational symmetry [56]. Hence, the ρ ln(2)-term arises naturally in any collision based kinetic model [81]. We emphasize that the symmetry penalty is due to a reduced product-formation rate rather than a decrease in stability of the complex.
Moreover, the symmetry penalty also appears in the standard databases for free energies of hybridized oligonucleotides [36,37]. These databases also add a constant initiation penalty to the total free energy β∆G tot (C). The initiation penalty is a constant multiplied by (n − 1), where n is the number of strands forming the complex. The initiation penalty accounts for the loss of system entropy due to the fusion of separate entities into one new complex. However, this penalty term can be set to zero through a rescaling of concentrations and therefore does not occur in our approach [56].

Appendix D. Distribution of Longer Motifs
In Section 3.3, we discussed a scenario without kinetic stalling (σ 1 = σ 2 = 1) and without energetic bias (∆γ = 0). There, the zebraness converged to Z = 0.5 in all individual realizations, pointing to an entirely random sequence pool (see Figure 3c). However, this result does not exclude a non-random distribution of longer (sub)motifs in the steady-state. To rule out correlations on larger scales, we analyze the distributions of (sub)motifs of lengths n > 2 in detail. Therefore, we compare the (sub)motif distributions P n from the simulated model dynamics to the (sub)motif distributions R n obtained from a true random process. To this end, we count the number of (sub)motifs of size n contained in all strand with L ≥ n from the simulation output at every time point. (A strand of length L ≥ n contains L − n + 1 (sub)motifs of size n.) For every time point, we also generate an equal number of random motifs of size n. We then analyze the evolution of P n and R n by means of the relative entropy D n (also called Kullback-Leibler divergence) with respect to the uniform distribution U n [144].
For the distributions P n and U n , the relative entropy D n (P n , U n ) is given by where the sum runs over all possible motifs m ∈ M n of the given length n. Using that, The relative entropy measures how much the distribution P n deviates from the distribution U n . The smaller D n (P n , U n ), the more similar are P n and U n . Since U n is the uniform distribution, i.e., the distribution with the largest entropy, we have that D n (P n , U n ) ≥ 0. Only a few long strands exist at early times and the occupation numbers for most motifs are zero. For this reason, we can not directly compare P n and R n by means of the relative entropy. The definition of relative entropy Equation (A27) requires that the relative frequency of any motif m ∈ M n is always finite for the second distribution. Therefore, D n (P n , R n ) would be ill-defined as long as are not all motifs present. Figure A2a-c show the evolution of D n (P n , U n ) (green lines) and D n (R n , U n ) (blue lines) for n = 4, 6, 8. The high values of the relative entropies at early times stem from small numbers of (sub)motifs. As the number of (sub) motifs grows with time, P n and R n become more similar to the uniform distribution, and the relative entropies decay. However, the decay of D n (P n , U n ) is slower, indicating that sequence correlations exist during early strand growth. Eventually, D n (P n , U n ) and D n (R n , U n ) converge to the same stationary value. This result implies that the motif distribution resulting from the model dynamics is entirely random. The insets in Figure A2 show the relative frequencies f n of the motifs sorted by abundance for the last time point. Simulation output and random process give similar results: Blue curves: data generated by the corresponding random process. At every time point, the number of (sub) motifs generated by the random process equals the number of (sub)motifs found in the simulation output. For small times, correlations in sequence lead to an increased relative entropy in the model dynamics. For large times, model dynamics and random processes yield similar results. The insets show the normalized frequency of (sub)motifs sorted by abundance in for the last time point.

Appendix E. Onset of Growth
In Sections 3.3-3.8, we have seen that the mean length L grows rapidly once the first new oligomers are formed from existing mononucletoides and dimers. Plotting the data with a logarithmic y-axis reveals that the growth dynamics is approximately exponential (see Figure S1). Moreover, increasing the strength of kinetic stalling (the stacking bias) shifts the onset of the rapid growth phase to later (earlier) times. The goal of this section is to derive an approximative formula for the onset of growth t.
To this end, we first focus on the scenario with ∆γ = 0 and σ 1 = σ 2 = 1 discussed in Section 3.3 (see Figure 3a). Formally, we define t as the abscissa of the intersection of the tangents to L(t) at t = 0 and the point where the growth is strongest (see dotted lines in Figure A3a).
To gain a microscopic picture of the growth dynamics, we look at the statistics of ligation events over time and distinguish between first-order and higher-order ligations. Firstorder ligations correspond to reactions where two monomers ligate to a dimer on a template of arbitrary length L ≥ 2. In contrast, higher-order ligations involve at most one monomer and lead to new strands with L ≥ 3. (We can further differentiate higher-order ligations into second and third-order ligations depending on whether the reaction comprises one or no monomer). Figure A3b reveals that t coincides approximately with the time point t high where higher-order ligations become more frequent than first-order ligations, i.e., for t = t high such that we have ν∆t are the numbers N f and N h of first-and higher-order ligations per time interval ∆t = 2.5 × 10 8 t 0 . Moreover, ν is a normalization constant. (In fact, Figure A3b reveals that t corresponds precisely to the time point, where third-order reactions become more abundant than first-order ligations).
To obtain an analytic estimate t est for t high , we consider the dynamical mean-field rate-equation for the evolution of the dimer concentration c 2 (t). This equation neglects the explicit sequence dependence and coarse-grains (de)hybridization and ligation into one effective extension. It is given bẏ The first term on the right-hand side of Equation (A31) describes the creation of a new dimer via a first-order ligation with a dimer serving as the template. k {1,1|2} is the rate constant for this effective extension process. We explain it below, together with the rate constants for the other extension processes. The second term also relates to the formation of a new dimer but using a template with L ≥ 3. The reaction is again a first-order ligation. The third term is a loss term accounting for the ligation of a monomer to a dimer, on a dimeric template. This reaction is a higher-order ligation. The fourth term is again a loss term describing a higher-order ligation of a dimer to a strand with L 1 ≥ 2 on a template with L 2 ≥ 2. The second to last term accounts for the loss of dimers due to cleavage with rate constant k cut . The last term represents the gain of dimers due to cleavage of strands of length L ≥ 3. There, a dimer can break apart at either side of the longer strand. Equation (A31) is valid if (1) concentrations are small enough such that the total strand concentration is approximately equivalent to the concentration of single strands, (2) complexes composed of more than three strands are negligible and, (3) the time scales for ligation and dehybridization are separated such that k lig k off . The assumptions (1)-(3) are satisfied (see Figure S24 and Table 2). We derive the effective rate constant for an extension k {1,1|2} as follows. First, we compute the average Boltzmann factor over the set of all complexes comprising exactly two monomers and one dimer C {1,1|2} . The Boltzmann factor associated with a specific complex C ∈ C {1,1|2} is determined by its total binding energy ∆G tot (C) (see Appendix C). The sequence-averaged Boltzmann factor then defines the sequence-averaged dissociation constant K {1,1|2} via . (A32) Recall that all concentrations are expressed as a multiple of the reference concentration c • = 1 mol/L. For that reason, the dissociation constant appears as a dimensionless quantity in Equation (A32). For complexes formed of exactly one dimer and two monomers, the total binding energy reduces to stacking energy associated with one nearest neighbor block such that This equation further reduces to Second, we multiply the inverse of the average dissociation constant with the ligation rate k lig . The result is Coarse-gaining (de)hybridization and ligation into an effective extension this way is valid as long as the ligation timescale is much slower than the dehybridization time scales such that there is enough time for the (de)hybridization dynamics to equilibrate. By our choice of the ligation rate (see Section 2.6), this premise is clearly fulfilled for monomers and dimers. The rate constant k {1,2|2} is obtained analogously. The other rate constants are less trivial. However, we will neglect them later on anyway.
For t < t high , mostly monomers and dimers populate the pool. We, therefore, assume that t high roughly matches the time point where the loss terms related to higher-order ligations balance the gain terms in Equation (A31). Equation (A31) has no analytic solution. We thus have to make (crude) approximations to obtain the estimate t est for t high . First, we neglect all terms that involve strands of length L > 2 or do not include at least one monomer. The resulting simplified equation then reads: (A36) Equation (A36) is a cubic ordinary differential equation since c 1 + 2c 2 is a constant. Hence, the simplified equation still does not allow for a closed analytic solution for c 2 (t), and further (crude) approximations will be necessary. For t → 0, the dimer concentration grows exponentially (see the dotted line in Figure A3), i.e., We now assume the monomer concentration to remain constant, i.e., c 1 (t) = c 1 (0) ≈ c tot . For t → ∞, the dimer concentration then approaches a stationary state concentration c 2 , which is given by with that, we estimate the time point for which the right-hand side of Equation (A36) vanishes by equating Equations (A37) and (A38). We obtain t est from Equation (A39) is an estimate for the time point t high where higher-order ligations become more frequent than first-order ligations. Comparing this estimate to the exact value extracted from simulation data in Figure A3b, we see that t est is only slightly smaller than t high . Hence, Equation (A39) yields a solid estimate for the transition from the first-order to the higher-order regime. Moreover, Figure A3a shows that t est matches the t with less then 10% error. We conclude that Equation (A39) is a useful approximation for the onset of the rapid growth. In addition, t est corresponds well to the time point where the dimer concentration c 2 reaches a maximum (see Figure A3c) underlining the validity of our initial assumptions. In Figure S3, we show that the prediction Equation (A39) is robust under parameter variations. We now turn to the more general case involving energetic bias ∆γ < 0 and kinetic stalling and σ 1 = σ 2 < 1. The energetic bias is automatically accounted for by averaging over all complex configurations. To include kinetic stalling, we introduce the effective sequence-averaged dissociation constant K {1,1|2} analogous to Equation (A32) as where we weight every term in the sum on the right-hand side with the overall effective stalling factor S(C) = Φ − (κ −1 (C), 1)Φ + (κ +1 (C), 1). (A41) S(C) considers both the + and the − monomer at the ligation site (see Equation (15)).  (A35)). Figures S5-S12 and S18-S23 show that the generalized approach also yields a reasonable estimate for the onset of growth.  Figure A4. A monomer hybridized adjacent to matching primer terminus (a) has a lower dehybridization rate than a monomer bound next to a non-complementary terminus (b), leading to thermodynamic stalling. Moreover, kinetic stalling reduces the ligation rate.
We now quantify the thermodynamic stalling contribution. To this end, we compare the dehybridization rates k off of the mononucleotides for the complementary and non-complementary primer termini. For simplicity, we assume an energetically unbiased scenario where ∆γ = δ = 0. The dehybridization rates depend exponentially on the hybridization energies Γ (c) and Γ (n) . Applying the procedure to compute hybridization energies described in Appendix B, Γ (c) and Γ (n) are given by The ratio of k Note that channel factors do not play a role here since they cancel out. Multiplying the kinetic and the thermodynamic stalling factors yields the combined stalling factor. For the stacking parameters used in the main text, we obtain ϑ ≈ 0.5. For σ 1 = 0.1, the combined stalling factor has a value of 0.05. This result aligns with overall stalling factors ranging between 0.1 and 0.003 measured in non-enzymatic primer extension experiments [32,38,39].
Primer extension experiments typically also consider the error fraction ω, which is defined as the ratio of the rate for an erroneous extension and the overall extension rate. We now quantify the error fraction that results from our model. To this end, we introduce the dehybridization rates k Using the same stacking parameters as above and σ 1 = 0.1, we obtain ω (n) = 0.72%. This result agrees with the experimentally observed error fraction of 0.8% in a binary DNA system containing C and G monomers [38]. The error fraction ω (n) for the extension of a primer ending with mismatch is derived analogously and reads ω (c) = 1 + 1 For σ 1 = 0.1, we find that ω (n) ≈ 3ω (c) . This observation is consistent with the experimental finding that non-complementarities at the primer terminus significantly increase the error fraction in the subsequent (stalled) extension process [32].