Toward stable replication of genomic information in pools of RNA molecules

In the VCG scenario, a circular genome is stored in a pool of oligomers, with each oligomer shorter than the genome it helps encode. Each oligomer bears a subsequence of the circular genome, such that, collectively, the oligomers encode the full genome (Figure 1A). As the spontaneous emergence of such VCG pools is expected to be rare (Chamanian and Higgs, 2022), our study focuses on the conditions under which an existing VCG pool can reliably replicate. We therefore begin with a known genome and an associated VCG pool, without addressing the question of origin. To set up our model of VCG dynamics, we specify (i) the circular genome used, (ii) the procedure by which the genome is mapped to a set of oligomers, and (iii) the chemical reactions governing the system’s evolution.

Figure 1

Download asset Open asset

We begin our analysis of the dynamics of VCG pools with an exemplary genome of length $L_{\mathrm{G}}=16\mathrm{n}\mathrm{t}$,

• 5 ′-AAAGAGGACACGGCAU-3′

• 3 ′-UUUCUCCUGUGCCGUA-5′.

This genome contains all possible monomers and dimers with equal frequency, ensuring that all motifs up to $L_{\mathrm{E}}=2\mathrm{n}\mathrm{t}$ are represented. Identifying a unique address on this genome requires at least three nucleotides. Therefore, the unique motif length is $L_{\mathrm{U}}=3\mathrm{n}\mathrm{t}$, and VCG oligomers need to be at least 3 nt long. Further below, we also explore genomes of different lengths $L_{\mathrm{G}}$, as well as varying characteristic sequence length scales $L_{\mathrm{E}}$ and $L_{\mathrm{U}}$ (genome construction detailed in the Methods section).

Based on the genome, we construct the initial oligomer pool. As a first step, we focus on a simple scenario in which the pool contains only monomers (serving as feedstock) and VCG oligomers of a single, defined length. The VCG pools are evolved in time using a stochastic simulation based on the Gillespie algorithm (Göppel et al., 2022; Laurent et al., 2024; Rosenberger et al., 2021). Since the Gillespie algorithm operates on the level of counts of molecules instead of concentrations, we must assign a volume to each system (in the range $1{\mu \mathrm{m}}^3$ to $10000{\text{μm}}^3$, see Methods). Besides the volume, we also need to choose the reaction rate constants appropriately: (i) The association time $t_0$ is the fundamental time unit in our kinetic model, and all other times are expressed relative to $t_0$. Experimentally determined association rate constants are typically around ${10}^6-{10}^7{\mathrm{M}}^{-1}{\mathrm{s}}^{-1}$(Ashwood et al., 2023; Braunlin and Bloomfield, 1991; Todisco et al., 2024a; Wetmur and Davidson, 1968). For the purpose of estimating absolute timescales, we therefore assume a constant association timescale of $t_0=(k_{\mathrm{o}\mathrm{n}}{c}^{\circ }{)}^{-1}\approx 1\text{μ}\mathrm{s}$ in the following. (ii) The timescale of dehybridization is computed via Equation 1 using the energy contribution $\gamma =-2.5k_{\mathrm{B}}T$ for a matching nearest-neighbor block and ${\gamma }_{\text{MM}}=25.0k_{\mathrm{B}}T$ in case of mismatches. The high energy penalty of nearest-neighbor blocks involving mismatches, ${\gamma }_{\text{MM}}$, is chosen to suppress the formation of mismatches, while the value of $\gamma$ roughly matches the average energy of all matching nearest-neighbor blocks given by the Turner energy model of RNA hybridization (Mathews et al., 2004). (iii) For templated ligation, we select a reaction rate constant of $k_{\mathrm{l}\mathrm{i}\mathrm{g}}={10}^{-12}{t}_0^{-1}$. This choice of $k_{\mathrm{l}\mathrm{i}\mathrm{g}}$ is consistent with ligation rates measured in enzyme-free template-directed primer extension experiments, which range from around ${10}^{-6}{\mathrm{s}}^{-1}$ $({10}^{-12}{t}_0^{-1})$ (Leveau et al., 2022; Sosson et al., 2019) to roughly ${10}^{-4}{\mathrm{s}}^{-1}$ $({10}^{-10}{t}_0^{-1})$ (Walton and Szostak, 2017), depending on the underlying activation chemistry. This indicates that, for sufficiently short oligomers (up to about 11 nt long), hybridization and dehybridization occur much faster than ligation, so binding equilibrates before ligation takes place.

Based on the ligation events observed in the simulation, we compute the observables introduced above. Due to the small ligation rate constant, it is computationally unfeasible to simulate the time evolution for more than a few ligation time units. Consequently, ligation events are scarce, which leads to high variances in the computed observables. We mitigate this issue by calculating the observables based on the concentration of complexes that are in a productive configuration, even if they do not undergo templated ligation within the time window of the simulation (Methods).

The observable $y$ (yield) quantifies the fraction of this total flux directed to producing VCG oligomers. Figure 2B shows how the yield depends on the concentration of VCG oligomers, $c_{\mathrm{V}}^{\mathrm{t}\mathrm{o}\mathrm{t}}$, at a fixed total monomer concentration, $c_{\mathrm{F}}^{\mathrm{t}\mathrm{o}\mathrm{t}}=0.1\mathrm{m}\mathrm{M}$. Here, the data points (with error bars) represent the simulation results with the different colors corresponding to different choices of initial VCG oligomer length, $L_{\mathrm{V}}$. We observe that the yield increases monotonically with the concentration of VCG oligomers and approaches 100% for high $c_{\mathrm{V}}^{\mathrm{t}\mathrm{o}\mathrm{t}}$. The concentration at which the pool reaches a yield of 50% depends on the oligomer length $L_{\mathrm{V}}$: Shorter VCG oligomers require higher concentrations for high yield. The concentration dependence of the yield can be rationalized by the types of templated ligations that are involved. For low VCG concentration, most templated ligation reactions are dimerizations (1+1) with the VCG oligomers merely acting as templates (Figure 2C). As dimers are shorter than the unique motif length $L_{\mathrm{U}}$, their formation does not contribute to the yield, which explains the low yield in the limit of small $c_{\mathrm{V}}^{\mathrm{t}\mathrm{o}\mathrm{t}}$. Conversely, for high VCG concentration, most of the templated ligations are F+V or V+V ligations, which produce oligomers of length $L\ge L_{\mathrm{U}}$, implying high yield.

Figure 2

Download asset Open asset

Figure 2C also shows that the relative contribution of V+V ligations increases with increasing $c_{\mathrm{V}}^{\mathrm{t}\mathrm{o}\mathrm{t}}$, with a large fraction of them producing incorrect products (denoted V+V,f). This reduces the fidelity of replication, $f$, and leads to a trade-off between fidelity and yield, which causes replication efficiency, $\eta =f\cdot y$, to reach a maximum at intermediate VCG concentrations (Figure 2D). Using hexamers as an example, Figure 2E illustrates why V+V ligations are prone to forming incorrect oligomers: In order to ensure that two oligomers only ligate if their sequences are adjacent to each other in the true circular genome, both oligomers need to share an overlap of at least $L_{\mathrm{U}}$ nucleotides with the template oligomer. Otherwise, the hybridization region is too short to identify the locus of the oligomer uniquely, and two oligomers from non-adjacent loci might ligate. The probability of forming incorrect products is a consequence of the combinatorics of possible subsequences in the VCG. Specifically, for a genome with $L_{\mathrm{G}}=16\text{nt}$, there are 32 different hexamers. For example, if the left educt hexamer only has one nucleotide of overlap with the template, there are $1/4\cdot 32=8$ possible educt oligomers. However, only one out of those eight hexamers is the correct partner for the right educt hexamer, implying an error probability of $1-1/8=7/8$ (first example in Figure 2E).

Characterizing replication efficiency via full simulations is computationally expensive. Depending on parameters, obtaining a single data point in Figure 2B-D can require hundreds of simulations, each lasting several days. To explore a broad parameter space more easily, we introduce an approximate adiabatic method that (i) assumes ligation is much slower than any hybridization or dehybridization event, and (ii) relies on a coarse-grained sequence-independent representation of oligomers. Details are provided in the Methods section. In brief, because ligation is rare, we first compute the equilibrium distribution of free and bound oligomers. Oligomers of the same length share a common concentration, and complex concentrations are determined via the mass action law using length-dependent dissociation constants. Combining the mass action law with a mass conservation constraint for each oligomer length allows us to compute the equilibrium concentrations of free VCG and feedstock strands, $c_{\mathrm{V}}^{\mathrm{e}\mathrm{q}}$ and $c_{\mathrm{F}}^{\mathrm{e}\mathrm{q}}$, given total concentrations $c_{\mathrm{V}}^{\mathrm{t}\mathrm{o}\mathrm{t}}$ and $c_{\mathrm{F}}^{\mathrm{t}\mathrm{o}\mathrm{t}}$. With these equilibrium values, we determine the concentrations of productive complexes and thus obtain the desired ligation‐based observables without running full stochastic simulations.

The results of the adiabatic approach agree well with the simulation data (Figure 2B-D), supporting that the replication efficiency depends non-monotonously on the concentration of VCG oligomers, with a maximum at intermediate concentration. While the available simulation data only allows for a qualitative characterization of this trend, the adiabatic approach enables a quantitative analysis. For instance, we use the adiabatic approach to determine the equilibrium concentration ratio at which replication efficiency is maximal as a function of the VCG oligomer length $L_{\mathrm{V}}$ (solid lines in Figure 2F). As expected from the qualitative trend observed in the simulation, pools containing longer oligomers reach their maximum for lower concentration of VCG oligomers. The shaded area indicates the range of VCG concentrations within which the pool’s efficiency deviates by no more than one percent from its optimum. We observe that this range of close-to-optimal VCG concentrations increases with $L_{\mathrm{V}}$. Thus, pools containing longer oligomers require less fine-tuning of the VCG concentration for replication with high efficiency.

In addition to the numerical results, we utilize the adiabatic approach to study the optimal equilibrium VCG concentration analytically (Appendix 1). We find that, for any choice of $L_{\mathrm{V}}$, replication efficiency reaches its maximum when the fractions of dimerization (1+1) reactions and erroneous V+V ligations are equal (Figure 2C for $L_{\mathrm{V}}=6\text{nt}$). This criterion can be used to derive a scaling law for the optimal equilibrium concentration ratio $c_{\mathrm{V}}^{\mathrm{e}\mathrm{q}}/c_{\mathrm{F}}^{\mathrm{e}\mathrm{q}}$ as a function of the oligomer length $L_{\mathrm{V}}$,

which is shown as dashed lines in Figure 2F (the length-scale ${\mathrm{\Lambda }}_{\mathrm{F}+\mathrm{F}}$ is defined in Appendix 1). The optimal equilibrium concentration ratio decreases exponentially with $L_{\mathrm{V}}$, while the hybridization energy $\gamma /2$ sets the inverse length scale of the exponential decay. Analytical estimate and numerical solution agree well, as long as the hybridization is weak and oligomers are sufficiently short. For strong binding and long oligomers, complexes involving more than three strands play a non-negligible role, but such complexes are neglected in the analytical approximation.

Figure 2G shows how the maximal replication efficiency depends on the VCG oligomer length. Consistent with the qualitative trend observed in Figure 2D, longer oligomers enable higher maxima in replication efficiency. Regardless of the choice of $\gamma$, replication efficiency reaches 100% if $L_{\mathrm{V}}$ is sufficiently high. Starting from Equation 2, we find the following approximation for the maximal replication efficiency attainable at a given oligomer length (dashed lines in Figure 2G),

where ${\eta }^{\circ }$ is a genome-dependent constant (Appendix 1). This equation can provide guidance for the construction of VCG pools with high replication efficiency: Given a target efficiency, the necessary oligomer length and hybridization energy, that is temperature, can be calculated. In Figure 2G, we determine the oligomer length necessary to achieve ${\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=95\mathrm{\%}$ for varying hybridization energies $\gamma$. At higher temperatures (weaker binding), VCG pools require longer oligomers to replicate with high efficiency.

Equation 3 is not restricted to the specific example genome of length $L_{\mathrm{G}}=16\mathrm{n}\mathrm{t}$, but applies more generally to genomes of arbitrary length. Any genome of length $L_{\mathrm{G}}$ can contain at most $2L_{\mathrm{G}}$ distinct motifs. Consequently, the minimum length required to specify a unique address on the genome equals $L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}=\lceil \frac{\ln (2L_{\mathrm{G}})}{\ln 4}\rceil$. By the same logic, the longest motif length for which all $4^L$ possible sequences can be exhaustively represented within the genome is $L_{\mathrm{E}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\lfloor \frac{\ln (2L_{\mathrm{G}})}{\ln 4}\rfloor$. Both of these characteristic lengths scale logarithmically with genome size. For genomes where $L_{\mathrm{E}}$ is taken to be maximal and $L_{\mathrm{U}}$ minimal, we find that the characteristic VCG oligomer length $L_{\mathrm{V}}^{\star }$ required for replication with high efficiency (${\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=95\mathrm{\%}$) also scales logarithmically with genome length (Figure 2H). Across genome lengths, the offset between $L_{\mathrm{U}}$ and $L_{\mathrm{V}}^{\star }$ is roughly constant.

In genomes with different choices of $L_{\mathrm{E}}$ and $L_{\mathrm{U}}$ $(L_{\mathrm{E}}<L_{\mathrm{E}}^{\mathrm{m}\mathrm{a}\mathrm{x}}\text{and}{L}_{\mathrm{U}}>L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}})$, the characteristic VCG oligomer length required for efficient replication, $L_{\mathrm{V}}^{\star }$, is primarily determined by the unique motif length $L_{\mathrm{U}}$. Specifically, the length of VCG oligomers, $L_{\mathrm{V}}$, must exceed $L_{\mathrm{U}}$, regardless of the value of $L_{\mathrm{E}}$. This is shown for genomes of length $L_{\mathrm{G}}=64\mathrm{n}\mathrm{t}$ in Appendix 2, where we generate genomes with specified length scales $L_{\mathrm{E}}$ and $L_{\mathrm{U}}$ using a Metropolis–Hastings algorithm and analyze their replication efficiency. Intuitively, the required VCG oligomer length $L_{\mathrm{V}}^{\star }$ is set by $L_{\mathrm{U}}$, because $L_{\mathrm{U}}$ defines the minimal hybridization region needed to ensure correct ligation via sequence-specific recognition. The precise offset between $L_{\mathrm{U}}$ and $L_{\mathrm{V}}^{\star }$ depends on genome-specific features, such as the frequency distribution of motifs with lengths between and $L_{\mathrm{U}}$. If this distribution is nearly uniform (noting that at least one motif must repeat, or else it would constitute a unique address), then $L_{\mathrm{V}}^{\star }$ will be close to $L_{\mathrm{U}}$ (Appendix 2—figure 1B). In contrast, strongly biased motif distributions require larger $L_{\mathrm{V}}^{\star }$ to achieve reliable replication (Appendix 2—figure 1A), though even in this case, $L_{\mathrm{V}}^{\star }$ typically exceeds by only a few nucleotides.

In the previous section, we characterized the behavior of pools containing VCG oligomers of a single length. We observed that V+V ligations are error-prone due to insufficient overlap between the educt strands and the template, whereas F+V ligations extend the VCG oligomers with high efficiency. The F+V ligations will gradually broaden the length distribution of the VCG pool, raising the question of how this broadening affects the replication behavior. In principle, introducing multiple oligomer lengths into the VCG pool might even improve the fidelity of V+V ligations, since a long VCG oligomer could serve as a template for the correct ligation of two shorter VCG oligomers.

To analyze this question quantitatively, we first consider the simple case of a VCG pool that only contains monomers, tetramers, and octamers. The concentration of monomers is set to $c(1)=0.1\mathrm{m}\mathrm{M}$, while the concentrations of the VCG oligomers are varied independently (Figure 3A). Replication efficiency reaches its maximum at $c(8)\approx 0.1\text{μ}\mathrm{M}$ and very low tetramer concentration, $c(4)\approx 7.4\mathrm{p}\mathrm{M}$, effectively resembling a single-length VCG pool containing only octamers. As shown in Figure 3B, the maximal efficiency is surrounded by a plateau of close-to-optimal efficiency. The octamer concentration can be varied by more than one order of magnitude without significant change in efficiency. Similarly, adding tetramers does not affect efficiency as long as the tetramer concentration does not exceed the octamer concentration. Figure 3C illustrates that the plateau of close-to-optimal efficiency coincides with the concentration regime where the ligation of a monomer to an octamer with another octamer acting as template, $\frac{1|8}{8}$, is the dominant ligation reaction. For high tetramer concentration and intermediate octamer concentration, templated ligation of tetramers on octamer templates, $\frac{4|4}{8}$, surpasses the contribution of $\frac{1|8}{8}$ ligations (green shaded area in Figure 3C). The $\frac{4|4}{8}$ reactions give rise to a ridge of increased efficiency, which, however, is small compared to the plateau of close-to-optimal efficiency (Figure 3B). Even though reactions of the type $\frac{4|4}{8}$ produce correct products in most cases, they compete with error-prone ligations like $\frac{4|8}{8}$ or $\frac{4|4}{4}$, which reduce fidelity. Increasing the binding affinity, that is lowering $\gamma$, enhances the contribution of $\frac{4|4}{8}$ ligations at the expense of $\frac{4|8}{8}$ ligations, resulting in an increased local maximum in replication efficiency (Figure 3—figure supplement 1). However, very strong hybridization energy, $\gamma =-5.0k_{\mathrm{B}}T$, is necessary for $\frac{4|4}{8}$ ligations to give rise to similar efficiency as $\frac{1|8}{8}$ ligations. Thus, in this example, high efficiency replication facilitated by the ligation of short VCG oligomers on long VCG oligomers is in theory possible, but requires unrealistically strong binding affinity. Moreover, this mechanism is only effective if educts and template differ significantly with respect to their length. If the involved VCG oligomers are too similar in size, templated ligation of two short oligomers tends to be error-prone, irrespective of the choice of binding affinity (Figure 3—figure supplement 2 and Figure 3—figure supplement 3). In that case, F+V ligations remain the most reliable replication mechanism.

Figure 3 with 3 supplements see all

Download asset Open asset

We observe similar behavior in VCG pools containing a range of oligomer lengths from $L_{\mathrm{V}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ to $L_{\mathrm{V}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$, rather than just tetramers and octamers. For a uniform VCG concentration profile (Figure 4A, dark gray), the maximal replication efficiency (at the optimal VCG concentration) is attained when F+V ligations involving long VCG oligomers dominate the templated ligation (see Figure 4D showing the contribution of different ligation types in pools with varying $L_{\mathrm{V}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ at fixed $L_{\mathrm{V}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=10\mathrm{n}\mathrm{t}$). Importantly, the maximal replication efficiency is always bounded by the efficiency of the longest VCG oligomer in the pool, independent of the presence or length of shorter oligomers (blue curve in Figure 4B; identical to the orange curve in Figure 2G). Including short VCG oligomers has minimal effect on the dominant ligation types and only slightly increases the proportion of unproductive F+F ligations. Consequently, reducing $L_{\mathrm{V}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ while keeping $L_{\mathrm{V}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$ fixed leads to only a modest decline in maximal efficiency. In contrast, decreasing $L_{\mathrm{V}}^{\max }$ while holding $L_{\mathrm{V}}^{\min }$ constant causes a substantial reduction in efficiency (Figure 4B), because the longest oligomer in the pool sets an upper bound on replication efficiency. Moreover, at low $L_{\mathrm{V}}^{\max }$, short and long oligomers become more similar in length, giving rise to a spectrum of erroneous V+V ligations that compete with the productive F+V ligations (Figure 4E).

Figure 4

Download asset Open asset

In a realistic prebiotic scenario, the concentration profile of the VCG pool would not be uniform. Depending on the mechanism producing the pool and its coupling to the non-equilibrium environment, it might have a concentration profile that decreases or increases exponentially with length. We use the parameter ${\kappa }_{\mathrm{V}}$ to control this exponential length dependence (Figure 4A): For negative ${\kappa }_{\mathrm{V}}$, the concentration increases as a function of length, while exponentially decaying length distributions have positive ${\kappa }_{\mathrm{V}}$. We find that replication efficiency is high if the concentration of long VCG oligomers exceeds or at least matches the concentration of short VCG oligomers (${\kappa }_{\mathrm{V}}\le 0$ in Figure 4C). In that case, replication efficiency is dominated by the long oligomers in the pool, since these form the most stable complexes. As the concentration of long oligomers is decreased further (${\kappa }_{\mathrm{V}}>0$), the higher stability of complexes formed by longer oligomers is eventually insufficient to compensate for the reduced concentration of long oligomers. Replication efficiency is then governed by short VCG oligomers, which exhibit lower replication efficiency. In the limit ${\kappa }_{\mathrm{V}}\to \mathrm{\infty }$, replication efficiency approaches the replication efficiency of a single-length VCG pool containing only oligomers of length $L_{\mathrm{V}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ (Figure 4C and Figure 2G).

So far, we focused on ensembles that contain solely monomers as feedstock. However, examining the influence of dimers on replication in VCG pools is of interest, since dinucleotides have proven to be interesting candidates for enzyme-free RNA copying (Leveau et al., 2022; Sosson et al., 2019; Walton and Szostak, 2016). For this reason, we study oligomer pools like those illustrated in Figure 5A: The ensemble contains monomers, dimers, and VCG oligomers of a single length, $L_{\mathrm{V}}$. As our default scenario, the dimer concentration is set to 10% of the monomer concentration, corresponding to ${\kappa }_{\mathrm{F}}=2.3$, but this ratio can be modified by changing ${\kappa }_{\mathrm{F}}$.

Figure 5

Download asset Open asset

Figure 5B compares the replication efficiency of a pool with ($L_{\mathrm{F}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=2\text{nt}$) and without dimers ($L_{\mathrm{F}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=1\text{nt}$). In both cases, pools that are rich in VCG oligomers exhibit low efficiency of replication. As erroneous V+V ligations are the dominant type of reaction in this limit, all pools achieve the same efficiency regardless of the presence of dimers. In contrast, when the pool is rich in feedstock (small $c_{\mathrm{V}}^{\mathrm{t}\mathrm{o}\mathrm{t}}/c_{\mathrm{F}}^{\mathrm{t}\mathrm{o}\mathrm{t}}$), pools with and without dimers behave differently: If only monomers are included, the efficiency approaches zero, as dimerizing monomers do not contribute to the yield, and thus not to the efficiency. However, the presence of dimers enables the ligation of monomers and dimers to form trimers and tetramers, which lead to a non-zero yield. Given the low VCG concentration, the ligations are likely to proceed using dimers as template. As a consequence, educt oligomers can only hybridize to the template with a single-nucleotide-long hybridization region, leading to frequent formation of incorrect products and, consequently, low efficiency.

Replication efficiency reaches its maximum at intermediate VCG concentrations, where replication is dominated by F+V ligations. Notably, the maximal attainable efficiency is significantly lower for pools with dimers than without, as dimers increase the number and the stability of complex configurations that can form incorrect products (Figure 5C). Without dimers, ligation products are only incorrect if the overlap between the VCG educt oligomer and template is shorter than the unique motif length, $L_{\mathrm{U}}$. With dimers, however, dangling end dimers can cause incorrect products even in case of long overlap of educt oligomer and template (right column in Figure 5D). The stability of the latter complexes depends on the length of the VCG oligomers, $L_{\mathrm{V}}$, whereas the stability of complexes facilitating incorrect monomer addition is independent of oligomer length (Figure 5D).

In the presence of dimers, the length-dependent stability of complexes allowing for correct and incorrect F+V ligations causes a competition, which sets an upper bound on the efficiency of replication (Appendix 3),

Here, we introduced effective association constants ${\mathcal{K}}^{\mathrm{a}}$, which depend differently on the VCG oligomer length, $L_{\mathrm{V}}$. While the effective association constant ${\mathcal{K}}_{1+\mathrm{V},\mathrm{f}}^{\mathrm{a}}$ of complexes enabling incorrect 1+V ligations is length-independent, the effective association constant for incorrect 2+V ligations, ${\mathcal{K}}_{2+\mathrm{V},\mathrm{f}}^{\mathrm{a}}$, scales exponentially with $L_{\mathrm{V}}$,

The effective association constants for correct 1+V and 2+V ligations also scale exponentially with the oligomer length (Appendix 3—figure 1).

In systems without dimers, that is ${\kappa }_{\text{F}}\to \mathrm{\infty }$, ${\eta }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ approaches 100%, which is consistent with the behavior observed in the previous sections. Conversely, in systems containing dimers, the maximal efficiency remains at a value below 100%, which depends on the concentration of dimers in the feedstock. Figure 5C shows the dependence of maximal replication efficiency on the length of VCG oligomers in pools containing monomers, dimers, and VCG oligomers. As $L_{\mathrm{V}}$ increases, ${\eta }_{\text{max}}$ converges towards the upper bound defined in Equation 4 (dashed line in Figure 5C).

In all scenarios considered so far, the efficiency of replication is limited by a common mechanism, regardless of the specifics of the VCG pool: Ligations involving two oligomers that hybridize to the template over a region shorter than $L_{\mathrm{U}}$ are prone to generate incorrect products. In previous sections, we minimized these erroneous ligations by fine-tuning the concentration and length of VCG oligomers. However, such control may become unnecessary if the typically error-prone ligation of two oligomers (V+V) is kinetically suppressed. Kinetic suppression can be an intrinsic property of the activation chemistry: Templated ligation of two oligomers can be several orders of magnitude slower than the extension of an oligomer by a single monomer (Ding et al., 2022; Prywes et al., 2016). As a result, V+V ligations are disfavored purely by their slower kinetics. In addition, it is conceivable that only monomers are chemically activated while longer oligomers remain inactive, which would further reduce the likelihood of erroneous ligations. This scenario has already been explored experimentally (Ding et al., 2023). In natural environments, it could occur, for instance, when activated monomers are produced externally and then diffuse into compartments containing the VCG but lacking internal activation pathways (Kriebisch et al., 2024; Toparlak et al., 2023).

Within our model, we capture the kinetic suppression by introducing two different rates of ligation, $k_{\text{lig,1}}$ for ligations involving a monomer and $k_{\text{lig,>1}}$ for ligations involving no monomer, allowing for kinetic discrimination between these processes. We explore the resulting replication efficiency in the limit of perfect kinetic discrimination ($k_{\text{lig,>1}}/k_{\text{lig,1}}\to 0$) where only monomers are reactive for ligation. We first consider a pool where the reactive monomers are mixed with VCG oligomers of a single length as well as non-reactive dimers (Figure 6A). We vary the concentration of VCG oligomers, but keep the feedstock concentrations constant. For small VCG concentrations, we observe low efficiencies (Figure 6B), as the ligation of two monomers, or one monomer and one dimer, is most likely. Conversely, high $c_{\mathrm{V}}^{\mathrm{t}\mathrm{o}\mathrm{t}}$ facilitates the formation of complexes in which VCG oligomers are extended by monomers, which implies high efficiency (Figure 6B). Note that, unlike in systems where all oligomers are reactive, replication efficiency does not decrease for high VCG concentration, as erroneous V+V ligations are impossible. Instead, perfect replication efficiency (100%) is reached for sufficiently high $L_{\mathrm{V}}$.

Figure 6

Download asset Open asset

While replication efficiency characterizes the relative amount of nucleotides used for the correct elongation of VCG oligomers, it is also interesting to analyze which fraction of VCG oligomers is in a monomer-extension-competent state,

Here, $c_{1+\mathrm{V}}^{\mathrm{e}\mathrm{q}}$ denotes the equilibrium concentration of all complexes enabling the addition of a monomer to any VCG oligomer. We find that $r_{1+\mathrm{V}}$ depends on the VCG concentration qualitatively in the same way as the efficiency: $r_{1+\mathrm{V}}$ is small for small VCG concentration but approaches a value $r_{1+\mathrm{V}}^{\mathrm{\infty }}$ asymptotically for high VCG concentration (Figure 6C). In this limit, almost all VCG oligomers are part of a duplex. Monomers can bind to these duplexes to form complexes allowing for 1+V ligations. Since almost all VCG oligomers are already part of a duplex, increasing $c_{\mathrm{V}}^{\mathrm{t}\mathrm{o}\mathrm{t}}$ further does not increase the fraction of VCG oligomers that can be extended by monomers. Instead, the asymptotic value is determined by the concentration of monomers and their binding affinity $K_{\mathrm{d}}(1)$ to an existing duplex (Appendix 4),

While the asymptotic value $r_{1+\mathrm{V}}^{\mathrm{\infty }}$ does not depend on the length of the VCG oligomers, the threshold VCG concentration at which $r_{1+\mathrm{V}}=r_{1+\mathrm{V}}^{\mathrm{\infty }}/2$ scales exponentially with $L_{\mathrm{V}}$ (Appendix 4),

This scaling implies that longer oligomers require exponentially lower VCG concentration to achieve a given ratio $r_{1+\mathrm{V}}$ (Figure 6C), as their greater length allows them to form more stable complexes. This observation implies that pools with longer oligomers will always be more productive than pools with shorter oligomers (at equal VCG concentration).

The behavior becomes more complex in pools containing VCG oligomers of multiple lengths, due to the competitive binding within such heterogeneous pools. To illustrate this, we examine an ensemble containing VCG oligomers ranging from $L_{\mathrm{V}}^{\mathrm{m}\mathrm{i}\mathrm{n}}=3\text{ nt}$ to $L_{\mathrm{V}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=9\text{ nt}$, along with the same feedstock as previously (reactive monomers and non-reactive dimers, see Figure 7A). For simplicity, we assume that the length distribution of VCG oligomers is uniform. We study the fraction of oligomers in a monomer-extension-competent state as a function of oligomer length,

Figure 7

Download asset Open asset

where $c_{1+V}^{\mathrm{e}\mathrm{q}}(L)$ denotes the equilibrium concentration of all complexes that allow monomer-extension of VCG oligomers of length $L$. At low VCG concentration, longer oligomers are more likely to be extended by monomers than shorter ones (Figure 7B). This behavior is intuitive, as longer oligomers tend to form more stable complexes, which lead to higher productivity. Surprisingly, increasing the VCG concentration reverses the length-dependence of the productivity, such that short oligomers are more likely to be extended by monomers than long ones (note the three crossings of the curves in Figure 7B). For example, 8-mers are more likely to undergo primer extension than 9-mers once the VCG concentration exceeds $\approx 0.5\text{μ}\mathrm{M}$. We derived a semi-analytical expression for the threshold VCG concentrations at which oligomers of two different lengths have equal productivity (dashed lines in Figure 7B, Appendix 5 and Appendix 6).

To understand the mechanism underlying this inversion of productivity, we analyze how different complex types contribute to $r_{1+\mathrm{V}}(L)$. We introduce

which denotes the fraction of oligomers of length $L_{\mathrm{E}}$ that are in a monomer-extension-competent complex configuration that uses an oligomer of length $L_{\mathrm{T}}$ as template. The term $c_{1+\mathrm{V}}^{\mathrm{e}\mathrm{q}}\left(\frac{1|L_{\mathrm{E}}}{L_{\mathrm{T}}}\right)$ includes the sum over all possible configurations of complexes with the given lengths. Focusing on 8-mers and 9-mers as an example, we observe that complexes utilizing the 9-mer as template are responsible for the inversion of productivity (top two curves in Figure 7C): As the VCG concentration increases, the fraction of monomer-extendable 8-mers eventually surpasses the fraction of monomer-extendable 9-mers, that is $r_{1+\mathrm{V}}\left(\frac{1|8}{9}\right)>r_{1+\mathrm{V}}\left(\frac{1|9}{9}\right)$ (Figure 7C). Two factors give rise to this feature: (i) Ternary complexes of type $\frac{1|8}{9}$ and $\frac{1|9}{9}$ have similar stability (Figure 7E), and (ii) the equilibrium concentration of free 8-mers is higher than that of 9-mers (Figure 7D). As a result, 8-mers are more likely than 9-mers to hybridize to an existing duplex $\frac{1}{9}$, and, given the stability of the complex $\frac{1|8}{9}$, 8-mers remain bound almost as stably as 9-mers. In summary, short oligomers sequester long oligomers as templates to enhance their monomer-extension rate, while long oligomers cannot make use of short oligomers as templates due to the relative instability of the corresponding complexes.

It is noteworthy that Ding et al. have already observed productivity inversion experimentally (Ding et al., 2023). In their study, they included activated monomers, activated imidazolium-bridged dinucleotides and oligomers up to a length of 12 nt and observed that the primer extension rate for short primers is higher than the extension rate of long primers. Even though our model differs from their setup in some aspects (e.g. different circular genome, no bridged dinucleotides), evaluating our model using parameters similar to those of the experimentally studied system predicts productivity inversion that qualitatively agrees with the experimental findings (Appendix 7). We therefore assume that the mechanism underlying productivity inversion described here also applies to the experimental observations.

In the Virtual Circular Genome (VCG) scenario, genomes are encoded in a pool of oligomers. The encoded genomes are assumed to be circular sequences of length $L_{\mathrm{G}}$, containing both the original sequence and its reverse complement. Each genome is characterized by two fundamental length scales that reflect different aspects of motif distribution along the sequence. The minimal unique motif length, $L_{\mathrm{U}}$, is defined as the shortest subsequence length for which all motifs of length $L\ge L_{\mathrm{U}}$ appear at most once in the genome. In contrast, the exhaustive coverage length, $L_{\mathrm{E}}$, denotes the largest motif length for which all $4^{L_{\mathrm{E}}}$ possible motifs are present within the genome. Since only $2L_{\mathrm{G}}$ distinct motifs can be encoded in a genome (including its complement), $L_{\mathrm{E}}$ cannot exceed

Similarly, for a motif to be uniquely addressable on the genome, its length must be at least

We note that $L_{\mathrm{E}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$ and $L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$ are essentially the same length scale (differing by at most one nucleotide).

The characteristic length scales $L_{\mathrm{E}}$ and $L_{\mathrm{U}}$ impose constraints on how motifs are distributed. For example, when $L=L_{\mathrm{U}}$, all motifs of length $L$ must appear at most once, while at least one motif of length $L-1$ must occur more than once. To quantify the motif distribution, we introduce the motif entropy,

where $f_i$ denotes the relative frequency of motif $i$ across the genome and its reverse complement. Motif entropy ranges from zero (a homogeneous sequence with only one motif) to a maximum value that depends on the subsequence length,

For motif length $L$ to qualify as the unique motif length $L_{\mathrm{U}}$, its entropy must be maximal, $S(L)=S^{\mathrm{m}\mathrm{a}\mathrm{x}}(L)$, while $S(L-1)$ must be smaller than its respective maximum, $S(L-1)<S^{\mathrm{m}\mathrm{a}\mathrm{x}}(L-1)$.

The correspondence between characteristic length scales and motif entropies provides a way to construct genome sequences with specified motif characteristics. By treating the entropy function as an effective ‘Hamiltonian’ $\mathcal{H}$, we can generate genome sequences through Metropolis–Hastings sampling. In our implementation, each update step in the Metropolis-Hastings algorithm involves either a single-nucleotide mutation or a cut-and-paste operation that relocates a segment of the genome to a new position (the cut-and-paste operation is 10 times more likely than the single nucleotide mutation). The acceptance criterion follows the standard Metropolis rule: modifications that reduce the Hamiltonian are always accepted, while increases in ‘energy’ are accepted with probability $\exp \left[-\beta (E_{\mathrm{o}\mathrm{l}\mathrm{d}}-E_{\mathrm{n}\mathrm{e}\mathrm{w}})\right]$. Here, ${\beta }^{-1}$ is an effective temperature chosen to be small compared to the typical energy to ensure convergence to the minimum. Simulations are either run until a predefined entropy target is reached or until the energy converges to a plateau.

To generate genomes with $L_{\mathrm{E}}=L_{\mathrm{E}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$ and $L_{\mathrm{U}}=L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$, we minimize the Hamiltonian

Starting from a random sequence, we perform 10,000 Metropolis–Hastings steps at an inverse temperature $\beta ={10}^{-5}$ to construct genome sequences of lengths $L_{\mathrm{G}}=16\text{nt}$ (listed in Results) and $L_{\mathrm{G}}=64\text{nt}$ (listed in Supplementary file 1). To explore genomes with $L_{\mathrm{E}}<L_{\mathrm{E}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$ and $L_{\mathrm{U}}>L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$, we initialize the algorithm from the maximally diverse genomes ($L_{\mathrm{E}}=L_{\mathrm{E}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$, $L_{\mathrm{U}}=L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$) and then reduce entropy across the range $L_{\mathrm{E}}<L<L_{\mathrm{U}}$. This is done by minimizing the Hamiltonian

via two different sampling protocols. In the first protocol, the simulation is terminated as soon as the genome reaches the desired values of $L_{\mathrm{E}}$ and $L_{\mathrm{U}}$. The resulting motif distributions on the intermediate length scales ($L_{\mathrm{E}}<L<L_{\mathrm{U}}$) remain close to uniform, with only minor biases sufficient to enforce the length-scale constraints. In the second protocol, entropy minimization continues beyond the point at which the target values are achieved, leading to more strongly biased motif distributions on intermediate length scales. These construction strategies allow us to systematically tune genome complexity and motif structure, enabling controlled investigations of how the characteristic length scales influence replication dynamics (see Appendix 2 for details).

We simulate the dynamics of VCG pools using a kinetic simulation that is based on the Gillespie algorithm. In the simulation, oligomers can hybridize to each other to form complexes or dehybridize from an existing complex. Moreover, two oligomers can undergo templated ligation if they are hybridized adjacent to each other on a third oligomer. At each time $t$, the state of the system is determined by a list of all single-stranded oligomers and complexes as well as their respective copy number. We refer to the state of the system at the time $t$ as the ensemble of compounds ${\mathcal{E}}_t$. Given the copy numbers, the rates $r_i$ of all possible chemical reactions $i\in \mathcal{I}$ can be computed. To evolve the system in time, we need to perform two steps: (i) We sample the waiting time until the next reaction, $\tau$, from an exponential distribution with mean $(\sum _{i\in \mathcal{I}}{r}_i{)}^{-1}$, and update the simulation time, $t\to t+\tau$. (ii) We pick which reaction to perform by sampling from a categorical distribution. Here, the probability to pick reaction $i$ equals $r_i/(\sum _{i\in \mathcal{I}}{r}_i)$. The copy numbers are updated according to the sampled reaction, yielding ${\mathcal{E}}_{t+\tau }$. Steps (i) and (ii) are repeated until the simulation time $t$ reaches the desired final time, $t_{\text{final}}$. A more detailed explanation of the kinetic simulation is presented in Göppel et al., 2022; Rosenberger et al., 2021.

Our goal is to compute observables characterizing replication in the VCG scenario based on the full kinetic simulation. In the following derivation, we focus on one particular observable (yield) for clarity. The results for other observables are stated directly, as their derivations follow analogously. Recall the definition of the yield introduced in the Results section,

As we are interested in the initial replication performance of the VCG, we compute the yield based on the ligation events that take place until the characteristic timescale of ligations ${\tau }_{\text{lig}}=k_{\text{lig}}^{-1}\approx {10}^{12}{t}_0$. In principle, we would like to compute the yield based on the templated ligation events that we observe in the simulation. Unfortunately, for reasonable system parameters, it is impossible to simulate the system long enough to observe sufficiently many ligation events to compute $y$ to reasonable accuracy. For example, for a VCG pool containing monomers at a total concentration of $c_{\mathrm{F}}^{\text{tot }}=0.1\mathrm{m}\mathrm{M}$ and VCG oligomers of length $L=8\mathrm{n}\mathrm{t}$ at a total concentration of $c_{\mathrm{V}}^{\mathrm{t}\mathrm{o}\mathrm{t}}=1\text{μ}\mathrm{M}$, it would take about 1700 hr of simulation time to reach $t=5\cdot {10}^{12}{t}_0$ (Figure 8). Multiple such runs would be needed to estimate the mean and the variance of the observables of interest, rendering this approach unfeasible.

Figure 8

Download asset Open asset

Instead, we compute the replication observables based on the copy number of complexes that could potentially perform a templated ligation, that is complexes in which two strands are hybridized adjacent to each other, such that they could form a covalent bond. We can show analytically that the number of productive complexes is a good approximation for the number of incorporated nucleotides: The number of incorporated nucleotides can be computed as the integral over the ligation flux, weighted by the number of nucleotides that are added in each templated ligation reaction,

Here, $N(\text{C})$ denotes the copy number of the complex C in the pool ${\mathcal{E}}_t$. $L_{\text{e,1}}$ and $L_{\text{e,2}}$ denote the lengths of the oligomers that undergo ligation, and $\mathbb{1}$ is an indicator function which enforces that only complexes in a ligation-competent configuration contribute to the reaction flux. As only a few ligation events are expected to happen until ${\tau }_{\text{lig}}$, it is reasonable to assume that the ensembles ${\mathcal{E}}_t$ do not change significantly during $t\in [0,{\tau }_{\text{lig}}]$. Therefore, the integration over time may be interpreted as a multiplication by ${\tau }_{\text{lig}}$,

where $⟨\dots ⟩$ denotes the average over realizations of the ensembles ${\mathcal{E}}_t$ within the time interval $t\in [{\tau }_{\text{eq}},{\tau }_{\text{lig}}]$. This average corresponds to the average number of complexes in a ligation-competent configuration. Note that, at this point, we made the additional assumption that no templated ligations are taking place between $[0,{\tau }_{\text{eq}}]$. This assumption is reasonable, as (i) the equilibration process is very short compared to the characteristic timescale of ligation, and (ii) the number of complexes that might allow for templated ligation during equilibration is lower than in equilibrium (we start the simulation with an ensemble of single-stranded oligomers). Both aspects imply that the rate of templated ligation is negligible during the interval $[0,{\tau }_{\mathrm{e}\mathrm{q}}]$.

In order to compute the average over different realizations of ensembles $\mathcal{E}$ (as required in Equation 6), we need to sample a set of uncorrelated ensembles that have reached the hybridization equilibrium, which can be done using the full kinetic simulation. The simulation starts with a pool containing only single-stranded oligomers and reaches the (de)hybridization equilibrium after a time ${\tau }_{\text{eq}}$. We identify this timescale of equilibration by fitting an exponential function to the total hybridization energy of all complexes in the system, $\mathrm{\Delta }{G}_{\text{tot}}$ (Figure 9A). In the set of ensembles used to evaluate the average in Equation 6, we only include ensembles for time $t>{\tau }_{\text{eq}}$ to ensure that the ensembles have reached (de)hybridization equilibrium. To ensure that the ensembles are uncorrelated, we require that the time between two ensembles that contribute to the average is at least ${\tau }_{\text{corr}}$. The correlation time, ${\tau }_{\text{corr}}$, is determined via an exponential fit to the autocorrelation function of $\mathrm{\Delta }{G}_{\text{tot}}$ (Figure 9B). Besides computing the expectation value (Equation 6), we are also interested in the ‘uncertainty’ of this expectation value, that is in the standard deviation of the sample mean ${\sigma }_{⟨X⟩}$. (We use $X$ as a short-hand notation for $\sum _{\text{C}\in \mathcal{E}}N(\text{C})\text{min}(L_{\text{e,1}},L_{\text{e,2}})\mathbb{1}(\mathrm{C}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{s}\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n})$). The standard deviation of the sample mean, ${\sigma }_{⟨X⟩}$, is related to the standard deviation of $X$, ${\sigma }_X$, by the number of samples, ${\sigma }_{⟨X⟩}=(N_s{)}^{-1/2}{\sigma }_X$. Moreover, based on the van-Kampen system size expansion, we expect the standard deviation of $X$ to be proportional to $V^{-1/2}$, such that ${\sigma }_{⟨X⟩}\propto (N_{s}V{)}^{-1/2}$.

Figure 9

Download asset Open asset

Using Equation 6 (as well as an analogous expression for the number of nucleotides that are incorporated in VCG oligomers), the yield can be expressed as

The additional condition $\mathbb{1}(L_{\text{e,1}}+L_{\text{e,2}}\ge L_{\mathrm{U}})$ in the numerator ensures that the product oligomer is long enough to be counted as a VCG oligomer, that is at least $L_{\mathrm{U}}$ nucleotides long. Analogously, the expression for the fidelity of replication reads

Multiplying fidelity and yield results in the efficiency of replication,

The ligation share of a particular type of templated ligation $s(\text{type})$, that is, the relative contribution of this templated-ligation type to the nucleotide extension flux, can be represented in a similar form as the other observables,

As all observables are expressed as the ratio of two expectation values, $\mathcal{Z}=⟨X⟩/⟨Y⟩$, we can compute the uncertainty of the observables via Gaussian error propagation,

Since the variances, ${\sigma }_{⟨X⟩}^2$ and ${\sigma }_{⟨Y⟩}^2$, as well as the covariance, ${\sigma }_{⟨X⟩,⟨Y⟩}^2$, are proportional to $(N_{s}V{)}^{-1}$, the standard deviation of the observable mean, ${\sigma }_Z$, scales with the inverse square root of the number of samples and the system volume, that is ${\sigma }_Z\propto (N_{s}V{)}^{-1/2}$. Therefore, the variance of the computed observable can be reduced by either increasing the system volume or increasing the number of samples used for averaging. Both approaches incur the same computational cost: (i) Increasing the number of samples, $N_s$, requires running the simulation for a longer duration, with the additional runtime scaling linearly with the number of samples. (ii) Similarly, the additional runtime needed due to increased system volume, $V$, also scales linearly with $V$ (Figure 8). One update step in the simulation always takes roughly the same amount of runtime, but the change in simulation time per update step depends on the total rate of all reactions in the system. The total rate is dominated by the association reactions, and their rate is proportional to the volume. Therefore, the change in simulation time per update step is proportional to $V^{-1}$. The runtime, which is necessary to reach the same simulation time in a system with volume $V$ as in a system with volume 1, is a factor of $V$ longer in the larger system. With this in mind, it makes no difference whether the variance is reduced by increasing the volume or the number of samples. For practical reasons (post-processing of the simulations is less memory- and time-consuming), we opt to choose a moderate number of samples, but slightly higher system volumes to compute the observables of interest. The simulation parameters (length of oligomers, concentrations, hybridization energy, volume, number of samples, characteristic timescales) used to obtain the results presented in Figure 2 are summarized in Table 1.

To characterize the replication performance of VCG pools across a broad range of system parameters, we developed an adiabatic approach that enables faster computation of replication observables than full kinetic simulations. In this method, we compute the equilibrium concentrations of complexes in productive configurations under the assumption that templated ligation is much slower than (de)hybridization. The approach relies on a coarse-grained representation of the oligomers in the pool, which we introduce in this section.

Duplexes

Two strands can form a duplex by hybridizing to each other. We refer to the bottom oligomer as ‘Crick’ strand C and to the top oligomer as ‘Watson’ strand W. A duplex is uniquely characterized by the lengths of the oligomers, $L_{\text{C}}$ and $L_{\text{W}}$, as well as their relative alignment (Figure 10A). The alignment index $i$ denotes the position of the Watson strand with respect to the Crick strand. As there needs to be at least one nucleotide of overlap between the strands for a duplex to exist, the alignment index needs to be in the interval $i\in [-(L_{\text{W}}-1),L_{\text{C}}-1]$. Using the alignment index, we can also determine if the duplex has a left (or right) dangling end. The corresponding indicator variables are called $d_l$ (or $d_r$),

${\displaystyle {\displaystyle {\displaystyle \begin{array}{l}{d}_l=\{\begin{array}{l}1\text{if}i=0,\\ 0\text{otherwise},\end{array}{d}_r=\{\begin{array}{l}1\text{if}i+L_{\text{W}}=L_{\text{C}},\\ 0\text{otherwise}.\end{array}\end{array}}}}$

Figure 10

Download asset Open asset

Moreover, the length of the hybridization region $L_{\text{o}}$ can be computed via

${\displaystyle {\displaystyle {\displaystyle L_{\text{o}}=\text{min}(L_{\text{C}},i+L_{\text{W}})-\text{max}(i,0).}}}$

The hybridization energy of the duplex depends on the length of the hybridization region as well as on the existence/absence of dangling ends. For a hybridization site of length $L_{\text{o}}$, there are $L_{\text{o}}-1$ nearest-neighbor energy blocks each of which contributes $\gamma$ to the energy. Moreover, each dangling end contributes $\frac{\gamma }{2}$ to the energy,

${\displaystyle {\displaystyle {\displaystyle E=\gamma (L_{\text{o}}-1)+\frac{\gamma }{2}(d_l+d_r).}}}$

To compute the combinatorial multiplicity for a duplex with fixed $L_{\text{C}}$, $L_{\text{W}}$ and alignment index $i$, we need to multiply the combinatorial multiplicity of the Crick strand by the number of possible hybridization partners. We assume that a hybridization partner is possible if its sequence is perfectly complementary to the lower strand within the hybridization region, whereas hybridization partners with mismatches are not accounted for. This is sensible as long as the energetic penalty for mismatches in the full kinetic simulation is sufficiently large to suppress mismatches. The number of possible hybridization partners is determined by the length of the overlap region $L_{\text{o}}$: If $L_{\text{o}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}$, the pool contains only one oligomer sequence that can act as hybridization partner by construction of the genome. For shorter hybridization regions, multiple hybridization partners might be possible. Their number is set by the combinatorial multiplicity of the Watson oligomer divided by the combinatorial multiplicity of the hybridization region,

${\displaystyle {\displaystyle {\displaystyle \begin{array}{l}\mathfrak{C}(L_{\text{C}},L_{\text{W}},i)=\{\begin{array}{ll}\frac{\mathfrak{C}(L_{\text{C}})\mathfrak{C}(L_{\text{W}})}{4^{L_{\text{o}}}}& \text{if}{L}_{\text{o}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}},\\ \mathfrak{C}(L_{\text{C}})& \text{if}{L}_{\text{o}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}.\end{array}\end{array}}}}$

To avoid double-counting, we only account for complexes in which the Crick strand is at least as long as the Watson strand, $L_{\text{C}}\ge L_{\text{W}}$, and multiply $\mathfrak{C}(L_{\text{W}},L_{\text{C}},i)$ by $1/2$ if $L_{\text{W}}=L_{\text{C}}$.

Quaternary complexes

The largest complexes to be accounted for in our coarse-grained adiabatic approach are quaternary complexes, that is complexes comprised of four strands. We need to distinguish three types of such complexes: (i) 3-1 quaternary complexes, (ii) left-tilted 2-2 quaternary complexes, and (iii) right-tilted 2-2 quaternary complexes. In 3-1 quaternary complexes, three Watson strands are hybridized to one Crick strand (Figure 11), whereas in 2-2 quaternary complexes, two Watson strands are hybridized to two Crick strands (Figure 12 and Figure 13).

Figure 11

Download asset Open asset

Figure 12

Download asset Open asset

Figure 13

Download asset Open asset

3-1 Quaternary complexes

Figure 11 depicts a typical 3-1 quaternary complex. Such a complex is uniquely characterized by the length of its oligomers, $L_{\text{C}},L_{\text{W,1}},L_{\text{W,2}},L_{\text{W,3}}$, as well as their relative position to each other denoted by the alignment indices $i,j$, and $k$. All positions within the ternary complex are measured relative to the left end of the C strand. Any W strand needs to have at least one nucleotide of overlap with the C strand, but two W strands must never occupy the same position on the C strand. Consequently, the alignment indices fall within the intervals,

${\displaystyle {\displaystyle {\displaystyle i\in [-(L_{\text{W,1}}-1),L_{\text{C}}-L_{\text{W,1}}-L_{\text{W,2}}-1],j\in [i+L_{\text{W,1}},L_{\text{C}}-L_{\text{W,2}}-1],\text{and}k\in [j+L_{\text{W,2}},L_{\text{C}}-1].}}}$

There are two dangling ends (left and right) and potentially two gaps between the W strands: one gap between W1 and W2 and another one between W2 and W3. The following boolean variables indicate the presence/absence of the respective dangling ends,

${\displaystyle {\displaystyle \begin{array}{rlrl}& d_l=\{\begin{array}{ll}0& \text{if}i=0,\\ 1& \text{otherwise,}\end{array}& & d_{m1}=\{\begin{array}{ll}0& \text{if}i+L_{\text{W,1}}=j,\\ 1& \text{otherwise},\end{array}\\ & d_{m2}=\{\begin{array}{ll}0& \text{if}i+L_{\text{W,2}}=k,\\ 1& \text{otherwise},\end{array}& & d_r=\{\begin{array}{ll}0& \text{if}k+L_{\text{W,3}}=L_{\text{C}},\\ 1& \text{otherwise}.\end{array}\end{array}}}$

The length of the hybridization regions is given by

${\displaystyle {\displaystyle {\displaystyle L_{\text{o,1}}=i+L_{\text{W,1}}-\text{max}(i,0),L_{\text{o,2}}=L_{\text{W,2}},\text{and}{L}_{\text{o,3}}=\text{min}(L_{\text{C}},k+L_{\text{W,3}})-k.}}}$

Following the same reasoning as in the case of ternary complexes, the energy equals

${\displaystyle {\displaystyle {\displaystyle E=\gamma (L_{\text{o,1}}+L_{\text{o,2}}+L_{\text{o,3}}-3)+\gamma (1-d_{m1})+\gamma (1-d_{m2})+\frac{\gamma }{2}(d_l+d_r+2d_{m1}+2d_{m2}).}}}$

Similarly, the combinatorial multiplicity of 3-1 quaternary complexes is constructed using the same reasoning as in the case of ternary complexes,

${\displaystyle {\displaystyle {\displaystyle \begin{array}{ll}\mathfrak{C}(L_{\text{C}},L_{\text{W,1}},L_{\text{W,2}},L_{\text{W,3}},i,j,k)=\mathfrak{C}(L_{\text{C}})& \times \left[\frac{\mathfrak{C}(L_{\text{W,1}})}{4^{L_{\text{o,1}}}}\mathbb{1}\left(L_{\text{o,1}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+\mathbb{1}\left(L_{\text{o,1}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right]\\ & \times \left[\frac{\mathfrak{C}(L_{\text{W,2}})}{4^{L_{\text{o,2}}}}\mathbb{1}\left(L_{\text{o,2}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+\mathbb{1}\left(L_{\text{o,2}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right]\\ & \times \left[\frac{\mathfrak{C}(L_{\text{W,3}})}{4^{L_{\text{o,3}}}}\mathbb{1}\left(L_{\text{o,3}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+\mathbb{1}\left(L_{\text{o,3}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right].\end{array}}}}$

As 3-1 quaternary complexes are not symmetric under rotation, no symmetry correction of the combinatorial multiplicity is necessary.

Left-tilted 2-2 quaternary complexes

A 2-2 quaternary complex is comprised of two C strands and two W strands. We call a 2-2 complex left-tilted if strand W1 is connected to strand W2 via strand C1 (Figure 12A). The lengths of the oligomers are called $L_{\text{W,1}},L_{\text{W,2}},L_{\text{C,1}}$, and $L_{\text{C,2}}$. The positions of the strands relative to each other are governed by the alignment indices. All positions are measured relative to the position of the left end of strand C1. The alignment indices may take on the following values,

${\displaystyle {\displaystyle {\displaystyle i\in [-(L_{\text{W,1}}-1),L_{\text{C},1}-L_{\text{W,1}}-1],j\in [i+L_{\text{W,1}},L_{\text{C,1}}],\text{and}k\in [L_{\text{C},2},j+L_{\text{W,2}}-1].}}}$

The complex can have dangling ends on the right and on the left end of the complex; the presence of these dangling ends is indicated by the boolean variables $d_l$ and $d_r$. Moreover, two gaps are possible: There might be a gap between strands W1 and W2, or a gap between C1 and C2. The respective boolean variables read

${\displaystyle {\displaystyle \begin{array}{rlrl}& d_l=\{\begin{array}{ll}0& \text{if }i=0,\\ 1& \text{otherwise},\end{array}& & d_{m1}=\{\begin{array}{ll}0& \text{if }i+L_{\mathrm{W},1}=j,\\ 1& \text{otherwise},\end{array}\\ & d_{m2}=\{\begin{array}{ll}0& \text{if }{L}_{\mathrm{C},1}=k,\\ 1& \text{otherwise},\end{array}& & d_r=\{\begin{array}{ll}0& \text{if }j+L_{\mathrm{W},2}=k+L_{\mathrm{C},2},\\ 1& \text{otherwise}.\end{array}\end{array}}}$

We refer to the hybridization region of strand W1 and C1 as overlap region 1, to the hybridization region of strand W2 and C1 as overlap region 2 and to the hybridization region of strand W2 and C2 as overlap region 3. The length of these hybridization regions is computed via

${\displaystyle {\displaystyle {\displaystyle L_{\text{o,1}}=i+L_{\text{W,1}}-\text{max}(i,0),L_{\text{o,2}}=L_{\text{C,1}}-j,\text{and}{L}_{\text{o,3}}=\text{min}(j+L_{\text{W,2}},k+L_{\text{C,2}})-k.}}}$

Given the length of the hybridization region as well as the presence/absence of dangling ends, we can compute the hybridization energy,

${\displaystyle {\displaystyle {\displaystyle E=\gamma (L_{\text{o,1}}+L_{\text{o,2}}+L_{\text{o,3}}-3)+\gamma (1-d_{m1})+\gamma (1-d_{m2})+\frac{\gamma }{2}(d_l+d_r+2d_{m1}+2d_{m2}).}}}$

The combinatorial multiplicity of a left-tilted 2-2 quaternary complex is constructed using the same reasoning as in the case of a 3-1 complex,

${\displaystyle {\displaystyle {\displaystyle \begin{array}{ll}\mathfrak{C}(L_{\text{C,1}},L_{\text{C,2}},L_{\text{W,1}},L_{\text{W,2}},i,j,k)=\mathfrak{C}(L_{\text{C}})& \times \left[\frac{\mathfrak{C}(L_{\text{W,1}})}{4^{L_{\text{o,1}}}}\mathbb{1}\left(L_{\text{o,1}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+\mathbb{1}\left(L_{\text{o,1}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right]\\ & \times \left[\frac{\mathfrak{C}(L_{\text{W,2}})}{4^{L_{\text{o,2}}}}\mathbb{1}\left(L_{\text{o,2}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+\mathbb{1}\left(L_{\text{o,2}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right]\\ & \times \left[\frac{\mathfrak{C}(L_{\text{C,1}})}{4^{L_{\text{o,3}}}}\mathbb{1}\left(L_{\text{o,3}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+\mathbb{1}\left(L_{\text{o,3}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right].\end{array}}}}$

To prevent double-counting the same quaternary complex, we include either the complex or its rotated representation in the container of possible complexes, but not both. If the complex is symmetric under rotation, we multiply the combinatorial multiplicity by $1/2$. Given a left-tilted 2-2 quaternary complex $(L_{\text{C,1}},L_{\text{C,2}},L_{\text{W,1}},L_{\text{W,2}},i,j,k)$, we can compute the corresponding rotated complex $(L_{\text{C,1}}^{\text{rot}},L_{\text{C,2}}^{\text{rot}},L_{\text{W,1}}^{\text{rot}},L_{\text{W,2}}^{\text{rot}},i^{\text{rot}},j^{\text{rot}}{,}^{\text{rot}})$ by applying a linear map. The mapping of the oligomer lengths is illustrated in Figure 12B. We see that strand C2 after rotation corresponds to strand W1 before rotation, implying that $L_{\text{C,2}}^{\text{rot}}=L_{\text{W,1}}$. The same reasoning can be applied to all strands leading to the map,

${\displaystyle {\displaystyle {\displaystyle L_{\text{C,1}}^{\text{rot}}=L_{\text{W,2}},L_{\text{C,2}}^{\text{rot}}=L_{\text{W,1}},L_{\text{W,1}}^{\text{rot}}=L_{\text{C,2}},\text{and}{L}_{\text{W,2}}^{\text{rot}}=L_{\text{C,1}}.}}}$

In order to compute the map of the alignment indices under rotation, we need to express the relative positions of the strands with respect to the position of strand C1 after rotation, which corresponded to W2 before rotation. For example, $|i^{\text{rot}}|$ corresponds to the number of nucleotides by which strand C2 (before rotation) protrudes beyond strand W2 (before rotation). Expressed in terms of variables before rotation, this distance may be written as $k+L_{\text{C,2}}-j{L}_{\text{W,2}}$. Analogous relations can be derived for all alignment indices,

${\displaystyle {\displaystyle {\displaystyle i^{\text{rot}}=j-k+L_{\text{W,2}}-L_{\text{C,2}},j^{\text{rot}}=j+L_{\text{W,2}}-L_{\text{C,1}},\text{and}{k}^{\text{rot}}=j-i+L_{\text{W,1}}-L_{\text{W,2}}.}}}$

Right-tilted 2-2 quaternary complexes

A 2-2 quaternary complex is called right-tilted if strand W1 is connected to strand W2 via strand C2 (Figure 13). As in the case of the left-tilted 2-2 complex, the oligomer lengths are again called $L_{\text{W,1}},L_{\text{W,2}},L_{\text{C,1}}$ and $L_{\text{C,2}}$, but the values of the alignment indices that are possible for the right-tilted quaternary complex differ from the ones of the left-tilted complex,

${\displaystyle {\displaystyle {\displaystyle i\in [L_{\text{C,1}}-L_{\text{W,1}}-1,L_{\text{C,1}}-1],j\in [L_{\text{C,1}},i+L_{\text{W,1}}-1],\text{and}k\in [i+L_{\text{W,1}},k+L_{\text{C,2}}-1].}}}$

Note that the range of $i$ is chosen such that at least one nucleotide of strand W1 always extends to the right beyond the end of strand C1, allowing for a hybridization region between strand C2 and W1. The boolean variables denoting the presence/absence of dangling ends read

${\displaystyle {\displaystyle \begin{array}{rlrl}& d_l=\{\begin{array}{ll}0& \text{if }i=0,\\ 1& \text{otherwise,}\end{array}& & d_{m1}=\{\begin{array}{ll}0& \text{if }i+L_{\mathrm{W},1}=j,\\ 1& \text{otherwise,}\end{array}\\ & d_{m2}=\{\begin{array}{ll}0& \text{if }{L}_{\mathrm{C},1}=k,\\ 1& \text{otherwise,}\end{array}& & d_r=\{\begin{array}{ll}0& \text{if }j+L_{\mathrm{W},2}=k+L_{\mathrm{C},2},\\ 1& \text{otherwise.}\end{array}\end{array}}}$

The length of the overlap regions is given by

${\displaystyle {\displaystyle {\displaystyle L_{\text{o,1}}=L_{\text{C,1}}-\text{max}(i,0),L_{\text{o,2}}=i+L_{\text{W,1}}-k,\text{and}{L}_{\text{o,3}}=\text{min}(k+L_{\text{C,2}},j+L_{\text{W,2}})-j.}}}$

Like in the case of left-tilted 2-2 quaternary complexes, the total hybridization energy is computed via

${\displaystyle {\displaystyle {\displaystyle E=\gamma (L_{\text{o,1}}+L_{\text{o,2}}+L_{\text{o,3}}-3)+\gamma (1-d_{m1})+\gamma (1-d_{m2})+\frac{\gamma }{2}(d_l+d_r+2d_{m1}+2d_{m2}),}}}$

and the combinatorial multiplicity via

${\displaystyle {\displaystyle {\displaystyle \begin{array}{ll}\mathfrak{C}(L_{\text{C,1}},L_{\text{C,2}},L_{\text{W,1}},L_{\text{W,2}},i,j,k)=\mathfrak{C}(L_{\text{C}})& \times \left[\frac{\mathfrak{C}(L_{\text{W,1}})}{4^{L_{\text{o,1}}}}\mathbb{1}\left(L_{\text{o,1}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+\mathbb{1}\left(L_{\text{o,1}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right]\\ & \times \left[\frac{\mathfrak{C}(L_{\text{W,2}})}{4^{L_{\text{o,2}}}}\mathbb{1}\left(L_{\text{o,2}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+\mathbb{1}\left(L_{\text{o,2}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right]\\ & \times \left[\frac{\mathfrak{C}(L_{\text{C,1}})}{4^{L_{\text{o,3}}}}\mathbb{1}\left(L_{\text{o,3}}<L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)+\mathbb{1}\left(L_{\text{o,3}}\ge L_{\mathrm{U}}^{\mathrm{m}\mathrm{i}\mathrm{n}}\right)\right].\end{array}}}}$

We include either the quaternary complex or its rotated representation in the list of possible complexes to avoid double-counting. Moreover, the combinatorial multiplicity is divided by 2 if the complex is symmetric under rotation. It turns out that the rotation map for the right-tilted 2-2 quaternary complex is identical to the one of the left-tilted 2-2 complex,

${\displaystyle {\displaystyle \begin{array}{ll}& L_{\text{C,1}}^{\text{rot}}=L_{\text{W,2}},L_{\text{C,2}}^{\text{rot}}=L_{\text{W,1}},L_{\text{W,1}}^{\text{rot}}=L_{\text{C,2}},L_{\text{W,2}}^{\text{rot}}=L_{\text{C,1}},\\ & i^{\text{rot}}=j-k+L_{\text{W,2}}-L_{\text{C,2}},j^{\text{rot}}=j+L_{\text{W,2}}-L_{\text{C,1}},\text{and}{k}^{\text{rot}}=j-i+L_{\text{W,1}}-L_{\text{W,2}}.\end{array}}}$