Blog

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,

${\displaystyle {\displaystyle y=\frac{\mathrm{\#}\mathrm{n}\mathrm{u}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{o}\mathrm{t}\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{V}\mathrm{C}\mathrm{G}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{g}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{s}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}{\tau }_{\text{lig}}}{\mathrm{\#}\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{n}\mathrm{u}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{o}\mathrm{t}\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}{\tau }_{\text{lig}}}.}}$

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.

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,

${\displaystyle {\displaystyle {\displaystyle (\mathrm{\#}\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{n}\mathrm{u}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{o}\mathrm{t}\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}{\tau }_{\text{lig}})={\int }_0^{{\tau }_{\text{lig}}}\mathrm{d}t\sum _{\text{C}\in {\mathcal{E}}_t}N(\text{C})\text{min}(L_{\text{e,1}},L_{\text{e,2}})\mathbb{1}(\text{C allows templated ligation}).}}}$

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}}$,

(6)

${\displaystyle {\displaystyle {\displaystyle (\mathrm{\#}\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{p}\mathrm{o}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{n}\mathrm{u}\mathrm{c}\mathrm{l}\mathrm{e}\mathrm{o}\mathrm{t}\mathrm{i}\mathrm{d}\mathrm{e}\mathrm{s}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{l}{\tau }_{\text{lig}})\approx {\tau }_{\text{lig}}⟨\sum _{\text{C}\in \mathcal{E}}N(\text{C})\text{min}(L_{\text{e,1}},L_{\text{e,2}})\mathbb{1}(\text{C allows templated ligation})⟩,}}}$

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}$.

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

${\displaystyle {\displaystyle {\displaystyle y=\frac{⟨\sum _{\text{C}\in \mathcal{E}}N(\text{C})\text{min}(L_{\text{e,1}},L_{\text{e,2}})\mathbb{1}(\text{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})\mathbb{1}(L_{\text{e,1}}+L_{\text{e,2}}\ge L_{\mathrm{U}})⟩}{⟨\sum _{\text{C}\in \mathcal{E}}N(\text{C})\text{min}(L_{\text{e,1}},L_{\text{e,2}})\mathbb{1}(\text{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 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

${\displaystyle {\displaystyle {\displaystyle f=\frac{⟨\sum _{\text{C}\in \mathcal{E}}N(\text{C})\text{min}(L_{\text{e,1}},L_{\text{e,2}})\mathbb{1}(\text{C allows templated ligation})\mathbb{1}(L_{\text{e,1}}+L_{\text{e,2}}\ge L_{\mathrm{U}})\mathbb{1}(\text{product correct})⟩}{⟨\sum _{\text{C}\in \mathcal{E}}N(\text{C})\text{min}(L_{\text{e,1}},L_{\text{e,2}})\mathbb{1}(\text{C allows templated ligation})\mathbb{1}(L_{\text{e,1}}+L_{\text{e,2}}\ge L_{\mathrm{U}})⟩}.}}}$

Multiplying fidelity and yield results in the efficiency of replication,

${\displaystyle {\displaystyle {\displaystyle \eta =\frac{⟨\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})\mathbb{1}(L_{\text{e,1}}+L_{\text{e,2}}\ge L_{\mathrm{U}})\mathbb{1}(\text{product correct})⟩}{⟨\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 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,

${\displaystyle {\displaystyle s(\text{type})=\frac{⟨\sum _{\text{C}\in \mathcal{E}}N(\text{C})\text{min}(L_{\text{e,1}},L_{\text{e,2}})\mathbb{1}(\text{C allows templated ligation of given type)}⟩}{⟨\sum _{\text{C}\in \mathcal{E}}N(\text{C})\text{min}(L_{\text{e,1}},L_{\text{e,2}})\mathbb{1}(\text{C allows templated ligation)}⟩}.}}$

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,

${\displaystyle {\displaystyle {\displaystyle {\sigma }_{\mathcal{Z}}=\sqrt{\frac{{\sigma }_{⟨X⟩}^2}{⟨Y{⟩}^2}+\frac{⟨X{⟩}^2{\sigma }_{⟨Y⟩}^2}{⟨Y{⟩}^4}-\frac{2⟨X⟩{\sigma }_{⟨X⟩,⟨Y⟩}^2}{⟨Y{⟩}^3}}.}}}$

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.