Next Article in Journal
Photobiomodulation Therapy Applied after 6 Months for the Management of a Severe Inferior Alveolar Nerve Injury
Next Article in Special Issue
A Mutation Threshold for Cooperative Takeover
Previous Article in Journal
House Dust Mite Precision Allergy Molecular Diagnosis (PAMD@) in the Th2-prone Atopic Dermatitis Endotype
Previous Article in Special Issue
Computational Analysis of a Prebiotic Amino Acid Synthesis with Reference to Extant Codon–Amino Acid Relationships
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Hitting Times of Some Critical Events in RNA Origins of Life

by
Caleb Deen Bastian
1,*,† and
Hershel Rabitz
1,2,†
1
Program in Applied and Computational Mathematics, Princeton University, Princeton, NJ 08544, USA
2
Department of Chemistry, Princeton University, Princeton, NJ 08544, USA
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Life 2021, 11(12), 1419; https://doi.org/10.3390/life11121419
Submission received: 1 November 2021 / Revised: 9 December 2021 / Accepted: 11 December 2021 / Published: 17 December 2021

Abstract

:
Can a replicase be found in the vast sequence space by random drift? We partially answer this question through a proof-of-concept study of the times of occurrence (hitting times) of some critical events in the origins of life for low-dimensional RNA sequences using a mathematical model and stochastic simulation studies from Python software. We parameterize fitness and similarity landscapes for polymerases and study a replicating population of sequences (randomly) participating in template-directed polymerization. Under the ansatz of localization where sequence proximity correlates with spatial proximity of sequences, we find that, for a replicating population of sequences, the hitting and establishment of a high-fidelity replicator depends critically on the polymerase fitness and sequence (spatial) similarity landscapes and on sequence dimension. Probability of hitting is dominated by landscape curvature, whereas hitting time is dominated by sequence dimension. Surface chemistries, compartmentalization, and decay increase hitting times. Compartmentalization by vesicles reveals a trade-off between vesicle formation rate and replicative mass, suggesting that compartmentalization is necessary to ensure sufficient concentration of precursors. Metabolism is thought to be necessary to replication by supplying precursors of nucleobase synthesis. We suggest that the dynamics of the search for a high-fidelity replicase evolved mostly during the final period and, upon hitting, would have been followed by genomic adaptation of genes and to compartmentalization and metabolism, effecting degree-of-freedom gains of replication channel control over domain and state to ensure the fidelity and safe operations of the primordial genetic communication system of life.

1. Introduction

The origins of life, abiogenesis, is a matter of high importance, for it gives insight into the distribution of life in the universe. We focus on the RNA world hypothesis, where life began with self-replicating RNA molecules that can evolve under Darwinian evolution, following necessary conditions of compartmentalization and metabolism, for geometry and synthesis of nucleobases from metabolic precursors, respectively. Self-replicating sets of RNA were proposed first by Tibor Ganti [1,2] and have been studied by many others [3,4,5,6]. This is an information-centric perspective on abiogenesis, representing the putative beginning of genomic Darwinian evolution. Information centrism interprets a living organism as an operating genetic communication system in some connected domain that encodes and decodes genomic state relative to a replication channel.
While genomics, epigenomics, and transcriptomics of modern-day organisms are based on DNA, RNA, and epigenetic marks such as DNA methylation, RNA origins in their purest form concern the dual-function of RNA as an informational polymer and ribozyme. This article is similar in spirit to works in the 1970s through the 1990s, including Manfred Eigen’s works on replicating sets of RNA [7,8]. Clues to the RNA world, among others, are found in the nucleotide moieties in acetyl coenzyme A and vitamin B12, the structure of the ribosome as a ribozyme [9] and moreover the centrality of RNA to the translation system, and the existence of viroids [10]. A putative canonical RNA origins sequence involves RNA dependent RNA polymerase (RdRp) ribozyme, whose first gene is perhaps the Hammerhead (HH) ribozyme, enabling rolling-circle amplification of the sequence [11]. This setup is unnecessary, as the first role of RNA may not have been template-assisted polymerization but instead based on mechanisms of RNA recombination and networking [12]. If we assume the RdRp sequence is size 200 nucleotides, then there are 4 200 10 120 sequences. Starting from some population of interacting RNA molecules, we are interested in the times of first occurrence of critical events. Evolution seems to have concluded this search for a high-fidelity replicator in a fairly short period of time, i.e., within 400 million years of the Earth having a stable hydrosphere [5]. RdRp’s are known to be very ancient enzymes, are necessary to all viruses with RNA genomes, and have been proposed to have originated from junctions of proto-tRNAs relative to the context of a primitive translation system [13]. Recent work has shown that replicative RNA and DNA polymerases have a common ancestor of a RdRp [14]. Directed evolution, selecting on polymerization, is a potential way of identifying such a ribozyme. Directed compartmentalized self-replicating systems (RNA or DNA polymerases) mimicking prebiotic evolution have been demonstrated whereby polymerases are selected on their ability to replicate their own encoding gene [15]. Directed evolution has identified a RdRp that can replicate its evolutionary ancestor, an RNA ligase ribozyme; however, at increased activity, it has reduced fidelity and cannot maintain the integrity of its information [16].
A subtlety to the RNA origins argument is that template-directed polymerization by its nature requires two copies of the RdRp sequence, one for the polymerase and another for the template. This makes the RNA origins search extend until two copies are discovered. Cross reactions with other species also may influence the search time for the high-fidelity replicators. The clay mineral montmorillonite, which is common on Earth, can catalyze RNA oligomerization [17]; however, the utility of montmorillonite in these activities is not thought to be sufficient for origins, having been extensively studied [18,19]. Interestingly, it has been proposed that clay not only promotes origins, but constitutes it, which then later gave rise to RNA-based replication [20]; such mineral life as a genetic communication system has a high mutation rate and is degenerate.
Theories of abiogenesis study metabolism [21], cellular compartmentalization [22,23,24], hydrothermal vent chemical gradient energy [25], or hot springs [26,27,28], and so on, to define geochemical settings suitable for origins. These settings are compatible with RNA origins. Compartmentalization leading to selection on random sequences has been explored by studying environment forcing in hot springs and their effects on sequence identification [26,27,28], where the hypothesis is that fluctuations in environment forcing through cycling of wet, dry, and moist phases of lipid-encapsulated sequences subject sequences to combinatorial selection and identify structural and catalytic functions from the initial system state of random sequences. These functions include metabolic activity, pore formation, and structural stabilization. We assume the prebiotic molecular inventories of RNA and its precursors are provided by meteorites [29] and/or by Miller–Urey processes [30], such as from formamide [31] or many phase synthesis [32,33]. For energy and environmental factors, we consider a variant of Darwin’s “warm little pond”, where the putative environment for RNA origins of life is an icy pond with geothermal activity, a hot spring, or perhaps a hydrothermal vent: ice and cold temperature facilitate complexing of single strands into double strands and polymerization [34], and heat (energy) facilitates dissociation of double strands into single strands. More information on abiotic sources of organic compounds, mechanisms of synthesis and function of macromolecules, energy sources, and environmental factors can be found in the literature [35].
RNA origins have attracted many modeling efforts and analyses [36,37]. The concept of a self-replicating set of RNA molecules was initially studied by Manfred Eigen [7], wherein he studied the error-threshold of the critical fidelity to the main information. Two-dimensional spatial modeling has been applied in which reactions occur locally with finite diffusion, suggesting a spatially localized stochastic transition [38], simulated using Gillespie’s stochastic simulation algorithm (SSA) [39]. Another model is of autocatalytic sets of collectively reproducing molecules, which has been developed in reflexively autocatalytic food-generated (RAF) theory [40]. Various physics-based analyses have been conducted, such as in light of Bayesian probability, thermodynamics, and critical phenomena [41]. Systems of quasi-species based on the principle of natural self-organization called hypercycles employ non-equilibrium auto-catalytic reactions [8]. Nonlinear kinetic models for polymerization have been used to study the emergence of self-sustaining sets of RNA molecules from monomeric nucleotides [42]. Theoretical analysis has been conducted into RNA origins. Attention has been drawn to an evolving population of dynamical systems and how dynamics affect the error threshold of early replicators and possibly towards compartmentalization conveying hypercycles [43]. String-replicator dynamics have been studied and properties suggested to be necessary to RNA origins, including the ability to operate a functional genetic communication system and ecological and evolutionary stability [44,45,46]. A variety of pre-RNA worlds have been suggested, with RNA being preceded or augmented by alternative informational polymers, such as other nucleic acids [47], beta amyloid [48], polycyclic aromatic hydrocarbons [49], lipids [24], peptides [50], and so on. It seems that pre-RNA worlds existed independent of the RNA world in the sense that they are not ancestral to the RNA world, and that these worlds may have had non-trivial interactions with the RNA world.
A key concept of stochastic systems is that of a hitting time: the time of the first occurrence of some event. We develop a simple mathematical model at the sequence level to represent the synthesis and function of RNA molecules in order to gain insight into the hitting times of various critical events of RNA origins of life. The idea is to study the surface of hitting times in terms of the structure of the system. The model lacks many features of realism, such as sequence size variability, finite sources of “food” (activated nucleotides in our context), limited diffusion rates, poor system mixing, and so on, in order to concentrate on the process as a search problem. The notion of fitness landscapes has been studied extensively in evolutionary biology [51]. Landscape topology has been considered in an Opti-Evo theory, which assumes sufficient environmental resources and argues that fitness landscapes do not contain “traps” and globally optimal sequences form a connected level-set [52].
We describe the model in Section 2, where we define a replicating reaction network, whose random realizations are constructed using SSA. We describe hitting times as key random variables of interest and characterize polymerization as a transition kernel. In Section 3, we conduct and discuss simulation studies based on SSA, where we analyze the structure of the polymerase measures and the input–output and survival behaviors of the hitting times given the parameters of the system. In Section 4, we end with conclusions.

2. Materials and Methods

We define M = { Adenine , Uracil , Guanine , Cytosine } , abbreviated as M = { A , U , G , C } , corresponding to the RNA bases. We study the space of sequences having length n, that is, E = M n with space of possibilities (a σ -algebra) E = 2 E , so that the pair ( E , E ) is the measurable space of all RNA sequences of length n with | E | = 4 n . We let ν be a probability measure (distribution) on ( E , E ) , giving the probability space triple ( E , E , ν ) . Appendix A gives an overview of ( E , E , ν ) . A related space, though not utilized in this article, is the space of all RNA sequences up to length n, E * = i = 1 n H i with E * = 2 E * and size | E * | = 4 3 ( | E | 1 ) . We denote the collection of non-negative E -measurable functions by E 0 .
We build a simplified mathematical model for the time-evolution of a population of interacting RNA molecules in solution. Let X t be the population at time t R 0 with initial population X 0 . X t is a multiset, that is, it is a set containing elements possibly with repeats. We assume that the system is well mixed and has access to an infinite source of activated nucleotides.
The complement of x E is denoted x c E , attained using the base-pairing A with U and C with G. Let
h ( x , y ) { 0 , 1 , , n } for x , y E
be the Hamming distance between x , y E as the number of positions in x and y where the nucleotides differ. We have h ( x , x ) = 0 and h ( x , x c ) = n .

2.1. Core Model

We describe the reaction network of the system below. We simulate trajectories of the system using the stochastic simulation algorithm (SSA) [39], to simulate exact trajectories for the evolution of stochastic reaction networks here. SSA forms a Markovian process, where the arrival of reactions follows a Poisson (point) process, and assumes that the reaction volume is well mixed and homogeneous, with all parts of the system accessible for reactions. Reactions across various disjoint volume elements of the system are dependent. We do not consider the effects of finite diffusion, which effects a length scale above which disjoint volume elements are effectively independent [38].

2.1.1. System

We model the population of sequences which can form double-stranded helices, dissociate, and replicate with mutation with replicator fitness and sequence specificity. We interpret each system element as a set, either containing one element—a single-stranded sequence—indicated as { x } —or two elements, single and complementary stranded sequences, indicated as { x } { x c } = { x , x c } , where x , x c E (double bracket notation indicates a collection, or set, of sets, i.e., { { x } , { y } , { z } , } ). We define system elements as sets
E ¯ { { x } : x E } F ¯ { { x , x c } : x E } G ¯ E ¯ F ¯
with respective σ -algebras, E ¯ = 2 E ¯ , F ¯ = 2 F ¯ and G ¯ = 2 G ¯ . The sizes are | E ¯ | = 4 n and | F ¯ | = 4 n / 2 , and | G ¯ | = 3 × 4 n / 2 . The set E ¯ contains single-stranded sequences, whereas the set F ¯ contains double-stranded sequences; the set G ¯ is the union of E ¯ and F ¯ , containing both single and double stranded sequences. These sets are necessary to track the various sequences (single and double stranded). The reaction network of the system is given by
x + x c k d s x x c
x x c k s s x + x c
x + y k r e p ( x , y ) x + y + y * c .
and expressed in terms of sets
{ x } , { x c } k d s { x , x c } { x , x c } k s s { x } , { x c } { x } , { y } k r e p ( x , y ) { x } , { y } , { y * c } .
Reaction (2) is double-strand formation from complementary single-strands with reaction rate k d s . Reaction (3) is the dissociation of double-strands into single-strands, caused by a heat source, with reaction rate k s s . Reaction (4) is template-directed polymerization of a single-strand (the template) by another single-strand (the polymerase), producing a single-strand complementary to the template with some fidelity, with reaction rate k r e p (which functionally depends on the polymerase and template sequences).
The polymerization reaction rate is defined by
k r e p ( x , y ) = a f ( x ) s ( x , y ) ( 0 , a ] for x , y E
where a > 0 is a positive constant, f : E ( 0 , 1 ] is the replicative fitness of x (as a polymerase) and s : E × E ( 0 , 1 ] is the similarity between x and y, a symmetric function. Finally, x replicates y, outputting a version y * c with mutations, where each nucleotide position has fidelity probability p : E ( 0 , 1 ] . Note that the similarity function can be either trivial/constant, i.e., s ( x , y ) = 1 , or non-trivial. For example, if we assume sequence similarity to correlate to spatial proximity of sequences, as assumed below, then s ( x , y ) is non-trivial.

2.1.2. High-Fidelity Set

To define a high-fidelity set, pick an arbitrary subset of sequences R E as high-fidelity replicators of size r = | R | . We define R two ways.
We define R using a product of non-empty random nucleotide subsets { A i M : i = 1 , , n } for each nucleotide position
R = A 1 × × A n
so that r = i = 1 n | A i | = 1 r 1 2 r 2 3 r 3 4 r 4 where r 2 = | { A i : | A i | = 2 } | , etc. and r 1 + r 2 + r 3 + r 4 = n . For simplicity, we assume | A i | { 1 , 4 } with fraction 4 being q ( 0 , 1 ) . Thus, R is a subset of E defined as a product space.
For another construction of R, we define a finite union of m random sequences R = { x 1 , , x m } .

2.1.3. Distance

We define the Hamming distance H between sequence x and high-fidelity sequence set R as
H ( x , R ) = min { h ( x , y ) : y R } { 0 , 1 , , n } for x E .

2.1.4. Fitness

We define “tent-pole” fitness f k of sequence x E and high-fidelity sequence set R for curvature parameter k R 0 as
f k ( x , R ) = exp [ k H ( x , R ) ] ( 0 , 1 ] for x E .
The maximums are the sequences of the high-fidelity sequence set R , which are the “points” or “poles” of the surface, with exponential decay into the remainder of the space in string distance. The strength of the decay is governed by parameter k, called the curvature parameter, which can be specified through the value of fitness at H ( x , R ) = n (sequence dimension),
k = log ( f k ( x , R ) ) n
Appendix B describes other fitness functions.

2.1.5. Similarity

We define two cases for the similarity function appearing in the reaction rate of template-directed polymerization. The first case is the trivial (constant) case where the
s b ( x , y ) = b ( 0 , 1 ] for ( x , y ) E × E .
This assumes that there is no mechanism by which sequence specificity is selected for, such as in the case that polymerases should evolve to generically well replicate sequences, including their own. This means that they will spend much of their time replicating other sequences.
For the second case of a non-trivial similarity function, we note that RNA origins of life are thought to be a spatially localized stochastic transition, where high-fidelity replicators are found concentrated in foci, following from the increased replicative mass of the replicators. Hence, we implicitly encode spatial information through a non-trivial similarity function, based on a distance function, which increases the replicative system mass for similar (and here nearby) high-fidelity replicators, that is, if they’re similar, then they’re likely proximal. In the following definition, we use the notation ∧ for the minimum of two numbers, x y = min { x , y } . Distance S is defined between two sequences x , y E as
S ( x , y ) = h ( x , y ) h ( x , y c ) h ( x c , y ) h ( x c , y c ) { 0 , 1 , , n } for x , y E .
Similarity s k of sequences x , y E for curvature parameter k R 0 is defined in terms of exponential decay as
s k ( x , y ) = exp [ k S ( x , y ) ) ] ( 0 , 1 ] for x , y E .
Presently, replicators can replicate other sequences well but not their own [34]. There may exist RdRps that are excellent polymerases and, in conjunction with RNA hammerhead ribozyme, engage in rolling circle amplification of the polymerase-hammerhead sequence (genome) so that the amplification process is self-cleaving. This results in a large increase in replicative mass due to the super-exponential growth in the population of the high-fidelity replicators within a small volume. In the context of SSA, the similarity function here is an ansatz for spatial locality.

2.1.6. Fidelity

Finally, polymerization fidelity probability for curvature k R 0 is defined as
p k ( x , R ) = f k ( x , R ) ( 0 , 1 ] for x E .
Note that f k ( x , R ) = p k ( x , R ) = 1 for high-fidelity sequences x R .
Note that fitness, similarity, and fidelity are defined for single-stranded sequences ( E , E ) .

2.1.7. Counting Representation

The process X = ( X t ) t R 0 is the time-evolution of the system. Recall that X t is a multiset. X t contains the individual single stranded molecules in the set E ¯ , i.e., sequences { x } E ¯ and double stranded molecules in the set { x , x c } F ¯ , with overall set G ¯ = E ¯ F ¯ having size m = | G ¯ | = 3 | E ¯ | / 2 . Note that, in the set representation, there is symmetry { x , x c } = { x c , x } , so that the size of the double-stranded set is equal to | F ¯ | = | E ¯ | / 2 . The system evolution X t induces a random counting measure N t on the overall space of single and double stranded sequences ( G ¯ , G ¯ ) as
N t ( A ) = x X t I A ( x ) for A G ¯
The total count, that is, the total number of molecules, is K t | X t | = N t ( G ¯ ) . We assume that the counter N t is maintained for all times t R 0 .

2.1.8. Reaction Rates

The total reaction rate is given by the sum of the individual reaction rates
k ˘ ( t ) = k ˘ d s ( t ) + k ˘ s s ( t ) + k ˘ r e p ( t ) for t R 0
where the reaction rate of double-strand formation from complementary single-strands is given by
k ˘ d s ( t ) = 1 2 { x } E ¯ k d s N t ( { x } ) N t ( { x c } ) ,
the reaction rate of dissociation of double-strands into single-strands is given by
k ˘ s s ( t ) = { x , x c } F ¯ k s s N t ( { x , x c } ) ,
and the reaction rate for template-directed polymerization of a single-strand by another single-strand (the polymerase) is given by
k ˘ r e p ( t ) = ( { x } , { y } ) E ¯ 2 k r e p ( x , y ) N t ( { x } ) ( N t ( { y } ) I ( x = y ) )
The first reaction rate is a sum over the single-strands of E ¯ with size 4 n . The second is a sum over double-strands F ¯ with size 4 n / 2 . The third reaction is a sum over the product space of single stranded sequences, E ¯ 2 , with 4 2 n = 16 n number of elements. Therefore, k ˘ ( t ) requires 4 n ( 3 2 + 4 n ) elements to be evaluated for every reaction. Clearly, direct representation on the full space is very expensive and impractical for even modest n. One obvious way to improve efficiency is not summing over the zero elements. We define sets
X 1 = { x set ( X t ) : | x | = 1 }
and
X 2 = { x set ( X t ) : | x | = 2 }
as the unique single and double stranded sequences of the system. Then, direct calculations of the reaction rates are
k ˘ d s ( t ) = 1 2 { x } X 1 k d s N t ( { x } ) N t ( { x c } )
and
k ˘ s s ( t ) = { x , x c } X 2 k s s N t ( { x , x c } )
and
k ˘ r e p ( t ) = ( { x } , { y } ) X 1 2 k r e p ( x , y ) N t ( { x } ) ( N t ( { y } ) I ( x = y ) ) .
For this approach, the replication rate has quadratic dependence on | X 1 | . Using the reaction rates, the system may be exactly simulated using SSA. The reaction at time t with rate k ˘ ( t ) occurs over time interval t Exponential ( 1 / k ˘ ( t ) ) . As the reaction rate k ˘ ( t ) increases with increasing number of molecules K t = N t ( G ¯ ) , the reaction rate increases and reaction duration t decreases over time. The natural consequence of increasing process intensity is that the system speeds up.
The quadratic dependence may still be too expensive for large simulations. Appendix E describes Monte Carlo approximation of the reaction rates.

2.2. Hitting Times

We define some hitting times. The initial population consists of I single-stranded sequences X 0 , i.e., | X 0 | = I . We define the hitting time τ for the time of the first replication event
τ rep = inf t R 0 : K t > I .
We define the hitting time τ for the appearance of sequences in the high-fidelity sequence set R ,
τ R = inf { t R 0 : I { 1 , 2 , } ( | R X t | ) = 1 } .
Put X R = R { { x , x c } : x R } and define the volume fraction of high-fidelity sequences R at time t as
V ( t ) = N t ( X R ) K t .
We define the hitting time τ where high-fidelity sequences of R emerge and reach a minimum volume fraction,
τ min = inf { t R 0 : t τ R , V ( t ) V ( s ) for τ R s t } .
This hitting time reflects that period wherein a high-fidelity replicator has been identified yet there exists no complementary high-fidelity sequence for amplification, hence the system diversity continues to increase, decreasing the concentrations of all extant sequences as more sequences are discovered. The minimum hitting time captures the duration of time the high-fidelity replicator exists by itself. We define the hitting time τ for the time high-fidelity sequences in R constitutes some volume fraction v ( 0 , 1 ] of the population,
τ v = inf t R 0 : V ( t ) v .
In practice for simulations, τ v is censored based on some total number of reactions, that is, if the volume fraction is not achieved by n reactions, τ v = because there is no arrival time.
For SSA, we specify a maximum number of reactions N to simulate. We have parameters θ Θ for τ , such as landscape curvature k, sequence dimension n, etc. Therefore, τ ( θ ) is right-censored with value at simulation time a, as some simulations will stop at time a with no arrival time. These are censoring events. For fixed θ , the τ ( θ ) is a random variable, due to the stochastic nature of SSA. Hence, for each parameter vector θ , we attain a set of M realizations of hitting time τ as
T ( θ ) = { τ i ( θ ) : i = 1 , , M } .
For convenience, we assume that the realizations T ( θ ) are ordered by non-censored followed by censored.

2.2.1. Functional Structure

For each parameter vector θ , we record two values: the number of hitting events in the hitting time set T
g ( θ ) = | { x T ( θ ) : x < } | { 0 , 1 , , M }
and the average of the hitting time τ
f ( θ ) = 1 g ( θ ) i = 1 g ( θ ) τ i ( θ ) if g ( θ ) > 0 0 if g ( θ ) = 0 R 0
To describe the functional structure of the average hitting time f ( θ ) , we require a classifier which determines whether or not there are zero hittings g ( θ ) = 0 and a regressor for the value of f ( θ ) for hittings g ( θ ) > 0 . We assume that the parameters θ = ( θ 1 , , θ n ) are randomly sampled according to distribution ν = i ν i and the hitting times recorded. High dimensional model representation (HDMR) may be attained for the classifier (as a probabilistic discriminative model) and the regressor of f ( θ ) . For the regressor, we have HDMR expansion
f ( θ 1 , , θ n ) = f 0 + i f i ( θ i ) + i < j f i j ( θ i , θ j ) + + f 1 n ( θ 1 , , θ 1 )
The HDMR component functions { f u } convey a global sensitivity analysis, where, defining variance term
σ u 2 = V a r f u = Θ f u 2 ( θ u ) ν ( d θ ) ,
we have a decomposition of variance
σ f 2 = V a r f = u { 1 , , n } : | u | > 0 σ u 2
The normalized terms S u = σ u 2 / σ f 2 are called sensitivity indices. Appendix F gives a brief description of global sensitivity analysis via HDMR.

2.2.2. Statistical Structure

A second analysis can be conducted on the hitting times T ( θ ) for the parameter vector θ using reliability theory. Put random hitting set T ( Θ ) { T ( θ ) : θ Θ } where Θ = { θ i } is an independency of parameter values. For each parameter vector θ Θ , we partition the hitting times T ( θ ) into C censored values with censor times C ( θ ) = { a i ( θ ) } and M C non-censored (hitting) values N ( θ ) = { x T ( θ ) : x < } . The likelihood is given by
L ( T ( Θ ) | ϑ ) = θ Θ x N ( θ ) f ( x | ϑ ) x C ( θ ) R ( x | ϑ ) ,
where f is the hitting time probability density function (‘failure density’) and R is censoring time distribution (‘reliability distribution’), and ϑ are the parameters of the density and distribution functions. Note that f and R each specify each other, so ϑ are the common parameters. Reliability definitions are given in Appendix G.
We show reliability quantities in Table 1 for the two-parameter ( α , β ) ( 0 , ) 2 Weibull ( α , β ) distribution and Cox proportional hazard’s model where γ is a vector of coefficients for the θ . We use the Python software lifelines for estimation of ϑ for the Weibull–Cox model from data [53]. The mean failure time E τ ( θ ) is given by
E τ ( θ ) = 0 t f ( t , θ | α , β , γ ) = β e γ · θ / α Γ ( 1 + 1 α )
We have second moment
0 t 2 f ( t , θ | α , β , γ ) = β 2 e 2 γ · θ / α Γ ( 1 + 2 α )
giving variance
V a r τ ( θ ) = β 2 e 2 γ · θ / α ( Γ ( 1 + 2 α ) Γ 2 ( 1 + 1 α ) )
Thus, if γ < 0 , then E τ ( θ ) and V a r τ ( θ ) exponentially increase in θ .

2.3. Surface Chemistries

The system given by (2) of polymerization with mutation requires two separate hitting events, one sequence in the high-fidelity set x R and either another sequence in the high-fidelity set x R or its complement x c E , in order for high-fidelity replicators to maximally engage in templated-directed polymerization and achieve some fraction of the population. This setup of RNA polymerase action, requiring two such events for the polymerase and template, makes the hitting times long. Basically, the same information must be discovered twice before it can be used, which is unsatisfactory. We idealize polymerase activity conveyed by a non-RNA species, here clay, with the parameter k c l a y p , as clay itself is not thought to be capable of polymerization but is capable of oligomerization of RNA. The reactions are given by
k c l a y o x
x k c l a y p x + x * c
where non-RNA polymerization has mutation with fidelity probability p ( 0 , 1 ] and x * c is the complement of x with mutation. The reaction rates are given by
k ˘ c l a y o ( t ) = k c l a y o
and
k ˘ c l a y p ( t ) = k c l a y p N t ( X 1 )
Therefore, upon the first hitting of the high-fidelity replicators R with sequence x R through (4) or (23), x gives two high-similarity single-stranded sequences x and x * c through (24), which then may participate in template-directed RNA polymerization (4).

2.4. Reactions as Measure-Kernel-Functions

All the reactions x y which involve substrate x may be represented using transition kernels, which form linear operators. At each iteration of SSA, a reaction type is chosen, followed by a transition to a particular domain ( X , X ) with distribution ν t , followed by mapping into a codomain ( Y , Y ) using transition probability kernel Q with distribution μ t = ν t Q . The notions of ν t and Q involve measure-kernel-functions. The probability of transition of y into B Y given x is given by Q
Q ( x , B ) = B Q ( x , d y ) = P ( y B : x )
Appendix C recalls some facts about Q.
We define kernels Q for RNA and non-RNA polymerization to provide insight into the reactions. Consider X t for some t R 0 . Recall that N t is the random counting measure of X t on single and double-stranded sequences ( E F , 2 E F ) . For RNA and non-RNA polymerization, we take ν t as a (random) probability measure for x in domain ( X , X ) and describe a transition probability kernel Q from x into y in codomain ( Y , Y ) .

2.4.1. RNA Polymerization

For RNA polymerization, we have that
{ x } , { y } { x } , { y } , { y * c }
which results in the creation of the single-stranded sequence { y * c } . The first dimension is the polymerase and the second is the template. We give a definition of the probability measure on the product space of sequences. Recall that ( E E ) 0 denotes the collection of non-negative E E -measurable functions.
Definition 1
(Measure on domain ν t ). Let ν t be a random probability measure on ( E × E , E E ) formed by random counting measure N t (12)
ν t { x , y } = k r e p ( x , y ) N t ( { x } ) ( N t ( { y } ) I ( x = y ) ) k ˘ r e p ( t ) for ( x , y ) E × E
with
ν t ( f ) = ( x , y ) E × E ν t { x , y } f ( x , y ) for f ( E E ) 0
We write ν t ( A ) = ν t I A for A E E .
Because the first two coordinates are preserved under the mapping, we focus on the new dimension as a transition from ( E × E , E E ) into ( E , E ) using transition probability kernel Q. In this case, Q is defined by a 16 n × 4 n matrix whose rows vectors (dimension 4 n ) are probability vectors. The structure of Q follows from the polymerase replication with mutation, whereby each nucleotide position has fidelity probability p : E ( 0 , 1 ] , which depends on the first dimension of E × E . We put p x = p ( x ) for sequence x E . Now, we state a simple fact on the binomial structure of the number of mutations made by a polymerase.
Theorem 1
(Mutation distribution). The number of mutations by polymerase x E on template y E is distributed
h ( y c , y * c ) Binomial ( n , 1 p x ) for p x ( 0 , 1 )
with mean n ( 1 p x ) and variance n p x ( 1 p x ) and
h ( y c , y * c ) Dirac ( 0 ) for p x = 1 .
Now, we partition E into level sets ( H i ( y ) ) by Hamming distance to the template complement y c ,
H i ( y ) = { x E : h ( y c , x ) = i } for i { 0 , , n } .
We define the transition kernel Q for RNA polymerization, where Q completely encodes RNA polymerization using Theorem 1.
Corollary 1
(Transition probability kernel Q). We have that the transition probability kernel Q for RNA polymerization is defined by
Q ( ( x , y ) , H i ( y ) ) = n i ( 1 p x ) i p x n - i for ( x , y ) E × E , i { 0 , , n } , p x ( 0 , 1 )
and
Q ( ( x , y ) , { z } ) = 1 | H i ( y ) | n i ( 1 p x ) i p x n i for ( x , y ) E × E , i = { 0 , , n } z H i ( y ) , p x ( 0 , 1 )
and
Q ( ( x , y ) , { y c } ) = 1 for ( x , y ) E × E , p x = 1 .
RNA polymerization is defined by Q using the binomial structure of polymerase mutation. A more sophisticated model could be defined as a sum of Bernoulli random variables with varying success probabilities in the Poisson binomial distribution. This could be used to take into account polymerase mutation that varies with nucleotide position. Another idea is taking into account schemata such as repeats which destabilize the polymerase [54].
Proposition 1
(Measure on codomain μ t ). μ t = ν t Q is a probability measure on ( E , E ) defined by
μ t ( f ) = E × E ν t ( d x , d y ) E Q ( ( x , y ) , d z ) f ( z ) for f E 0
It is multiplication of ν t as a 16 n dimension row vector with 16 n × 4 n dimension matrix Q, giving a 4 n dimension row vector ν t Q . We write μ t ( A ) = μ t I A for A E .
Define the partition ( H i ) of E as
H i { x E : min { H ( x , R ) , H ( x c , R ) } = i } for i { 0 , , n } .
Then, μ t ( H i ) for i { 0 , , n } is the distribution on sequences by distance to R, i.e.,
μ t ( H i ) = x E μ t { x } I H i ( x ) for i { 0 , , n }
contains the instantaneous information of RNA polymerization.
A more general model for replication is where polymerase activity is tied to geometry, i.e., compartmentalization/spatial confinement, with state space ( C , C ) and to metabolic state ( M , M ) . In this telling, the polymerase reaction rate could be tied to the degree of spatial confinement and the source of the activated nucleotides from metabolic precursors. Then, the polymerase state-space is ( C × M × E × E , C M E E ) with law ν t and the polymerase transition kernel Q c m e is defined as the mapping from ( C × M × E × E , C M E E ) into ( E , E ) . Thus, the law on the input–output space-state ( C × M × E × E × E , C M E E E ) is given by μ t = ν t × Q c e m , or in differential notation,
μ t ( d c , d m , d x , d y , d z ) = ν t ( d c , d m , d x , d y ) Q c e m ( ( c , m , x , y ) , d z )

2.4.2. Non-RNA Polymerization

If there exists some kind of non-RNA polymerase activity, we have that the mapping
{ x } { x } , { x * c }
which we regard as a mapping from ( E , E ) into ( E , E ) . Let ν t be a probability measure on ( E , E ) defined by
ν t { x } = N t ( { x } ) N t ( E ) for x E
Similar to RNA polymerization, for fidelity probability p ( 0 , 1 ] , we have Q as the 4 n × 4 n matrix defined by
Q ( x , H i ( x ) ) = n i ( 1 p ) i p n i for x E , i { 0 , , n } , p ( 0 , 1 )
and
Q ( x , { z } ) = 1 | H i ( x ) | n i ( 1 p ) i p n i for x E , i = { 0 , , n } , z H i ( x ) , p ( 0 , 1 )
and
Q ( x , { x c } ) = 1 for x E , p = 1 .
μ t = ν t Q is a probability measure on ( E , E ) defined by
μ t ( f ) = E ν t ( d x ) E Q ( x , d y ) f ( y ) for f E 0
Note that, for SSA, Q is fixed over the simulation, whereas the probability measure ν t depends on time. That is, the reactions are chosen according to the reaction rates, and the reactions each use respective Q. The ν t is formed using a random counting measure, so ν t is random. This approach generalizes in the obvious way to all the reactions.

2.5. Decay

The RNA sequences have finite lifetimes in reality. This comes from a variety of sources, including radiation, pH, intrinsic molecular stability, etc. We assume double-stranded RNA is stable, whereas single-stranded RNA is not. Therefore, we create a reaction for decay of single-stranded RNA into constitutive nucleotides
x k
with reaction rate
k ˘ ( t ) = k N t ( X 1 )

2.6. Compartmentalization

It is thought that compartmentalization plays a role in RNA origins of life, giving foci of reproducing sequences [23,55]. This is somewhat anticipated by the similarity function s : E × E ( 0 , 1 ] , where sequences are more likely to copy similar sequences than less similar ones, due to an underlying spatial localization. Explicit spatial effects may be modeled by assuming each x E is marked with a position on a bounded subset of the real line ( [ T , T ] , B [ T , T ] ) ( R , B R ) . We think of this as a one-dimensional projection of the three-dimensional system. Additional species can be introduced, such as lipids, with reactions forming a vesicle M (vesiculation), which encloses some A = [ r , s ] [ T , T ] . We assume the lipids interact with the single stranded sequences in A to form vesicles as
x k m i c M ( A )
with reaction rate
k m i c ( t ) = k m i c N t ( X 1 ) .
Hence, vesiculation is coupled to the population of sequences by design so that it evolves on roughly the same time-scale as sequence activities. Note that vesicles can enclose one another, i.e., M ( A ) and M ( B ) where A B or B A , but cannot cross, i.e., for all vesicles at locations A , , B we have that A B { A , B , } . For example, suppose one vesicle A = [ 0 , 1 ] encloses another two, B = [ 1 3 , 1 2 ] and C = [ 2 3 , 3 4 ] . Then, A \ ( B C ) = [ 0 , 1 3 ) ( 1 2 , 2 3 ) ( 3 4 , 1 ] . Although A \ ( B C ) is disconnected in one dimension, the intervals are physically connected in three dimensions, where vesicles are spheres. We identify each vesicle to a union of disjoint intervals, disjoint across the vesicles.
We posit that compartmentalization precedes the hitting of a high-fidelity replicator through ensuring necessary concentration of RNA and a stable environment. Generally, we identify compartmentalization state to ( C , C ) with probability measure ν . Let Q c be a transition probability kernel from ( C , C ) into ( E , E ) , encoding the transition from compartmentalization coordinates to RNA sequences. The product space ( C × E , C E ) has law μ = ν × Q c . Upon hitting a high-fidelity replicator and achieving Darwinian evolution to acquire information, e.g., genes, the sequences are assumed to become adapted to compartmentalization coordinates ( C , C ) through the transition probability kernel Q c from ( C × E , C E ) into ( C , C ) , so that μ = ν × Q c × Q c is the law on the full space ( C × E × C , C E C ) . In this telling, compartmentalization precedes RNA activity, and, upon hitting high-fidelity replicators that can maintain their information, is followed by genomic adaptation.

2.7. Metabolism

We identify metabolism reaction-state to the measurable space ( M , M ) with probability measure ν . Let Q m be a transition kernel from ( M , M ) into ( E , E ) , positing that metabolism precedes replication. For example, certain metabolic state may be precursors to the synthesis of RNA. Consider product space ( M × E , M E ) with measure μ = ν × Q m . Now we suppose that, upon achieving Darwinian evolution in replicators, the replicators will eventually become adapted to ( M , M ) . Hence, we interpret ( M , M ) as a mark-space of ( M × E , M E ) , representing genomic adaptation. Let Q m be a transition kernel from ( M × E , M E ) into ( M , M ) . Then, μ = ν × Q m × Q m is a probability measure on ( M × E × M , M E M ) , where
μ ( f ) = M ν ( d x ) E Q m ( x , d y ) M Q m ( ( x , y ) , d z ) f ( x , y , z ) for f ( M E M ) 0
or
μ ( d x , d y , d z ) = ν ( d x ) Q m ( x , d y ) Q m ( ( x , y ) , d z ) .
Therefore, metabolism-first followed by replication and genomic adaptation is encoded by the structures of Q m and Q m . We do not specify these transition kernels in this article but mention that they are richly textured.

2.8. Reaction Overview

The reactions of the system having decay and clay and their reaction orders are shown in Table 2. There is one zero-order reaction, three first-order reactions, and two second-order reactions. Additionally, we show reactions and orders for compartmentalization and metabolism.

3. Results

Consider some initial population of I random sequences X 0 . The population over time is given by X t with associated random counting measure N t on ( G ¯ , 2 G ¯ ) . Recall parameters
θ = ( n , q , k , l , m , p , k , k s s , k d s , k r e p , k c l a y o , k c l a y p )
for sequence dimension n, high-fidelity sequence set size q, fitness degree k, similarity degree l, fidelity degree m, clay fidelity probability p, RNA decay rate k , double-strand dissociation rate k s s , double-strand formation rate k d s , RNA replication rate k r e p , and clay oligomerization and polymerization rates k c l a y o and k c l a y p . These parameters are summarized in Table 3.
The following is the description of how the parameter values were specified and to what they biologically correspond. The sequence dimension n is chosen from { 3 , 4 , 5 } . The fitness and similarity functions are chosen by setting the value of the range of the curvature parameters k and l from one (inside the high-fidelity manifold) to some small values, such as over an exponential grid. For example, when i = 0.1 , the fitness of sequences that are maximally dissimilar have 10% of the fitness of the high-fidelity sequences. We range the grid from 0.1 to 0.001 for fitness and similarity. The RNA fidelity parameter m = 0.25 is chosen such that the high-fidelity sequences have value one and the lowest fidelity sequences have value 0.25, equal to random chance. The clay fidelity parameter is set to an optimistically high value of 0.9 for clay studies. The double-strand dissociation and formation rates k s s and k d s are set to unity as a baseline. In comparison, the RNA replication rate is set to a large value, 10, whereby replication is the dominant reaction. The decay parameter k is set to some uniform random value in (0, 1). The clay RNA oligomerization rate is set to unit, and ‘clay’ RNA polymerization rate is set to a uniform random value in (0, 20).
With the parameters governing the reaction rates, different values of these parameters confer different regimes for the system.

3.1. Stability: ODEs

We characterize the zeros of the vector field f from ODE system (A1) and use the eigenvalues of the Jacobian (A2) to determine their stability.
Theorem 2.
The ODE system (A1) for R = { x } , x E , has a single unstable fixed-point at [ x ] = 1 and [ y ] = 0 for y G \ x .
Proof. 
Solving f = 0 gives a single solution [ x ] = 1 and [ y ] = 0 for y G \ x . For this solution, the eigenvalues of the Jacobian contain no zero values and positive values. Therefore, the solution is unstable. □
It follows from Theorem 2 that, for all other initial conditions, the system has no equilibria.
Corollary 2
(Unbounded). For all initial conditions X 0 such that I = | X 0 | > 1 , the system is unbounded.
This confirms the obvious: the system, a replicating network with no death, is almost always an increasing system.

3.2. Simulation Reaction State

We are interested in the behavior of temporal probability measures, ν t (25) on the sequence product space and μ t = ν t Q (27) on the sequence space, for RNA polymerization. These reveal the instantaneous information of the system. The structure of μ t reveals the state of polymerization and is a leading indicator of the population concentrations over time.

3.2.1. Core Model with “Tent” Functions, Probable Hitting P ( τ v ( θ ) < ) 1

Take sequence dimension n = 3 and fitness and similarity curvature parameters k = l = log ( 0.01 ) / n and fidelity parameter m = log ( 0.25 ) / n . Set rates for double-strand dissociation and formation k s s = k d s = 1 and polymerization rate k r e p = 10 and use the “tent” function for fitness, similarity, and fidelity. Take random initial population X 0 with initial population size I = | X 0 | = 10 and random singleton R = { { x } } ( q = 0 ). We simulate 5000 reactions, simulation censored at hitting time τ v for volume fraction v = 0.25 . Take partition of the sequence space by Hamming distance to the high-fidelity manifold ( H i ) (28) of sequence space E. In Figure 1, we plot measures of a typical realization of the system X t on the partition ( H i ) of sequence concentration (Figure 1a), growth (Figure 1b), and polymerase sequence output μ t (Figure 1c). Some quantities are plotted on log-log scale, whereas others are plotted on a linear-log scale. These results show that the concentrations are relatively stable for most time, until the high-fidelity manifold is hit. Then, the concentration of high-fidelity replicators rapidly increases to exceed 25%. Similarly, Figure 1b shows the growth curves on a log-log scale, where the high-fidelity manifold rapidly increases near the end of the simulation. Figure 1c shows the structure of the RNA sequence polymerization output temporal probability measure μ t . Low probability is assigned to polymerization of high-fidelity replicators for most of the reaction time, followed by a large increase near the end of the simulation, where high-fidelity replicators dominate with 56% probability. Therefore, the RNA sequence polymerization output temporal probability measure μ t is a leading indicator of the concentration curve, i.e., at simulation end-time, concentration of high-fidelity replicators is 25% and polymerization output is 56%.

3.2.2. Core Model with “Tent” Functions, Improbable Hitting P ( τ v ( θ ) < ) 0

We use the same configuration as Section 3.2.1 except for setting fitness and similarity curvature parameters to k = l = log ( 0.1 ) / n . In Figure 2, we plot measures of a typical realization of the system X t on sequence partition by Hamming distance to the high-fidelity manifold ( H i ) of concentration (Figure 2a), growth (Figure 2b), and μ t (Figure 2c). The behavior has completely changed: the high-fidelity group ends the simulation with around 6% concentration, only steadily increasing, and never hits. The polymerase output μ t shows 6%. This indicates that the concentration of high-fidelity replicators is unlikely to increase further, as the population is generally in equilibrium with the polymerase output.

3.2.3. Core Model with Linear Functions, Improbable Hitting P ( τ v ( θ ) < ) 0

The same configurations for Section 3.2.1 are used, except the fitness, similarity, and fidelity functions are linear. Similar to the “tent” functions, we specify the terminus landscape curvature for fitness and similarity k = l . Then, the fitness function for RNA polymerization is given by
f k ( x , R ) = 1 + k 1 n H ( x , R ) for x E
and
s k ( x , y ) = 1 + k 1 n S ( x , y ) for x , y E
We put fitness and similarity landscape curvature k = l = 0.01 for fitness and similarity and v = 0.25 for hitting volume fraction of high-fidelity replicators. We simulate X t for 5000 reactions. We find that the probability of hitting is near zero P ( τ 0.25 ( θ ) < ) 0 . In Figure A1, we plot measures of a typical realization of X t on ( H i ) of concentration (Figure A1a), growth (Figure A1b), and μ t (Figure A1c). The simulation ends with high-fidelity concentration of ∼5% and polymerase output of ∼4%. Therefore, the concentration of high-fidelity replicators will continue to decrease. Linear surfaces are not sufficient to achieve hitting times τ v ( θ ) < for high-fidelity replicator volume fraction v = 0.25 , in contrast to the nonlinear “tent” functions.

3.2.4. Expanded Model with “Tent” Functions, Probable Hitting P ( τ v ( θ ) < ) 1

We consider a similar model to the previous subsections and expand it with clay oligomerization rate (of RNA) k c l a y o , clay polymerization rate (of RNA) k c l a y p , and clay polymerization fidelity p. Therefore, the full set of variables is given by θ = ( n , k s s , k d s , k , k c l a y o , k c l a y p , p ) . The value of fitness/similarity landscape curvature k , l and clay RNA polymerization raet k c l a y p are set such that the replicative mass of each is initialized to 10. This means that RNA and clay polymerization have the same reaction mass at the beginning of the simulation. We set sequence dimension n = 3 , fitness/similarity landscape curvature k = l = log ( 0.01 ) / n , clay RNA polymerization fidelity p = 0.9 , and double-strand dissociation and formation rates k s s = k d s = 1 . This is a high hitting regime, i.e., the probability of hitting is close to one P ( τ v ( θ ) < ) 1 . In Figure A2, we plot measures of a typical realization of X t on ( H i ) and additionally the probability of reactions over time. High-fidelity replicators ended the simulation with 25% concentration (Figure A2a) and RNA polymerase output ∼62% (Figure A2c), indicating that the concentration of high-fidelity replicators will continue to increase. All species exhibit superexponential growth (Figure A2b). Clay polymerization decreases in contribution over time, whereas RNA polymerization increases substantially over time, and RNA double-strand reactions are small and stable (Figure A2d).

3.3. Hitting Times: Functional and Survival Analysis

We study various models in order of increasing complexity. We examine the hitting time surface τ v ( θ ) in the parameters θ Θ , including probability of hitting P ( τ v ( θ ) < ) . We begin with the core model with no decay or clay.

3.3.1. Core Model, τ v ( θ ) for v = 0.1 with θ = ( n , k ) and “Tent” Functions

We are interested in the structure of the hitting time τ v ( θ ) of (19) as a function of the parameter vector θ . We use the Weibull–Cox proportional hazard’s model of Table 1 for the hitting time τ v for volume fraction v = 0.1 . Let θ = ( n , k ) with sequence dimension n { 3 , 4 } and fitness and similarity parameters k = l = log ( i ) / n for i { 0.1 , 0.05 , 0.01 , 0.005 , 0.001 } . Set m = log ( 0.25 ) / n for fidelity probability parameters. For each value of sequence dimension n, take random initial population X 0 with initial population size I = | X 0 | = 10 and random singleton R = { { x } } for the high-fidelity manifold and fix these for the fitness/similarity landscape curvature parmeters k = l . We fix the double-strand dissociation and formation rate parameters k s s = k d s = 1 and set RNA polymerization rate k r e p such that the overall RNA polymerization rate is given by k ˘ r e p ( 0 ) = 10 and use the “tent” function for fitness, similarity, and fidelity. We take 10 realizations of hitting time τ v ( θ ) for each parameter vector θ Θ and allocate 5000 reactions. This gives 100 independent hitting times and up to 500,000 reactions. The times are comparable because the system is initialized to the same replication mass.
For the simulations, 66 hitting times are finite. The coefficients positively contribute to hitting, where γ n 0.97 and γ k 13.29 , both with p-values less than 0.005. Therefore, hittings are strongly positively influenced by the parameters of the fitness and similarity functions and less so by the dimension. Plots of the coefficients and survival and cumulative hazard curves are given in Figure A3. High survival is found for k large and low survival for k small. Cumulative hazard is highest for k small.
We estimate HDMR of the classifier (whether or not hitting time is finite) using all 100 samples. The results are shown in Figure A5 and Table 4. The HDMR explains roughly 80% of variance. S k 0.69 and S n 0.06 , so hitting probability is strongly influenced by fitness landscape curvature k and less so by sequence length n. The component functions f k and f n for fitness landscape curvature and sequence dimension are strictly decreasing, where larger fitness landscape curvature parameter k results in decreasing hitting probability. These results are consistent with the survival analysis.
We estimate HDMR of the regressor (hitting time) using the 66 simulations with finite hitting time. The results are shown in Figure A4 and Table 5. The HDMR explains roughly 60% of variance. S n 0.57 and S k 0.04 , so sequence dimension dominates the hitting time. Both HDMR component functions f n and f k for sequence dimension and fitness landscape curvature are increasing. The HDMR results reveal that conditioning on hitting reverses the roles of sequence dimension n and fitness landscape curvature k.

3.3.2. Clay and Decay Model, τ v ( θ ) for v = 0.1 with θ = ( n , k , k , k c l a y p , p ) and “Tent” Functions

We expand the model to include clay and decay. We take parameter vector
θ = ( n , k , k , f c l a y , p c l a y ) Θ Θ = { 3 , 4 } × { log ( i ) / n : i = 0.1 , 0.05 , 0.01 , 0.005 , 0.001 } × ( 0 , 1 ) × ( 0 , 1 ) × ( 0 , 1 )
with double-strand dissociation and formation and clay oligomerization reaction rates k s s = k d s = k c l a y o = 1 . For each parameter vector θ Θ , (i) we choose singleton high-fidelity replicator manifold R = { x } for some RNA sequence x E and choose random initial population of RNA molecules X 0 such that the initial population size is 10, I = | X 0 | = 10 , and where the initial population does not intersect the high-fidelity manifold X 0 R = , that is, the initial population does not reside on the high-fidelity manifold; (ii) we initialize the replicative mass of the system such that the initial overall RNA polymerization reaction rate is given by k ˘ r e p ( 0 ) = ( 1 f c l a y ) 20 and the initial overall clay RNA polymerization reaction rate k ˘ c l a y p ( 0 ) = f c l a y 20 ; (iii) we sample the hitting times τ v ( θ ) for volume fraction of high-fidelity replicators v = 0.10 a total of M = 10 times, each censored by 5000 reactions, giving hitting time set T ( θ ) of (20). We attain input–output data set as D = { ( θ i , T ( θ i ) ) : i = 1 , , 240 } . This gives a total of 2400 simulations.
For the simulations, 1546 hitting times are finite. The results of fitting the Weibull–Cox model are shown below in Table 6 and Figure A6. The curvature parameter k again significantly dominates with a large positive value. All the remaining parameters have values less than one. Sequence dimension n again is a relatively small positive contributor. The clay fraction f c l a y is small and positive and replication fidelity parameter p is not significant.
We estimate HDMR of the classifier (whether or not hitting time is finite) using all 2400 samples. Component functions and sensitivity indices are shown below in Figure A7 and Table 7. First-order HDMR captures 74% of explained variance, and second-order captures 4%. Curvature dominates hitting probability with large sensitivity index S k 67 % . The HDMR component function in landscape curvature f k is a decreasing function, where small values increase and large values decrease hitting probability. Sequence dimension n has sensitivity index S n 2 % , and the HDMR component function f n is decreasing, where high dimension decreases the probability of hitting. Clay parameter sensitivity index is small S f c l a y 2 % , and the HDMR component function for fractional clay RNA polymerization rate, f f c l a y , is decreasing, where low-to-medium clay fractions increase and high-clay fractions decrease probability of hitting. The HDMR results are consistent with the Weibull–Cox model.
We estimate HDMR of the regressor (hitting time) using the 1546 simulations with finite hitting time. Component functions and sensitivity indices are shown below in Figure A8 and Table 8. First-order HDMR captures 33% of explained variance, and second-order captures 7%. In stark contrast to the contributions to the classifier, the parameters k and n are insignificant to hitting time. Instead, the largest sensitivity index is S f c l a y 20 % . The HDMR component function for fractional clay RNA polymerization rate, f f c l a y , is an increasing function, where small f c l a y decreases and large f c l a y increases the hitting time. This suggests that high clay-fractions representing first-order reactions increase the hitting time, as clay polymerization has less replicative mass than RNA polymerization, i.e., things go faster with RNA polymerization. The second largest sensitivity index is decay S k 11 % . Decay is an increasing function, with sharp increase in hitting times nearby one, i.e., things go slower with large decay resulting in increased hitting time.

3.4. Compartmentalization

Compartmentalization has a direct effect on the calculation of the reaction rates, specifically replication, by computing only a subset of the reactions in X 1 2 . Put
X t ( A ) = { x X t : l ( x ) A } .
For vesicle region A M , we have that
k ˘ r e p ( t , A ) = ( x , y ) X 1 2 k r e p ( x , y ) M t ( { x } × A ) ( M t ( { y } × A ) I ( x = y ) ) for t R 0 , A [ T , T ] = ( x , y ) X t 2 ( A ) k r e p ( x , y ) M t ( { x } × A ) ( M t ( { y } × A ) I ( x = y ) )
and total replicative mass
k ˘ r e p ( t ) = A M k ˘ r e p ( t , A ) for t R 0
As M increases in size over time, the number of partitions grows, and
A M | X t 2 ( A ) | | X 1 2 | .
Therefore, the replicative mass will be reduced with M , and the system evolves less quickly. This suggests that there is a trade-off between the degree of compartmentalization and the replicative mass of the system.

4. Discussion and Conclusions

Origins of life is a fascinating problem. The wonderful complexity of extant life follows from origins. The distribution of life in the universe is tied to origins.
In this article, we have attempted to peek into the problem by concentrating on the RNA world hypothesis, studying hitting times of high-fidelity replicators. We develop fitness, similarity, and fidelity functions as landscapes for a mathematical model of replicating RNA molecules at the sequence level and observe hitting times through simulation studies. We draw attention to the distinction between the probability of hitting P ( τ ( θ ) < ) and the hitting time τ ( θ ) < .
In terms of mathematical set-up, we interpret the reactions as measure-kernel-functions. Each reaction is identified to and fully encoded by a probability transition kernel. The reactions take place in some domain, whereby all molecules may interact. We note that, in reality, molecules have limited diffusion, and this effectively breaks the reaction domain into independent subdomains above some length scale, i.e., molecules are more likely to react with their neighbors. Therefore, we assume our reaction volume is sufficiently small such that all molecules may participate in the reactions. We use for modeling purposes the ansatz that sequence distance is correlated to spatial proximity, where similar sequences are proximal, using a non-trivial similarity function s : E × E ( 0 , 1 ] .
Theorem 2 and its Corollary 2 show that the system (without decay) is unbounded and strictly increases. This formally shows the system to be a growth process. Next, we illustrate findings about the core system with probable hitting (Section 3.2.1). In particular, we see that the temporal image measure μ t = ν t Q , which describes the polymerization output, is a leading indicator of high-fidelity sequence concentration. Polymerization output and high-fidelity replicators super-exponentially increase near the end of the simulation. Next, the fitness and sequence curvature parameters are set at a higher value which confers reduced fidelity (Section 3.2.2). This reveals that hitting is never achieved and that the polymerization output is in equilibrium with population composition. Hence, the probability of hitting is strongly influenced by landscape curvatures. Next, linear curvature is utilized for fitness and similarity and results in no hitting (Section 3.2.3). This reveals that nonlinear curvature is necessary to achieve hitting of high-fidelity replicators. Next, we expand the model with non-RNA (‘clay’) based polymerization and find that such activity decreases over time, in contrast to RNA polymerization, which greatly increases and dominates other reactions over time the system (Section 3.2.4).
For functional and survival analysis of the hitting times, we study the core model, whereby hitting times are strongly positively influenced by the fitness and similar functions yet are not impacted significantly by sequence dimension (Section 3.3.1). In particular, survival analysis reveals low fitness curvature confers low survival (high hitting), whereas high fitness curvature confers high survival (low hitting). HDMR analysis shows that hitting probability is strongly influenced by fitness curvature and much less so by sequence dimension, supporting the survival analysis. HDMR analysis of hitting time shows reversed roles for sequence dimension and landscape curvature, where sequence dimension dominates hitting time, with curvature playing a far less significant role. This gives the finding that hitting probability is driven by curvature, whereas hitting time is driven by sequence dimension. Next, we perform functional and survival analysis of the core model augmented with ‘clay and decay’ dynamics (Section 3.3.2). Survival analysis shows similar results to the core model, where curvature dominates survival (no hitting), with sequence dimension playing a significantly reduced role. HDMR analysis of hitting probability shows that curvature dominates hitting probability, similar to the core model, whereas sequence dimension again plays a significantly reduced role. HDMR analysis of hitting time reveals that the presence of ‘clay and decay’ significantly increase hitting time, with curvature and sequence dimension playing insignificant roles. These results are consistent in that clay polymerization has less replicative mass than RNA polymerization, where RNA polymerization is a faster reaction.
Overall, we find that nonlinear landscapes are necessary for hitting: linear landscapes are insufficient. For nonlinear landscapes, we find that the probability of hitting is dominated by curvature and that hitting times are dominated by sequence dimension. These results suggest that the landscapes in nature are nonlinear with high curvature, and that the hitting time for high-fidelity replicators is an increasing function of sequence dimension. When clay and decay are added to the model, hitting probability is again dominated by curvature, and clay and decay are relatively insignificant. This reflects that clay and decay are low order reactions. They increase hitting times.
For replication and compartmentalization, we suggest that compartmentalization, while a necessary condition, slows overall system dynamics with increasing vesiculation rate. Essentially, as compartmentalization increases, there is a corresponding reduction in absolute replicative system mass, as certain reactions among elements are no longer possible, being physically sequestered. While the timescale of a simulation is tied to the replicative system mass of the system, there is variability in replicative mass across compartments. Some compartments contain large genomic and metabolic populations. It favors the search for the high-fidelity replicator by there being a distribution on compartmental ‘fitness’ such as resource concentrations so that the high-fitness compartments drive replicative system mass. Compartmentalization is identified to the measure ν on coordinates in ( C , C ) , for which coordinates are “marked” by sequences through the transition probability kernel Q c , and followed by genomic adaptation via the transition probability kernel Q c .
Metabolism is thought to be identified to production of precursors to RNA synthesis, leading to replication identification, followed by genomic adaptation to metabolic state. Metabolism is thus defined through the transition kernels Q m and Q m .
The independence of the transition kernels can be scrutinized, and it is possible that general transition kernels on the full product spaces across location, genomic, and metabolic states are necessary to satisfactorily explain RNA origins, i.e., all three functions may have co-evolved. This notion is suggested in the hot springs hypothesis for origins, where compartmentalization is hypothesized to furnish necessary conditions to genomic and metabolic state [26,27,28]. In this telling, the base measurable space of interest is ( F , F ) = ( C × E × M , C E M ) with measure ν . Then, identification of a high-fidelity replicator is described through the transition probability kernel Q c e m from ( F , F ) into ( F , F ) in “marking” the base measurable space with state for genomic replication; finally, genomic adaptation is conveyed through the kernel Q c e m from ( F × F , F F ) into ( F , F ) . Hence, RNA origins of life has law ν × Q c e m × Q c e m on the product space ( F × F × F , F F F ) , reflecting the steps of replicator identification and adaptation through the definitions transition kernels Q c e m and Q c e m .
A putative “genesis machine” here is a mapping from the base (initial) measurable space ( F , F ) into the product space of identified high-fidelity replicators undergoing adaptation, i.e., ( F × F × F , F F F ) . More generally the base space could additionally contain amino acid sequence space ( P , P ) . Such a machine is fully specified through the definitions of the distribution ν on the base space and the transition kernels Q c e m p and Q c e m p (pre and post genes, respectively). Because the stages of transition occur purely through random drift, an experiment performed by such a machine would take an unacceptably long period of time to complete. Experimental demonstration can be contemplated by augmenting the base measurable space with a control space ( X , X ) to accelerate dynamics, using for instance closed-loop shaped radiation to address molecular degrees of freedom in their appropriate frequency domains, (open-loop) catalysts, temperature, geometry, selection, concentration through centrifugal force, etc., resulting in the new (four-dimensional) base space ( F ˜ , F ˜ ) = ( C × E × M × P × X , C E M P X ) . Then, the transition kernels Q ˜ c e m p and Q ˜ c e m p become mappings from ( F ˜ , F ˜ ) into ( F ˜ , F ˜ ) and from ( F ˜ × F ˜ , F ˜ × F ˜ ) into ( F ˜ , F ˜ ) , respectively. The general design of a genesis machine is the definitions of the augmented base space ( F ˜ , F ˜ ) , its distribution ν ˜ , and the augmented transition kernels Q ˜ c e m p and Q ˜ c e m p , giving law μ ˜ = ν ˜ × Q ˜ c e m p × Q ˜ c e m p on the full 15-dimensional product space ( F ˜ × F ˜ × F ˜ , F ˜ F ˜ F ˜ ) , written in differential notation
μ ˜ ( d x , d y , d z ) = ν ˜ ( d x ) Q ˜ c e m p ( x , d y ) Q ˜ c e m p ( ( x , y ) , d z )
Origins could be experimentally demonstrated using a sequence of adaptive control fields in ( X , X ) , cycling through the transitions, and a detection system for online identification of the system whose elements belong to the product measurable space. The full design and estimated operating timescale for such a machine needs further research to assess practical feasibility. We call the creation of primordial life (primordia) by the continuous causal efforts of a genesis machine given initial prebiotic conditions artebiogenesis, where arte- is Latin and means “from skill.” The primordia are not necessarily those that occurred in nature. Primordia and their genesis represent non-trivial system trajectories across the transition to the earliest life in sterile environments and belong to a manifold of primordial lifeforms, each having characteristic geochemical setting.
The probability measure μ t = ν t Q has additional utility to integrate ‘test’ functions or queries about the system. If we let f E 0 be a fitness function, then the fitness value J ( μ t ) = μ t ( f ) is the expected value of the fitness function with respect to the probability measure μ t . In OptiEvo theory, J ( μ t ) is studied as a function of the population X t on ( E , E ) [52]. OptiEvo assumes that the set of all probability measures { μ t } is convex and that X t has sufficient flexibility such that J ( μ t ) may be explored around μ t . Then, OptiEvo predicts that J ( μ t ) has global maxima on ( E , E ) and that these form a connected level-set of sequences with the same fitness value. Both predictions are consistent with our model. The first prediction is consistent with zero distance in fitness and similarity functions for high-fidelity sequences. The second prediction is consistent with the high-fidelity set being a singleton or a product-space construction. A contention is whether X t has sufficient flexibility in exploring J ( μ t ) around μ t . This has direct bearing on the structure of Q: if X t is inflexible, then Q is constrained to certain subspaces of ( E , E ) , i.e., not all transitions are possible.
In future work, the model could be extended to the space of sequences of lengths up to n, ( E * , E * ) or even the space of sequences of all lengths, where distance and similarity functions would utilize a more general string distance metric, e.g., Levenshtein distance. We note that the size of ( E * , E * ) is not much larger than ( E , E ) . Alternative similarity functions could be explored, such as the trivial case of constant similarity, e.g., s ( x , y ) = 1 for ( x , y ) E × E . A limitation of this article is the restriction to short sequences due to computational efficiency. The numerical size of the sequences is mathematically low-dimensional and does not correspond to actual functions of RNA molecules. Other parameter sets can be explored for example using experimentally derived values for reaction rates, so that the timescales are calibrated. Future research could see the simulation software rewritten for a high-performance computing environment, enabling much longer, e.g., length 10–1,000, sequences to be studied. Polymerase fitness can be made empirical using known RdRp sequences as members of the high-fidelity manifold. Another area of future work could be studying the aforementioned transition kernels Q c , Q c , Q m , and Q m . More general models for polymerization transition kernel based on the structure of the Poisson-binomial distribution could be employed. It would be interesting to study lipid-RNA and metabolism-RNA interactions and equip the system with the ability to append nucleotides to their sequences to form functional genes, such as storing useful information for the replication channel, perhaps a Hammerhead ribozyme to convey rolling-circle amplification. We note that transition kernels here generally lack amino acid state and are pure-RNA. An area to explore is the notion of the transition kernel into the space of high-fidelity replicators to depend on amino acid sequences and then to elaborate the system to contain a primitive translation system and examine various hitting times. Additional reactions can be introduced as operations on pairs of sequences, such as concatenation, and others for sequence splitting, and so on, with corresponding transition kernels enabling RNA networking and recombination dynamics.

Author Contributions

All authors contributed to the study conception and design. Material preparation, prepared software, data generation and analysis were performed by C.D.B. The first draft of the manuscript was written by C.D.B. and all authors commented on previous versions of the manuscript. Grant funding was provided through H.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Science Foundation Grant No. CHE-1763198.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Please see https://github.com/calebbastian/originoflife (accessed on 31 October 2021) for the Python software and example script of usage.

Acknowledgments

The authors thank the late Freeman Dyson for the discussions in 2018 of these ideas at the Institute for Advanced Study in Princeton, New Jersey, as well as anonymous reviewers who gave critical comments that substantially improved the quality of this article.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Discrete Probability Space

We define some concepts related to the space ( E , E ) . The discrete probability measure ν on ( E , E ) is defined by
ν ( A ) = ν I A = x E ν { x } I A ( x ) for A E
where ν { x } is the probability mass at the point x E and
I A ( x ) = 1 if x A 0 otherwise
For the collection of non-negative E -measurable functions E 0 , we have
ν ( f ) = x E ν { x } f ( x ) for f E 0 .

Appendix B. Other Fitness Functions

Another fitness function can be defined using polynomials, such as lines, quadratics, etc, in terms of k N > 0
f k ( x , R ) = 1 H ( x , R ) n + 1 k ( 0 , 1 ] for x E
or
f k ( x , R ) = 1 H ( x , R ) n + 1 k ( 0 , 1 ] for x E
However, another surface is using a sigmoid function. Put
E ( x ) = 1 1 + exp [ x ] for x R
We have fitness for k ( 0 , )
f k ( x , R ) = 1 E H ( x , R ) n 2 k 1 E n 2 k ( 0 , 1 ] for x E

Appendix C. Measure-Kernel-Function

We recall a few facts about transition kernel Q. Q defines a function
Q f ( x ) = Y Q ( x , d y ) f ( y ) for x E
that is in X 0 for every function f Y 0 .
For every probability measure ν t on ( E , E ) and at time t 0 , the quantity μ t = v t Q defines a probability measure on ( Y , Y ) as
μ t ( A ) = X ν t ( d x ) Q ( x , A ) for A Y .
For every probability measure ν t on ( X , X ) and function f Y 0 , we have that
μ t ( f ) = ( ν t Q ) f = X ν t ( d x ) Y Q ( x , d y ) f ( y ) .
Here, the spaces are discrete, i.e.,
ν t ( A ) ν t ( I A ) = x X ν t { x } I A ( x ) for A X
where ν t { x } is the probability mass at the point x X at time t 0 .

Appendix C.1. Reactions as Measure-Kernel-Functions

We index the reaction types on ( Z , Z ) = ( N > 0 , 2 N > 0 ) . Let η t be the probability measure on ( Z , Z ) formed from the normalized reaction rates. Let Q * be the transition kernel from ( Z , Z ) into ( X , X ) . Then, η t Q * = ν t is the distribution on ( X , X ) and μ t = ν t Q is the distribution on ( Y , Y ) . Then, a reaction is the mapping ( Z , Z ) ( X , X ) ( Y , Y ) .

Appendix C.2. Deterministic Model

Consider the core model defined in (2)–(4) with reaction rates k ˘ d s ( t ) , k ˘ s s ( t ) , and k ˘ r e p ( t ) . We are neglecting clay and decay for the moment. All the reactions impact ( E , E ) . For double-strand formation, the input space is ( X , X ) = ( E × E , E E ) and the output space is ( Y , Y ) = ( F , F ) . For double-strand dissociation, the input and output spaces are swapped. For polymerization, ( X , X ) = ( E × , E E ) and ( Y , Y ) = ( E , E ) . Hence, ( E , E ) is positively impacted by polymerization, positively impacted by double-strand dissociation, and negatively impacted by double-strand formation. Put [ x ] = N t ( { { x } } ) and [ x , x c ] = N t ( { { x , x c } } ) . Recall the ν t , Q, and μ t = ν t Q for the reactions, e.g., ν t r e p , ν t d s , Q r e p , etc.
For x E , we have the system of m = 3 2 4 n deterministic nonlinear ordinary differential equations (ODEs) in m variables as mean-field equations
d [ x ] d t = f x = k ˘ r e p ( t ) μ t r e p { { x } } + k ˘ s s ( t ) μ t s s { { x } , { x c } } k ˘ d s ( t ) μ t d s { { x , x c } } = k ˘ r e p ( t ) ( y , z ) E × E ν t r e p { { y } , { z } } Q r e p ( ( y , z ) , { { x } } ) + k s s [ x , x c ] k d s [ x ] [ x c ] = ( y , z ) E × E k r e p ( y , z ) [ y ] ( [ z ] I ( y = z ) ) Q r e p ( ( y , z ) , { { x } } ) + k s s [ x , x c ] k d s [ x ] [ x c ] d [ x , x c ] d t = f x x c = k d s [ x ] [ x c ] k s s [ x , x c ]
The fixed points of f are the equilibria of the system, i.e., f ( x ) = 0 for x R 0 m . The Jacobian of the system is
J = d f 1 d x 1 d f 1 d x m d f m d x 1 d f m d x m
The eigenvalues of the Jacobian reveal the stability of the fixed points. If all the eigenvalues of the Jacobian evaluated at the fixed point have negative real parts, then the fixed point is stable. If none of the eigenvalues are zero and at least one of the eigenvalues has a positive real part, then the fixed point is unstable. If at least one eigenvalue is zero, then the fixed point can be either stable or unstable.

Appendix D. Hitting Cardinality

We index the N reactions of X t with arrival times { T i } in ( R ¯ 0 , B R ¯ 0 ) , where R ¯ = R 0 { } . Define the hitting reaction ϖ ( θ ) N in terms of hitting time τ ( θ ) R ¯ 0 as
ϖ ( θ ) = i = 1 N I [ 0 , τ ( θ ) ] ( T i )
ϖ ( θ ) is right-censored at N reactions. If τ ( θ ) = , then ϖ ( θ ) = N . If τ ( θ ) < , then ϖ ( θ ) < N .

Reaction Cardinality

In this section, we study, instead of the hitting time τ v ( θ ) , the hitting reaction number ϖ v ( θ ) for v = 0.1 . We study the core model with parameter vector θ = ( n , k ) . We uniformly sample sequence length n { 3 , 4 , 5 } and fitness/similarity landscape curvature k = l { 0.0001 , 0.0005 , 0.001 , 0.005 , 0.01 , 0.05 , 0.1 } and form input–output data D = { ( θ i , ϖ ( θ i ) ) : i = 1 , , 1200 } . We set double-strand dissociation and formation rates k s s = k d s = 1 and initialize the overall RNA polymerization rate k r e p such that the initial replicative mass is 10. For each parameter vector θ Θ , we randomly choose the initial population of sequences X 0 with initial population size I = | X 0 | = 10 and singleton high-fidelity sequence manifold R such that the initial population of sequences does not intersect the high-fidelity replicator manifold, X 0 R = .
D contains 781 hitting events. We denote these D * . We form a first-order HDMR on ϖ ( θ ) using D * . The HDMR truth-plot, component functions, and sensitivity indices are shown below in Figure A9. First-order HDMR captures approximately 31% of variance. Both component functions f n and f k have similar sizes, with f n somewhat larger than f k . The HDMR component function for sequence dimension f n is essentially an increasing linear function of sequence length n, and component function for landscape curvature f k is generally an increasing function of the fitness/similarity landscape curvature. These results suggest that sequence dimension and curvature influence the hitting reaction. Larger sequences and flatter curvature increase the hitting reaction.
Table A1. HDMR sensitivity indices of ϖ v ( θ ) < for core model and v = 0.1 .
Table A1. HDMR sensitivity indices of ϖ v ( θ ) < for core model and v = 0.1 .
θ S θ
Sequence dimension n0.1619
Curvature k0.1464
0.3083

Appendix E. Approximate Reaction Rates

One approach to reducing the computational complexity of k ˘ ( t ) is to approximate the sums using Monte Carlo. Define random variables X 1 Uniform ( X 1 ) and X 2 Uniform ( X 2 ) . Let { X 1 i } and { X 2 i } be independencies of such random variables. Given N random samples of X 1 , the first reaction rate becomes
k ˘ d s ( t | N ) = | X 1 | 2 N i = 1 N k d s N t ( X 1 i ) N t ( X 1 i c )
whose expected value is approximated using M realizations,
E k ˘ d s ( t | M , N ) | X 1 | 2 M N j = 1 M i = 1 N k d s N t ( X 1 i j ) N t ( X 1 i j c )
requiring a total of M N evaluations. In a similar manner, the second reaction rate is
E k ˘ s s ( t | M , N ) | X 2 | M N j = 1 M i = 1 N k s s N t ( X 2 i j X 2 i j c )
and, putting X 1 * Uniform ( X 1 ) , we have the third reaction
E k ˘ r e p ( t | M , N ) | X 1 | 2 M N j = 1 M i = 1 N k r e p ( X 1 i j , X 1 i j * ) N t ( X 1 i j ) ( N t ( X 1 i j * ) I ( X 1 i j = X 1 i j * ) )
We refer to SSA simulation with Monte Carlo approximate reaction rates as Monte Carlo Approximate SSA, or MCASSA.

Appendix F. High Dimensional Model Representation

Suppose we have a real-valued square-integrable function f ( x ) L 2 ( E , E , ν ) with E = R n , E = B R n , ν = i ν i , and x = ( x 1 , , x n ) E . The { ν i } may be diffuse (continuous) and/or discrete. Put B = { 1 , , n } . We would like to decompose f into orthogonal function subspaces { f u : u B } (projections) in such a way that each projection on an input subspace f u maximizes variance and across subspaces retrieves total variance, i.e., f = u B f u and V a r f = u B V a r f u . The solution to this problem in the retrieval of { f u : u B } is known as high dimensional model representation (HDMR) or functional ANOVA expansion and for f is written as
f ( x 1 , , x n ) = f 0 + i f i ( x i ) + i < j f i j ( x i , x j ) + + f 1 n ( x 1 , , x n )
where the { f u : u B } are called component functions. For independent inputs, the component functions are mutually orthogonal and, aside from the constant component function f 0 = E f (order zero), have zero mean E f u = 0 for all non-empty ( 2 n 1 ) subspaces, where
V a r f u = E f u 2 ( x u ) ν ( d x ) < for u B
and
V a r f = E ( f ( x ) f 0 ) 2 ν ( d x ) = u B V a r f u < .
A key principle of HDMR is that the expansion for most f may be truncated at low order T n in a T-order HDMR,
f ( x ) f T ( x ) = u B : | u | T f u ( x u ) for T n .
HDMR is often used in global sensitivity analysis to assess input–output correlations at various orders, where the variances are normalized to define sensitivity indices
S u = V a r f u V a r f for u B .
When the inputs are correlated ν i ν i , then the component functions may still be uniquely recovered under hierarchical orthogonality, the variance decomposes
V a r f = u , v B C o v ( f u , f v ) ,
where
C o v ( f u , f v ) = E f u ( x u ) f v ( x v ) ν ( d x ) for u , v B
and the sensitivity indices generalize to structural and correlative sensitivity indices [56], defined respectively as
S u a = V a r f u V a r f for u B
and
S u b = v B : u v C o v ( f u , f v ) V a r f for u B
with total sensitivity index
S u = S u a + S u b for u B
where u B S u = 1 . We use the total sensitivity index as a measure of variable importance and the component functions as profiles of output dependence on the input subspaces.

Appendix G. Reliability Definitions

Given a failure distribution f and reliability (survival) distribution R, we give some relations: the cumulative failure distribution is defined as
F ( t | ϑ ) = 0 t f ( s | ϑ ) d s
where
R ( t | ϑ ) + F ( t | ϑ ) = 1 ,
the hazard rate h ( t | ϑ ) is defined as
h ( t | ϑ ) = f ( t | ϑ ) R ( t | ϑ ) ,
the cumulative hazard is defined as
H ( t | ϑ ) = 0 t h ( s | ϑ ) d s ,
and we have reliability expressed in terms of the cumulative hazard
R ( t | ϑ ) = e H ( t | ϑ ) .
Another useful quantity is the mean residual life
μ ( t | ϑ ) = t R ( s | ϑ ) d s R ( t | ϑ ) .

Appendix H. Additional Figures

Appendix H.1. Linear Landscape

Figure A1. Linear landscape: Measures of system population X t until hitting time τ v for high-fidelity replicator volume fraction v = 0.25 with sequence dimension n = 3 , fitness/similarity curvature k = l = log ( 0.01 ) / n , initial population size, I = | X 0 | = 10 , singleton high-fidelity replicator R = { { x } } , with linear fitness and similarity functions. (a) Concentration of RNA sequences by Hamming distance to high-fidelity replicator; (b) population size of RNA sequences by Hamming distance to high-fidelity replicator; (c) polymerase RNA sequence output by Hamming distance to high-fidelity replicator.
Figure A1. Linear landscape: Measures of system population X t until hitting time τ v for high-fidelity replicator volume fraction v = 0.25 with sequence dimension n = 3 , fitness/similarity curvature k = l = log ( 0.01 ) / n , initial population size, I = | X 0 | = 10 , singleton high-fidelity replicator R = { { x } } , with linear fitness and similarity functions. (a) Concentration of RNA sequences by Hamming distance to high-fidelity replicator; (b) population size of RNA sequences by Hamming distance to high-fidelity replicator; (c) polymerase RNA sequence output by Hamming distance to high-fidelity replicator.
Life 11 01419 g0a1

Appendix H.2. Core Model with Clay

Figure A2. Core model with clay: Measures of system population X t until hitting time τ v ( θ ) for high-fidelity replicator volume fraction v = 0.25 with sequence dimension n = 3 , fitness/similarity curvature k = log ( 0.01 ) / n , initial population size I = | X 0 | = 10 , singleton high-fidelity replicator R = { { x } } , double strand separation and formation rates reaction rate k s s = k d s = 1 , clay replication fidelity probability p = 0.9 , and RNA polymerization rate k r e p and clay polymerization rate k c l a y p chosen such that the replicative mass of each is 10. (a) Concentration of RNA sequences by Hamming distance to high-fidelity replicator; (b) population size of RNA sequences by Hamming distance to high-fidelity replicator; (c) polymerase RNA sequence output by Hamming distance to high-fidelity replicator; (d) probability of reactions over time.
Figure A2. Core model with clay: Measures of system population X t until hitting time τ v ( θ ) for high-fidelity replicator volume fraction v = 0.25 with sequence dimension n = 3 , fitness/similarity curvature k = log ( 0.01 ) / n , initial population size I = | X 0 | = 10 , singleton high-fidelity replicator R = { { x } } , double strand separation and formation rates reaction rate k s s = k d s = 1 , clay replication fidelity probability p = 0.9 , and RNA polymerization rate k r e p and clay polymerization rate k c l a y p chosen such that the replicative mass of each is 10. (a) Concentration of RNA sequences by Hamming distance to high-fidelity replicator; (b) population size of RNA sequences by Hamming distance to high-fidelity replicator; (c) polymerase RNA sequence output by Hamming distance to high-fidelity replicator; (d) probability of reactions over time.
Life 11 01419 g0a2

Appendix H.3. Hitting/Survival Analysis

Figure A3. Survival analysis of hitting time τ v for high-fidelity replicator volume fraction v = 0.1 for core model. (a) Coefficients of the Cox proportional hazard survival model; (b) survival curves in sequence dimension n = L and landscape curvature k = l ; (c) cumulative hazard in sequence dimension n = L and landscape curvature k = l .
Figure A3. Survival analysis of hitting time τ v for high-fidelity replicator volume fraction v = 0.1 for core model. (a) Coefficients of the Cox proportional hazard survival model; (b) survival curves in sequence dimension n = L and landscape curvature k = l ; (c) cumulative hazard in sequence dimension n = L and landscape curvature k = l .
Life 11 01419 g0a3
Figure A4. First-order HDMR analysis of hitting time τ v ( θ ) < for core model. (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence length n = L , f n ( n ) , in sequence length for hitting time, of hitting time; (c) HDMR component function for landscape curvature, f k ( k ) , in landscape curvature, of hitting time. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Figure A4. First-order HDMR analysis of hitting time τ v ( θ ) < for core model. (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence length n = L , f n ( n ) , in sequence length for hitting time, of hitting time; (c) HDMR component function for landscape curvature, f k ( k ) , in landscape curvature, of hitting time. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Life 11 01419 g0a4
Figure A5. First-order HDMR analysis of hitting probability P ( τ v ( θ ) < ) for core model. (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence length, f n ( n ) , in sequence length, of hitting probability; (c) HDMR component function for landscape curvature, f k ( k ) , in landscape curvature, of hitting probability. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Figure A5. First-order HDMR analysis of hitting probability P ( τ v ( θ ) < ) for core model. (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence length, f n ( n ) , in sequence length, of hitting probability; (c) HDMR component function for landscape curvature, f k ( k ) , in landscape curvature, of hitting probability. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Life 11 01419 g0a5
Figure A6. Survival analysis of hitting time τ v for volume fraction v = 0.1 for expanded model (clay and decay). Coefficients of the Cox proportional hazards survival model.
Figure A6. Survival analysis of hitting time τ v for volume fraction v = 0.1 for expanded model (clay and decay). Coefficients of the Cox proportional hazards survival model.
Life 11 01419 g0a6
Figure A7. First-order HDMR analysis of hitting probability P ( τ v ( θ ) < ) for expanded model (clay and decay). (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence length n = L , f n ( n ) , in sequence length, of hitting probability; (c) HDMR component function for curvature, f k ( k ) , in curvature, of hitting probability; (d) HDMR component function for clay fitness, f f c l a y ( f c l a y ) , in clay fitness, of hitting probability. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Figure A7. First-order HDMR analysis of hitting probability P ( τ v ( θ ) < ) for expanded model (clay and decay). (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence length n = L , f n ( n ) , in sequence length, of hitting probability; (c) HDMR component function for curvature, f k ( k ) , in curvature, of hitting probability; (d) HDMR component function for clay fitness, f f c l a y ( f c l a y ) , in clay fitness, of hitting probability. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Life 11 01419 g0a7
Figure A8. First-order HDMR analysis of hitting time τ v ( θ ) < for expanded model (clay and decay). (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence length, f n ( n ) , in sequence length, of hitting time; (c) HDMR component function for curvature, f k ( k ) , in curvature, of hitting time; (d) HDMR component function for clay-fraction, f f c l a y ( f c l a y ) , in clay fraction, of hitting time; (e) HDMR component function for decay rate, f k ( k ) , in decay rate, of hitting time. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Figure A8. First-order HDMR analysis of hitting time τ v ( θ ) < for expanded model (clay and decay). (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence length, f n ( n ) , in sequence length, of hitting time; (c) HDMR component function for curvature, f k ( k ) , in curvature, of hitting time; (d) HDMR component function for clay-fraction, f f c l a y ( f c l a y ) , in clay fraction, of hitting time; (e) HDMR component function for decay rate, f k ( k ) , in decay rate, of hitting time. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Life 11 01419 g0a8
Figure A9. First-order HDMR analysis of hitting cardinality ϖ v ( θ ) < for core model and volume fraction v = 0.1 . (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence dimension, f n ( n ) , in sequence dimension, of hitting cardinality; (c) HDMR component function for curvature, f k ( k ) , in curvature, of hitting cardinality. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Figure A9. First-order HDMR analysis of hitting cardinality ϖ v ( θ ) < for core model and volume fraction v = 0.1 . (a) Hexagonal-bin truth plot; (b) HDMR component function for sequence dimension, f n ( n ) , in sequence dimension, of hitting cardinality; (c) HDMR component function for curvature, f k ( k ) , in curvature, of hitting cardinality. The color function is from blue (negative) to white (zero) to red (positive). The black dots represent standard deviation of the error.
Life 11 01419 g0a9

References

  1. Ganti, T. The Principles of Life; Oxford University Press: Oxford, UK, 2003. [Google Scholar]
  2. Ganti, T. Chemoton Theory; Kluwer Academic/Plenum Publishers: Dordrecht, The Netherlands, 2003. [Google Scholar]
  3. Orgel, L.E. Evolution of the genetic apparatus. J. Mol. Biol. 1968, 38, 381–393. [Google Scholar] [CrossRef]
  4. Gilbert, W. Origin of life: The RNA world. Nature 1986, 319, 618. [Google Scholar] [CrossRef]
  5. Joyce, G.F. The antiquity of RNA-based evolution. Nature 2002, 418, 214–221. [Google Scholar] [CrossRef] [PubMed]
  6. White, H.B. Coenzymes as fossils of an earlier metabolic state. J. Mol. Evol. 1976, 7, 101–104. [Google Scholar] [CrossRef]
  7. Eigen, M. Selforganization of matter and the evolution of biological macromolecules. Naturwissenschaften 1971, 58, 465–523. [Google Scholar] [CrossRef] [PubMed]
  8. Eigen, M.; Schuster, P. A principle of natural self-organization. Naturwissenschaften 1977, 64, 541–565. [Google Scholar] [CrossRef]
  9. Cech, T.R. The Ribosome Is a Ribozyme. Science 2000, 289, 878. [Google Scholar] [CrossRef]
  10. Diener, T.O. Potato spindle tuber “virus”. IV. A replicating, low molecular weight RNA. Virology 1971, 45, 411–428. [Google Scholar] [CrossRef]
  11. Tupper, A.S.; Higgs, P.G. Rolling-circle and strand-displacement mechanisms for non-enzymatic RNA replication at the time of the origin of life. J. Theor. Biol. 2021, 527, 110822. [Google Scholar] [CrossRef]
  12. Vaidya, N.; Manapat, M.L.; Chen, I.A.; Xulvi-Brunet, R.; Hayden, E.J.; Lehman, N. Spontaneous network formation among cooperative RNA replicators. Nature 2012, 491, 72–77. [Google Scholar] [CrossRef] [PubMed]
  13. de Farias, S.T.; dos Santos Junior, A.P.; Rêgo, T.G.; José, M.V. Origin and Evolution of RNA-Dependent RNA Polymerase. Front. Genet. 2017, 8, 125. [Google Scholar] [CrossRef] [PubMed]
  14. Koonin, E.V.; Krupovic, M.; Ishino, S.; Ishino, Y. The replication machinery of LUCA: Common origin of DNA replication and transcription. BMC Biol. 2020, 18, 61. [Google Scholar] [CrossRef]
  15. Ghadessy, F.J.; Ong, J.L.; Holliger, P. Directed evolution of polymerase function by compartmentalized self-replication. Proc. Natl. Acad. Sci. USA 2001, 98, 4552. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Tjhung, K.F.; Shokhirev, M.N.; Horning, D.P.; Joyce, G.F. An RNA polymerase ribozyme that synthesizes its own ancestor. Proc. Natl. Acad. Sci. USA 2020, 117, 2906. [Google Scholar] [CrossRef]
  17. Ertem, G.; Ferris, J.P. Template-Directed Synthesis Using the Heterogeneous Templates Produced by Montmorillonite Catalysis. A Possible Bridge Between the Prebiotic and RNA Worlds. J. Am. Chem. Soc. 1997, 119, 7197–7201. [Google Scholar] [CrossRef]
  18. Acevedo, O.L.; Orgel, L.E. Non-enzymatic transcription of an oligodeoxynucleotide 14 residues long. J. Mol. Biol. 1987, 197, 187–193. [Google Scholar] [CrossRef]
  19. Szostak, J.W. The eightfold path to non-enzymatic RNA replication. J. Syst. Chem. 2012, 3, 2. [Google Scholar] [CrossRef] [Green Version]
  20. Cairns-Smith, A.G. Genetic Takeover and the Mineral Origins of Life; Cambridge University Press: Cambridge, UK, 1987. [Google Scholar]
  21. Dyson, F. Origins of Life; Cambridge University Press: Cambridge, UK, 1999. [Google Scholar] [CrossRef] [Green Version]
  22. Sakuma, Y.; Imai, M. From vesicles to protocells: The roles of amphiphilic molecules. Life 2015, 5, 651–675. [Google Scholar] [CrossRef]
  23. Szostak, J.W.; Bartel, D.P.; Luisi, P.L. Synthesizing life. Nature 2001, 409, 387–390. [Google Scholar] [CrossRef]
  24. Segré, D.; Ben-Eli, D.; Deamer, D.W.; Lancet, D. The Lipid World. Orig. Life Evol. Biosph. 2001, 31, 119–145. [Google Scholar] [CrossRef]
  25. Martin, W.; Baross, J.; Kelley, D.; Russell, M.J. Hydrothermal vents and the origin of life. Nat. Rev. Microbiol. 2008, 6, 805–814. [Google Scholar] [CrossRef] [PubMed]
  26. Damer, B.; Deamer, D. The Hot Spring Hypothesis for an Origin of Life. Astrobiology 2019, 20, 429–452. [Google Scholar] [CrossRef] [Green Version]
  27. Damer, B.; Deamer, D. Coupled phases and combinatorial selection in fluctuating hydrothermal pools: A scenario to guide experimental approaches to the origin of cellular life. Life 2015, 5, 872–887. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Deamer, D.; Damer, B.; Kompanichenko, V. Hydrothermal Chemistry and the Origin of Cellular Life. Astrobiology 2019, 19, 1523–1537. [Google Scholar] [CrossRef]
  29. Kvenvolden, K.; Lawless, J.; Pering, K.; Peterson, E.; Flores, J.; Ponnamperuma, C.; Kaplan, I.R.; Moore, C. Evidence for Extraterrestrial Amino-acids and Hydrocarbons in the Murchison Meteorite. Nature 1970, 228, 923–926. [Google Scholar] [CrossRef] [PubMed]
  30. Miller, S.L.; Urey, H.C. Organic Compound Synthesis on the Primitive Earth. Science 1959, 130, 245. [Google Scholar] [CrossRef]
  31. Pino, S.; Sponer, J.E.; Costanzo, G.; Saladino, R.; Mauro, E.D. From formamide to RNA, the path is tenuous but continuous. Life 2015, 5, 372–384. [Google Scholar] [CrossRef] [Green Version]
  32. Powner, M.W.; Sutherland, J.D. Prebiotic chemistry: A new modus operandi. Phil. Trans. R. Soc. B Biol. Sci. 2011, 366, 2870–2877. [Google Scholar] [CrossRef] [Green Version]
  33. Powner, M.W.; Gerland, B.; Sutherland, J.D. Synthesis of activated pyrimidine ribonucleotides in prebiotically plausible conditions. Nature 2009, 459, 239–242. [Google Scholar] [CrossRef]
  34. Attwater, J.; Wochner, A.; Holliger, P. In-ice evolution of RNA polymerase ribozyme activity. Nat. Chem. 2013, 5, 1011–1018. [Google Scholar] [CrossRef] [Green Version]
  35. Hays, L. NASA Astrobiology Strategy; Technical report; National Aeronautics and Space Administration: Washington, DC, USA, 2015.
  36. Coveney, P.V.; Swadling, J.B.; Wattis, J.A.D.; Greenwell, H.C. Theory, modelling and simulation in origins of life studies. Chem. Soc. Rev. 2012, 41, 5430–5446. [Google Scholar] [CrossRef]
  37. Lanier, K.A.; Williams, L.D. The Origin of Life: Models and Data. J. Mol. Evol. 2017, 84, 85–92. [Google Scholar] [CrossRef] [Green Version]
  38. Wu, M.; Higgs, P.G. The origin of life is a spatially localized stochastic transition. Biol. Direct 2012, 7, 42. [Google Scholar] [CrossRef] [Green Version]
  39. Gillespie, D.T. Exact stochastic simulation of coupled chemical reactions. J. Phys. Chem. 1977, 81, 2340–2361. [Google Scholar] [CrossRef]
  40. Hordijk, W.; Steel, M. Detecting autocatalytic, self-sustaining sets in chemical reaction systems. J. Theor. Biol. 2004, 227, 451–461. [Google Scholar] [CrossRef] [Green Version]
  41. Walker, S.I. Origins of life: A problem for physics, a key issues review. Rep. Prog. Phys. 2017, 80, 092601. [Google Scholar] [CrossRef]
  42. Wattis, J.A.D.; Coveney, P.V. The Origin of the RNA World: A Kinetic Model. J. Phys. Chem. B 1999, 103, 4231–4250. [Google Scholar] [CrossRef] [Green Version]
  43. Kun, Á.; Szilágyi, A.; Könnyű, B.; Boza, G.; Zachar, I.; Szathmáry, E. The dynamics of the RNA world: Insights and challenges. Ann. N. Y. Acad. Sci. 2015, 1341, 75–95. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  44. Szilágyi, A.; Zachar, I.; Scheuring, I.; Kun, Á.; Könnyű, B.; Czárán, T. Ecology and Evolution in the RNA World Dynamics and Stability of Prebiotic Replicator Systems. Life 2017, 7, 48. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  45. Scheuring, I.; Szilágyi, A. Diversity, stability, and evolvability in models of early evolution. Curr. Opin. Syst. Biol. 2019, 13, 115–121. [Google Scholar] [CrossRef] [Green Version]
  46. Takeuchi, N.; Hogeweg, P. Evolutionary dynamics of RNA-like replicator systems: A bioinformatic approach to the origin of life. Phys. Life Rev. 2012, 9, 219–263. [Google Scholar] [CrossRef] [Green Version]
  47. Orgel, L. A Simpler Nucleic Acid. Science 2000, 290, 1306. [Google Scholar] [CrossRef] [PubMed]
  48. Maury, C.P.J. Origin of life. Primordial genetics: Information transfer in a pre-RNA world based on self-replicating beta-sheet amyloid conformers. J. Theor. Biol. 2015, 382, 292–297. [Google Scholar] [CrossRef] [Green Version]
  49. Ehrenfreund, P.; Rasmussen, S.; Cleaves, J.; Chen, L. Experimentally Tracing the Key Steps in the Origin of Life: The Aromatic World. Astrobiology 2006, 6, 490–520. [Google Scholar] [CrossRef]
  50. Kunin, V. A System of Two Polymerases—A Model for the Origin of Life. Orig. Life Evol. Biosph. 2000, 30, 459–466. [Google Scholar] [CrossRef] [PubMed]
  51. Wright, S. Evolution and the Genetics of Populations, Volume 1; The University of Chicago Press: Chicago, IL, USA, 1984. [Google Scholar]
  52. Feng, X.; Pechen, A.; Jha, A.; Wu, R.; Rabitz, H. Global optimality of fitness landscapes in evolution. Chem. Sci. 2012, 3, 900–906. [Google Scholar] [CrossRef]
  53. Davidson-Pilon, C. lifelines: Survival analysis in Python. J. Open Source Softw. 2019, 4, 1317. [Google Scholar] [CrossRef] [Green Version]
  54. Bornholt, J.; Lopez, R.; Carmean, D.M.; Ceze, L.; Seelig, G.; Strauss, K. A DNA-Based Archival Storage System. In Proceedings of the Twenty-First International Conference on Architectural Support for Programming Languages and Operating Systems, Association for Computing Machinery, New York, NY, USA, 19–23 April 2016; pp. 637–649. [Google Scholar] [CrossRef] [Green Version]
  55. Kitadai, N.; Maruyama, S. Origins of building blocks of life: A review. Geosci. Front. 2018, 9, 1117–1153. [Google Scholar] [CrossRef]
  56. Li, G.; Rabitz, H. General formulation of HDMR component functions with independent and correlated variables. J. Math. Chem. 2012, 50, 99–130. [Google Scholar] [CrossRef]
Figure 1. Measures of system population X t until hitting time τ v for high-fidelity replicator volume fraction v = 0.25 with sequence dimension n = 3 , fitness/similarity curvature l = k = log ( 0.01 ) / n , initial population size I = | X 0 | = 10 , singleton high-fidelity replicator R = { { x } } , with “tent” fitness and similarity functions. (a) concentration of RNA sequences by Hamming distance to high-fidelity replicator; (b) population size of RNA sequences by Hamming distance to high-fidelity replicator; (c) polymerase RNA sequence output by Hamming distance to high-fidelity replicator.
Figure 1. Measures of system population X t until hitting time τ v for high-fidelity replicator volume fraction v = 0.25 with sequence dimension n = 3 , fitness/similarity curvature l = k = log ( 0.01 ) / n , initial population size I = | X 0 | = 10 , singleton high-fidelity replicator R = { { x } } , with “tent” fitness and similarity functions. (a) concentration of RNA sequences by Hamming distance to high-fidelity replicator; (b) population size of RNA sequences by Hamming distance to high-fidelity replicator; (c) polymerase RNA sequence output by Hamming distance to high-fidelity replicator.
Life 11 01419 g001
Figure 2. Measures of system population X t until hitting time τ v for high-fidelity replicator volume fraction v = 0.25 with sequence dimension n = 3 , fitness/similarity curvature l = k = log ( 0.1 ) / n , initial population size I = | X 0 | = 10 , singleton high-fidelity replicator R = { { x } } , with “tent” fitness and similarity functions. (a) concentration of RNA sequences by Hamming distance to high-fidelity replicator; (b) population size of RNA sequences by Hamming distance to high-fidelity replicator; (c) polymerase RNA sequence output by Hamming distance to high-fidelity replicator.
Figure 2. Measures of system population X t until hitting time τ v for high-fidelity replicator volume fraction v = 0.25 with sequence dimension n = 3 , fitness/similarity curvature l = k = log ( 0.1 ) / n , initial population size I = | X 0 | = 10 , singleton high-fidelity replicator R = { { x } } , with “tent” fitness and similarity functions. (a) concentration of RNA sequences by Hamming distance to high-fidelity replicator; (b) population size of RNA sequences by Hamming distance to high-fidelity replicator; (c) polymerase RNA sequence output by Hamming distance to high-fidelity replicator.
Life 11 01419 g002
Table 1. Weibull reliability model.
Table 1. Weibull reliability model.
NameQuantityBaselineQuantityProportional Hazards
Failure density f ( t | α , β ) α t ( t β ) α e ( t β ) α f ( t , θ | α , β , γ ) α t ( t β ) α e γ · θ e ( t β ) α e γ · θ
Failure distribution F ( t | α , β ) 1 e ( t β ) α F ( t , θ | α , β , γ ) 1 e ( t β ) α e γ · θ
Reliability distribution R ( t | α , β ) e ( t β ) α R ( t , θ | α , β , γ ) e ( t β ) α e γ · θ
Cumulative hazard H ( t | α , β ) ( t β ) α H ( t , θ | α , β , γ ) ( t β ) α e γ · θ
Hazard rate h ( t | α , β ) α t ( t β ) α h ( t , θ | α , β , γ ) α t ( t β ) α e γ · θ
Table 2. Reactions and orders.
Table 2. Reactions and orders.
ReactionOrder
RNA double strand formation2
RNA double strand dissociation1
RNA polymerization2
RNA decay1
Clay polymerization1
Clay oligomerization0
Compartmentalization1
Metabolism to replication1
Metabolism & replication to metabolism2
Table 3. Model parameters.
Table 3. Model parameters.
θ NameDomainValue(s)
nsequence dimension N > 0 { 3 , 4 , 5 }
kRNA fitness parameter R 0 { log ( i ) / n : i = 0.1 , 0.05 , 0.01 , 0.005 , 0.001 }
lRNA similarity parameter R 0 { log ( i ) / n : i = 0.1 , 0.05 , 0.01 , 0.005 , 0.001 }
mRNA fidelity parameter R 0 log ( 0.25 ) / n
pclay fidelity probability ( 0 , 1 ] 0.9
k s s double-strand dissociation rate R 0 1
k d s double-strand formation rate R 0 1
k r e p RNA replication rate R 0 10
k RNA decay rate R 0 (0, 1)
k c l a y o clay RNA oligomerization rate R 0 1
k c l a y p clay RNA polymerization rate R 0 (0, 20)
Table 4. HDMR sensitivity indices of hitting probability P ( τ v ( θ ) < ) for the core model.
Table 4. HDMR sensitivity indices of hitting probability P ( τ v ( θ ) < ) for the core model.
θ S θ
Sequence length n0.06
Curvature k0.69
0.75
Table 5. HDMR sensitivity indices of hitting time τ v ( θ ) for the core model.
Table 5. HDMR sensitivity indices of hitting time τ v ( θ ) for the core model.
θ S θ
Sequence length n0.57
Curvature k0.04
0.61
Table 6. Weibull–Cox model parameters for hitting times of clay and decay model.
Table 6. Weibull–Cox model parameters for hitting times of clay and decay model.
θ NameCoefficient γ θ p-Value
nsequence dimension0.54<0.005
k = l RNA fitness parameter27.87<0.005
pclay fidelity probability−0.080.35
k RNA decay rate0.80<0.005
f c l a y fraction clay RNA polymerization rate0.96<0.005
Table 7. HDMR sensitivity indices of hitting probability P ( τ v ( θ ) < ) for expanded model (clay and decay).
Table 7. HDMR sensitivity indices of hitting probability P ( τ v ( θ ) < ) for expanded model (clay and decay).
θ S θ
Sequence length n0.0213
Curvature k0.6732
Decay rate k 0.0120
Clay fidelity p0.0114
Fraction clay RNA polymerization rate f c l a y 0.0219
0.7399
Table 8. HDMR sensitivity indices of τ v ( θ ) < for expanded model (clay and decay).
Table 8. HDMR sensitivity indices of τ v ( θ ) < for expanded model (clay and decay).
θ S θ
Sequence dimension n0.0013
Curvature k0.0030
Decay k 0.1129
Clay fidelity p0.0152
Fraction clay RNA polymerization rate f c l a y 0.2014
0.3339
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bastian, C.D.; Rabitz, H. Hitting Times of Some Critical Events in RNA Origins of Life. Life 2021, 11, 1419. https://doi.org/10.3390/life11121419

AMA Style

Bastian CD, Rabitz H. Hitting Times of Some Critical Events in RNA Origins of Life. Life. 2021; 11(12):1419. https://doi.org/10.3390/life11121419

Chicago/Turabian Style

Bastian, Caleb Deen, and Hershel Rabitz. 2021. "Hitting Times of Some Critical Events in RNA Origins of Life" Life 11, no. 12: 1419. https://doi.org/10.3390/life11121419

APA Style

Bastian, C. D., & Rabitz, H. (2021). Hitting Times of Some Critical Events in RNA Origins of Life. Life, 11(12), 1419. https://doi.org/10.3390/life11121419

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop