Nonequilibrium Entropic Bounds for Darwinian Replicators

Life evolved on our planet by means of a combination of Darwinian selection and innovations leading to higher levels of complexity. The emergence and selection of replicating entities is a central problem in prebiotic evolution. Theoretical models have shown how populations of different types of replicating entities exclude or coexist with other classes of replicators. Models are typically kinetic, based on standard replicator equations. On the other hand, the presence of thermodynamical constraints for these systems remain an open question. This is largely due to the lack of a general theory of statistical methods for systems far from equilibrium. Nonetheless, a first approach to this problem has been put forward in a series of novel developements falling under the rubric of the extended second law of thermodynamics. The work presented here is twofold: firstly, we review this theoretical framework and provide a brief description of the three fundamental replicator types in prebiotic evolution: parabolic, malthusian and hyperbolic. Secondly, we employ these previously mentioned techinques to explore how replicators are constrained by thermodynamics. Finally, we comment and discuss where further research should be focused on.


Introduction
Biology follows the laws of physics, and yet it remains distinctive from many standard physical systems in a number of ways. In the first place, life's self-replicating mechanisms stand as a major difficulty when approaching it from a simple physical setup. On the other hand, life too differs from physics in its computational nature: all living forms conduct some sort of computation as a crucial component of their adaptive potential [1]. The success of life over chemistry is largely associated to the emergence of prebiotic molecular mechanisms that, in turn, allowed for a template-based landscape to become dominant over the whole biosphere. How this took place is one of the most fundamental questions in science [2][3][4].
Life forms are out-of-equilibrium structures capable to employ available matter, energy and information to propagate some type of identity. Most theoretical approaches to the evolution of replicators have been grounded on a kinetic description. Under such framework, interactions between (typically molecular) agents are represented by nonlinear differential equations, known as replicator equations [5]. They provide a deterministic view of Darwinian dynamics. However, as pointed out by Smith and Morowitz, "the abstraction of the replicator, which reduces Darwinian dynamics to its essentials, also de-emphasizes the chemical nature of life" [6]. The same can be concluded in relation with the lack of a thermodynamical context. Despite early efforts towards the development of a physics of evolutionary dynamics [6][7][8][9] a more satisfactory formalism has yet to emerge. In particular, life propagation processes require an entropy production and balance equations can be defined [9][10][11]. However, a more general non-equilibrium statistical physics approach suitable for the problem of self-replication has been missing until recently [12][13][14][15][16][17][18][19]. How can this novel approach apply to the fundamental problem of replicator dynamics in the eary stages of Life on Earth? Beyond the self-replicating potential of cells and molecules, several replication strategies are at work in living systems, also involving multiple scales [20][21][22]. The basic growth dynamics followed by each class has remarkably different consequences for selection. The simplest class is the Mathusian (exponential) growth dynamics exhibited by cellular systems growing under unlimited resources. Two other types of replicators are observed in Nature. One is associated to the emergence of cooperation dynamics, with different classes of replicators helping each other and forming a mutualistic assembly [23]. The second is related to a template-based replication mechanism that we can identify in living systems as the standard mechanism of nucleic acid replication. This mechanism has been shown to lead to the "survival of everyone": it provides a mechanism capable of sustaining very diverse populations of replicators [24][25][26].
From the physics perspective, these systems involve large number of internal degrees of freedom interacting in an out-of-equilibrium context. In turn, this interplay in the microscopic level leads to a macroscopic emergent (coarse-grained) dynamics. A thermodynamical connection between these two levels can be made following the statistical physics methods cited above. The work presented here is an attempt to delineate these fundamental thermodynamical constraints for the three elemental types of prebiotic replicators.

Entropic Bounds for Replicators
Let us begin by reviewing the theoretical framework upon which the analysis of the problem will unfold [15,[17][18][19]. Here, we outline a simplified version of theoretical basis behind this non-equilibrium approach. We also comment on the generalizations of the so-called extended second law [19]. Then, we summarize the elemental classes of replicators and their essential aspects [21], together with a series of implications regarding selection and adaptation. Finally, we lay out an approach to the question of how non-equilibrium thermodynamical bounds arise in these types of systems and how such constraints might have affected early evolutionary scenarios.

The Extended Second Law
Consider a classical time-evolving system described by its microscopical trajectory in the phase space x(t) ∈ Ω plus a set of controlled parameters λ(t) evolving in a time interval t ∈ [0, τ] that act like external drivers for any given trajectory. Assume that the system remains in contact with a heat bath at temperature T = 1/β throughout the entire trajectory. Denote the transition probability from a miscroscopical state x to y in the time interval by π [x → y]. Now, if we slice time as t i+1 − t i = , with t n = τ = n and t 0 = 0, then, for sufficiently small , the microscopical reversibility condition implies [13,14]: where the superscript * denotes momentum-reversed microstates, and Q b i→i+1 denotes the heat exchange in going from from states x(t i ) to x(t i+1 ) as measured from the heat bath. Heuristically, (1) is interpreted as the composed detailed balance condition on each time-slice of the trajectory x(t) (see Figure 1a). This can be represented by the functional relation: Next, let us introduce two macrostates which can be interpreted as two disjoint sections of the phase space, A, B ⊂ Ω (see Figure 1b). Let us introduce notation for macrostate bounded trajectories in Ω by defining the set of forward trajectories , the set of possible trajectories subject to condition that the initial microstate is in A and the final must be in B. Then, construct the formal coarse-grained transition rate from A to B as while, equivalently, denote x * τ = {x * (τ − t), t ∈ [0, τ] | x * (τ) ∈ B ∧ x * (0) ∈ A} as the set of reversed macrostate bounded trajectories, driven by the reverse protocolλ(τ − t) (details on the derivation can be found in [19]), and compute the inverse coarse-grained transition rate from B to A as Here onwards, let use bracket notation · to denote averages over forward paths x τ . Under this theoretical framework, it can be shown [17,19] that the following relation must hold: where we have defined the path-dependant observable: with p τ (x(τ)) and p 0 (x(0)) standing for the probability of landing at a certain x(τ) ∈ B at time t = τ and departing from x(0) ∈ A at time t = 0. Notice that (6) is a functional that depends on the boundary conditions of the trajectory x(t). Let us define, as a functional observable over the sample of forward paths x τ . On the one hand, a first order expansion on (5) imposes the following boundaries to the fraction of the coarse-grained transition rates: This results implicitely allude to the Landauer bounds on heat production for bit erasure [27][28][29]. Inequality (8) constraints the irreversibility of the macroscopic process A → B with respect to the average generalized entropy produced internally, ∆H, and externally (into the bath), βQ b , and it is dubbed the Extended/Bayesian Second Law (ESL) [17,19]. One interpretation is that macroscopic irreversibility increases the minimum dissipated energy during the process A → B. Interestingly, expression (8) formalizes a bound on entropy production in relation to the coarse-grained properties of the process, such as the macroscipic transition rates. This result is of particular interest since, under many experimental circumstances, these are the only measurable quantities for a given system. We will come back to this point in the following sections.
On the other hand, a general perturbative analysis using the cumulant expansion [30] onto (5) leads to where ω l stands for the l−th cumulant of the distribution of βW [x(t)]. In fact, (9) allows for a more sophisticated view of where, formally with the subscript c indicating cumulant expressions. Combining Equations (5) and (9), it can be shown that Φ τ ≥ 0. Indeed, Φ τ is a measure the fluctuations of the distribution associated to observable W [x(t)]. Thus, Equation (10) represents an extended fluctuation-dissipation theorem, where the LHS reflects the macroscopic (coarse-grained) irreversibility property and the RHS a balance between dissipated work and fluctuations over the x τ sample. This result is of particular interest when a system is arranged such that a choice between two macroscopical end-states is forced. In such cases, fluctuation discrepancies might break symmetry thus favoring certain macroscopical transitions or supressing others [18].
Moreover, these theoretical results can be generalized to less constrained versions of the ESL where no equlibrium trajectory end-points are required plus the system needs not to be at a fixed temperature, eventhough there is still contact with a heat bath (cf. [19]). Under this generalized lens, relations (5), (8)-(10) are formally equivalent, only now the space of possibilities over which averages are taken is constrained by the implemented coarse-grain. On the other hand, this implies that the operators in (7) are too redefined owing to the coarse-graining imposed in the system.
In the following sections we will revisit the paradigm of prebiotic replicators, and focus on how to minimally embed this problem into the formalism discussed above. Subsequently, we will argue how these entropic constraints may have coupled to prebiotic selection and added preassure to in an evolutionary context.

Replicators & Reproducers
Several fundamental replication strategies are at play in living systems. These strategies are present in multiple scales, from molecular replicators to cells and beyond. Each class of replicating agent is characterized by a kinetic pattern, which dynamics entail distinct selective implications.
Here, we will focus on three characteristic replicator classes [20,21].
Simple replicators: commonly known as Malthusian agents, correspond to systems whereby a single component A is capable of making a copy of itself by using the available resources, namely E, generating a certain waste product, W. Schematically, Assuming a large repository of resources, the kinetics of this process can be reduced to a linear dynamical equation (see Table 1). Systems following this mechanism obey exponential growth laws.
Hyperbolic replicators: one of the most relevant novelties in evolution [31,32] is the concept of autocatalysis. This mechanism is a precursor of self-replicating entities that largely define the nature of living structures. It has been put forward by several authors [33][34][35][36] as a central process in the chemistry of prebiotic systems involving the emergence of cooperative agents (see Figure 2a).
Again, under well-mixed and unlimited resource conditions, the hyperbolic replicator kinetics is reduced to a second order equation (see Table 1). Autocatalytic growth is characterized by displaying a finite-time singularity at t c = 1/hx 0 [21]. a b Figure 2. Hyperbolic and parabolic replicators. In (a) we display a simplified scheme of an experimental implementation of a catalytic set of ribozymes forming a cooperative loop. Here each component of the system helps the next to replicate. Dashed lines indicate weaker catalytic links (modified from [37]).
The parabolic system outlined in (b) is based on complementary (template) peptide chains involving a ligation mechanism (adapted from [38]).
Parabolic replicators: this type of replicator arises from a combination of molecular reactions. In particular, oligonucleotides are known to exhibit such behaviour [26,[39][40][41]. The minimal scheme where this particular dynamics is observed consists of the set of processes (see Figure 2b).
which, under conditions a b c is reduced to a parabolic lawẋ = ρ √ x, where x denotes the total concentration of the molecular component A regardless of the configuration, it being either associated (AA) or dissociated (A) (see Appendix A). Parameter ρ = c √ 2b/a. Table 1. Summary of the minimal expressions for the kinetics of the three replicator classes discussed above. We have denoted as x the gross concentration of replicating molecules A, independently of the configuration.

Replicator Class Reaction Scheme Effective Dynamics
Simple

Coarse-Grained Dynamics of Replicators
The dynamics of the three types of replicators discussed above are taking place on the macroscopic level. Molecular replicators encapsulate a whole system rich in complexity and structure, thus the measurable transition rates, such as g, h or ρ above, are emergent features of the interplay of the many internal degrees of freedom of the system. However, the statistical properties of these phenomena are non-ergodic, since replicating is constrained by an initial and a final coarse-grained states. As discussed in Section 2.1, averages reflecting the macroscopic transition rates are taken over a section of the space of possibilities, specifically over the subset of possible microscopical trajectories with an initial number of replicators n − 1 and a final number n (given a time scale τ), as detailed below.
To begin with, suppose that a system is composed of a fixed number of molecular templates or chains, N, which can either be internally ordered such that they behave as a replicators (A), namely active chains, or simply act as substrate (E), namely inactive chains. The goal here is to define an unambiguous coarse-graining measure capable of distinguishing two meaningful macroscopic states of the system. To do so, we will consider three such systems which replicators' act accordingly with the three replicator classes summarized in Table 1. We will also suppose that all replicators undergo equivalent decay processes. This assumption is taken so that we are able to probe the thermodynamical bounds purely for the processes involving replication. For simplicity, we use open systems (source flowing in) but finite (fixed total number of particles).
Following a markovian approach [42,43], each set of reaction rules allows defining transition probabilities and a master equation that in general will read: which gives the probability P(n, t) of observing n active chains at time t. Here the ω(i|j) terms introduce the transition probabilities associated to each rule, duely determined by the corresponding Malthusian, hyperbolic and parabolic cases. The three urn-like systems analysied here are chemostat models since, when an element (replicator) decays, it is replaced by newly available source particles E (see Appendix B for details). In summary, dP(n, t) dt = bc 2a Notice that (16)-(18) are non-equilibrium macroscopic representations of the replicating dynamics. Here, the internal interactions that produce the effective behaviour described by the previous set of equations are all integrated out into its corresponding coupling constants. Thus, within this macroscopical framework we shall define the phase space subsets: • A-state in which the system contains a total number of n − 1 active chains. • B-state in which the system contains a total amount of n active chains.
Let us focus on the explicit bounds given by the LHS in expression (5). We first introduce notation for these lower entropic bounds, where the subscript r ∈ {s, h, p} indicates the replicator type (simple, hyperbolic and parabolic respectively), while x := n/N in each case. Therefore, considering that the transition rates Π τ (A → B) and Π τ (B → A) for the defined coarse-grained states A and B correspond to the prefactors in each master equation above, where we have defined α := b/2a. Finally, introduce notation ∆LEB(r|r ) := LEB r (x) − LEB r (x) in order to compare each replicator type. Hence, for h and p against s we derive while, ∆LEB(h|p) = ∆LEB(h|s) − ∆LEB(p|s). Notice that, since all replicators decay mechanism has been chosen to be equivalent (see Appendix B), then relative bounds ∆LEB(r|r ) are δ−independent. Figure 3a-f show various curves (22) and (23) against the density value x.
Focusing on the limiting cases where the lower bounds between distinct replicators coincide, ∆LEB(r|r ) = 0, it is possible to derive the density values for which the LEB for replicator r exceeds that of replicator r and viceversa. This is an interesting exercise since minimal entropy production can provide a guideline for thermodynamically advantageous processes. Bare in mind that exploring LEBs does not include the full picture, as fluctuations can shift the average dissipared energy and unbalance irreversibility as discussed above (cf. [18]).
Thus, let us define the LEB crossover density x rr from r-LEB dominance to r -LEB dominance, or, simply, ∆LEB(r|r ) x rr = 0. Working with reduced variablesh := h/g andc := c/g we derive x rr = x rr (h,c) following (22) and (23): where the equation for x ph , the density value where LEB dominance shifts from parabolic hyperbolic is given in an implicit form (Algebraic analysis shows that the equation for x ph contains a single real root.). On the other hand, 0 < x rr < 1 must be held, as it stands for a density variable. These considerations allow for a construction of a diagram (h,c) where space is divided into sections characterised by the replicator-types that display a dominant LEB. For instance, forh,c < 1 the simple replicator's lower entropic production bound is always larger than the other two types, we denote this sector of the phase space by S (red shaded region in Figure 3). Most regions, however, will display dominance of entropy production by one type of replicators for a range of densities, and shift dominance over another type for another range of x values (see Figure 3b,d-f).
The lines separating sections of LEB dominance are given by the following set of inequalities, all derived from the results above: with the associated functions Notice that, in several patches of the space of parameters depicted in Figure 3, LEB dominance is dependent on specific density values. Also, ∆LEB(r|r ) functions behave such that LEB dominance always appears ordered as P, S and H, respectively. This ordered sequence can be understood as an indication of an underlying thermodynamical constraint for these pre-biotic replicating systems. Finally, notice that this analysis has been performed with fixed value of α. Nonetheless, shifting the values of this internal parameter does not substantially modify the structure of the phase space given in Figure 3, in fact, its topological arrangement will remain invariant.
Hence, from macroscopical considerations involving both coarse-grained values for the coupling constants {g, h, c} and internal parameter α, we are able to derive a phase space compartmentalisation that allows a classification based on the lower (generalized) entropy production bounds for each replicator type. A qualitative tendency emerges from this picture: the parabolic replicator generates more entropy at low densities while so does the hyperbolic at high x values, leaving the simple replicator in between.

Discussion
A significant gap in our understanding of evolution, particularly in relation with early events and simple living systems, stems from the lack of a physical theory incorporating a thermodynamic description of replication dynamics. Self-replication stands as the one characteristic feature of living matter and its singular character was early appreciated by theoreticians when comparing cells and machines [44,45]. This work was an important step towards an understanding of the logic and computational nature of self-replicating agents. But a physical equivalent addressing the fundamental physics bounds to replication has been missing.
Recent work has addressed this problem revealing a powerful connection between entropy production and the transition probabilities underlying a stochastic, microscopic description [17,46]. Such connection can be efficiently exploited to analyse, under the coarse-graining described above, the general tendency of a Darwinian replicator to replicate itself. In this way, it is possible in particular to compare the efficiency of different classes of replicators by looking at their relative lower entropy bounds.
Instead of a direct comparison of the systems' measurable replication rates, this framework focuses on how, via a coarse-graining procedure, these parameters are resulting from the interplay of the many internal degrees of freedom. This technique ultimately leads to the estimation of the lower entropic bounds for each replicator. We interpret these non-equilibrium thermodynamic bounds as a consistent way of comparing and evaluating the likelihood of observing different classes of replicators. This is summarised in the phase diagram shown in Figure 3 where the relative dominance of each class is indicated. Notice that the analysis above does not involve competition between the replicator classes. All computations for the entropic bounds are done by considering the replicators to be evolving separately (see Appendix B for details).
Congruent approaches have been recently put forward following an equivalent theoretical formalism studying the non-equilibrium costs of production and destruction of polymers [47]. On the other hand, the present approach ought to be regarded as a minimal theoretical setup, and a number of issues can be raised. For instance, the fact that prebiotic systems might have exploited physical environments where sharp gradients are present, as it occurs with water-air interfaces [48]. Further developments in non-equilibrium statistical physics are needed in order to tackle these types of heterogeneities.
Even at this level of description, we can see how the coarse-graining predicts what to expect for the constraints operating on the classes of replicators in early evolutionary stages. The diagrams reveal the threshold conditions that would allow particular types of replicators to thrive or coexist in a competing scenario. In some domains only Malthusian dynamics are thermodynamically dominant, while, in others, parabolic replicators seem to be more efficient at generating entropy. Also, in some regions, a combination of parabolic and hyperbolic (cooperator) agents would share dominance. Overall, there is a robust characterisation of dominance related to the density of the system, revealing a preferential order as we move from low to high densities. Future work should be aimed at the construction of theoretical microscopic models such that coarse-graining operations can be unambiguously defined and subsequent operations may be computed in order to obtain the emergent transition rates. This would yield a deeper understanding of both the coarse-graining process and how some biological systems seem to be able to operate at the edge of what is possible. Such an approach can lead to novel insights into the problem of how major evolutionary transitions (which are often tied to the emergence of novel classes of replicators) occur.

Acknowledgments:
The authors thank Thomas Ouldridge for personal communication, and the Complex Systems Lab members for fruitful discussions. Special thanks to Gerda Taro for her inspiring ideas. This work was supported by the Botín Foundation by Banco Santander through its Santander Universities Global Division, the Spanish Ministry of Economy and Competitiveness, grant FIS2015-67616-P and the Santa Fe Institute.
Author Contributions: Jordi Piñero & Ricard Solé concieved, designed, perfomed the mathematical analysis, and wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A
The argument for the effective kinetic law for the parabolic replicator goes as following: let y = [AA] (concentration of associated molecules), and z = [A] (concentration of dissociated molecules). Thus, define x = 2y + z as the gross stoichiometric concentration of molecules of type A, regardless of configuration. Assuming that the time-scale of the replication reaction (here moduled by ratio c) is much larger than the association/dissociation processes, then, by focusing on the dynamics of replication, we can assume balanced equilibrium Then, analysing the dynamics of the replication reaction, which is modulated by the parameter c, while the kinetics for the gross concentration x is obtained by using (A1) as but, as a b, then the equilibrium of the association/dissociation reaction is very much unbalanced in favour of the associated molecular configuration AA, which implies that x ≈ 2y. Thus, we conclude that the kinetics for x is given by Truncating at leading terms in (b/a), we derive with ρ = c 2b a .

Appendix B
Consider a well-mixed urn filled with N elements that can be characterised as dead or alive. Notice that such a characterisation embodies some kind of coarse-grained measure, since we are deliberately ignoring (integrating) all internal degrees of freedom for each element. Denote by n < N the number of active (alive) elements in the urn at a given time t. In the following sections we will derive the coarse-grained (mesoscopic) time-dependent dynamics. Hence, for each replicator type, let us construct a master equation of the form while restricting the dynamics to a first-step process and introducing a natural (single-particle) decay process modulated by parameter δ that will be equivalent to all replicating motifs. It is important to remark that these systems are implicitly open. This is because every time an active element turns into an inactive one what really is happening it is flowing out of the system and letting new source (E) flow in. In this sense, the proposed models are analogous to chemostats.

Appendix B.1
Beginning with the simple replicator, introduce the following rules (see Figure A1):

1.
Pick an element of the urn at random.

2.
If active, with probability g, pick a second element at random and (if not active) activate.

3.
Pick an element at random again.

4.
If active, with probability δ, deactivate. b a c Figure A1. A summary of the rules of replication in an urn model. Active chains are drawn as filled balls and inactive chains are white balls. (a) represents the action of selecting an active chain replicating following the simple replicator mechanism; (b) shows the replicating process of a hypercyclic replicator; (c) corresponds to the decay which, for the purpose of this work, is supposed to act equivalently in each replicator-type.
For the simple replicator, using the notation on Table 1, plus adding a single-particle decay process, Then, for sufficiently large N, it is possible to approach the dynamics by dP(n, t) dt Consider now the dynamics of hyperbolic replicators. Following the rules summarized in Table 1, introduce the algorithm (see Figure A1):

1.
Pick an element of the urn at random.

2.
If active, pick a second element at random. 3.
If active, with probability h, pick a third element at random and (if not active) activate.

4.
Pick an element at random again. 5.
Hence, restricting the dynamics to a first-step process, we may deduce the following transition probabilities (cf. [49,50]) which, for N 1, lead to the master equation

Appendix B.3
Finally, let us derive the macroscopical dynamics for a parabolic replicator by implementing the following set of rules on an urn of N elements (see Figure A2): 1. Pick an element of the urn at random. If active, then: (i) if in associated state (AA) then, with probability a, dissociate and iterate; (ii) if dissociated, pick a second element and, if active, with probability b, associate. Iterate this process until equilibrium is reached for association/dissociation reaction. 2. Pick an element of the urn at random. If active, pick a second element at random, if empty, with probability c, replicate. 3. Pick an element of the urn at random. If active, with probability δ, deactivate. a b c Figure A2. This diagram shows how the urn model of parabolic replicators is implemented. (a,b) correspond to the rapid association/dissociation reactions, which are supposed to equilibrate in much shorter time-scales than the replicating process, which is shown in (c), i.e., τ 0 τ 1 . The process of equilibration (left box) is iterated a large number of times before the loop goes into the replication process (right box).
The situation for the parabolic replicator is a peculiar one, for one thing, it involves two characteristic time-scales, a rapid one, concerning the association/dissociation process (see Appendix A), and the replication process. In order to approximate the transition rates let us define k as the number of associated pairs, AA, and m as the number of dissociated active elements in the urn, A. Let N be the total number of elements in the urn, including associated, dissociated and deactivated elements. Let us denote by n the total number of active elements, regardless of configuration, then, n = 2k + m. Now, assuming rapid equilibration of the association/dissociation reaction in (14), which can be related to the number n by where we neglect the negative root, as it is non-physical. Hence, it is now possible to construct a master equation for the first-step process of replication as in (14), with which, for N 1, and using (A12) lead to dP(n, t) − δ n N [P(n, t) − P(n + 1, t)] ,