Thermodynamics of Darwinian evolution in molecular replicators

Abstract

We consider the relationship between thermodynamics, fitness, and Darwinian evolution in autocatalytic molecular replicators. We uncover a thermodynamic bound that relates fitness, replication rate, and the Gibbs free energy dissipated per copy. This bound applies to a broad range of systems, including elementary and non-elementary autocatalytic reactions, polymer-based replicators, and certain kinds of autocatalytic sets. In addition, we show that the critical selection coefficient (the minimal fitness difference visible to selection) is bounded by the Gibbs free energy dissipated per replication. Our results imply fundamental thermodynamic bounds on the strength of selection in molecular evolution, complementary to other bounds that arise from finite population sizes and error thresholds. These bounds may be relevant for understanding thermodynamic constraints faced by early replicators at the origin of life. We illustrate our approach on several examples, including a classic model of replicators in a chemostat.

I Introduction

Recent work has uncovered fundamental bounds on the thermodynamic costs of various biomolecular processes, including chemical sensing (mehta2012energetic, ; barato_efficiency_2014, ; govern2014energy, ), copying of polymer-stored information (andrieuxNonequilibriumGenerationInformation2008, ; ouldridge2017fundamental, ; poultonNonequilibriumCorrelationsMinimal2019, ; sartoriKineticEnergeticDiscrimination2013, ), and growth and replication (kondo_growth_2011, ; england2013statistical, ; himeoka_entropy_2014, ; virgo2016complex, ; saakianNonlinearStochasticDynamics2016, ; bishopStochasticBistabilityBifurcation2010, ; pinero_nonequilibrium_2018, ; corominas2019thermodynamics, ). These results are derived from general principles of nonequilibrium thermodynamics — such as flux-force relations and fluctuation theorems (beard2007relationship, ; jarzynski_equalities_2011, ; seifert2012stochastic, ; kondepudiModernThermodynamicsHeat2015, ) — that relate the dynamical and thermodynamic properties of nonequilibrium processes. Due to their generality, these bounds shed light on universal thermodynamic constraints on lifelike systems, including modern and protobiological organisms, synthetic life, and even possible non-terrestrial lifeforms.

One of the most important properties of living systems is that they undergo Darwinian evolution. Generally speaking, Darwinian evolution refers to a process in which organisms with higher fitness are able to outcompete organisms with lower fitness. Darwinian evolution can be exhibited by chemical systems, such as individual replicating molecules or networks of molecules (eigenSelforganizationMatterEvolution1971, ; schusterDynamicsEvolutionaryOptimization1985, ; leeAutocatalyticNetworksTransition1997, ; feistelPhysicsSelfOrganization2011, ; vasasEvolutionGenes2012a, ; takeuchiEvolutionaryDynamicsRNAlike2012, ; vaidyaSpontaneousNetworkFormation2012, ; szilagyi2017ecology, ; peng2020ecological, ; ameta2021self, ; mizuuchiEvolutionaryTransitionSingle2022, ; ametaDarwinianPropertiesTheir2021, ). Synthesizing such systems has been a major focus of research on the origin of life, given that the emergence of Darwinian evolution is considered a crucial point in the transition from nonliving to living matter (bernstein1983darwinian, ; yarusGettingRNAWorld2011, ; vasasEvolutionGenes2012a, ; ameta2021self, ; jeancolasThresholdsOriginLife2020, ).

Here, we consider the relationship between thermodynamics and Darwinian evolution in minimal chemical systems. This relationship may be particularly relevant for understanding thermodynamic constraints on the origin of life (pascal2013towards, ; jeancolasThresholdsOriginLife2020, ; dangerConditionsMimickingNatural2020, ). We consider a nonequilibrium reaction volume which contains autocatalytic replicators. These replicators may copy themselves via elementary reactions, or via more complex multi-step replication mechanisms. We also consider certain types of collectively autocatalytic sets, where replication involves a uniform cycle of cross-catalytic reactions. Our setup includes many types of molecular replicators previously considered in the literature, including self-complementary and complementary templates, polymer-based replicators, and autocatalytic small molecules (patzke2007self, ; bissetteMechanismsAutocatalysis2013a, ; lincolnSelfsustainedReplicationRNA2009, ; leeSelfreplicatingPeptide1996, ). It also encompasses several classic models of molecular replicators, including the chemostat model (schusterDynamicsEvolutionaryOptimization1985, ; feistelPhysicsSelfOrganization2011, ) and Eigen’s quasispecies model (eigenSelforganizationMatterEvolution1971, ).

Each replicator is associated with three quantities. The first is the Gibbs free energy of replication, the free energy dissipated when one copy of the replicator is synthesized, which we write as $\sigma$ in dimensionless units of $k_{B}T$-per-copy. The second is the replication rate $\rho$, the number of copies per replicator and unit time under current conditions. Lastly, the fitness $f$ is the maximum possible replication rate, as achieved in the limit of small concentrations. In Section IV.1, we show that fitness determines a replicator’s ability to invade a given population (metzHowShouldWe1992, ) and to sustain a high dilution rate. The introduction of this operational definition of fitness for molecular replicators is an important contribution of our work.

In Section IV.2, we derive a thermodynamic bound that relates Gibbs free energy of replication, the replication rate, and the fitness as

$$\sigma\geq-\ln\left(1-\frac{\rho}{f}\right).$$ (1)

Thus, the Gibbs free energy of replication increases as the replication rate $\rho$ approaches that $f$, the maximum possible. Equivalently said, the Gibbs free energy of replication limits how close the replication rate can approach the maximum.

We also derive a thermodynamic bound on the strength of Darwinian evolution. Observe that a higher fitness replicator is not always able to outcompete a lower fitness one, as this depends on the difference in fitness as well as various environmental and demographic factors (gillespie2004population, ). The strength of Darwinian evolution can be quantified in terms of the smallest fitness difference that can affect evolutionary outcomes in a given population and environment. This is the so-called critical selection coefficient $s$, which sets a “resolution limit” for selection beyond which fitness differences are indiscernible. As an example, the strength of Darwinian evolution in finite populations is limited by stochasticity, and a fitter mutant will fixate with significantly higher probability only when $s\gg 1/N_{e}$, where $N_{e}$ is the effective population size (ewensMathematicalPopulationGenetics2004, ). As another example, Eigen’s work on the “error threshold” shows that the strength of Darwinian evolution is limited by the mutation rate $\mu$, so that a fitter mutant can dominate the population only when $s>\mu$ (Eq. II-45, eigenSelforganizationMatterEvolution1971, ; smithMajorTransitionsEvolution, ).

We use the inequality (1) to derive a thermodynamic bound on the critical selection coefficient. We suppose that Darwinian evolution is sufficiently strong so that a replicator with fitness $f$ is present in the steady state of a flow reactor, while another replicator type with fitness $f^{\prime}<f$ is driven to extinction. We show that $s=1-f^{\prime}/f$, the selection coefficient between the two replicators, obeys

$$s\geq e^{-\sigma^{*}},$$ (2)

where $\sigma^{*}$ is the Gibbs energy dissipated by the fitter replicator in steady state. This bound on the strength of Darwinian evolution applies even in the case of infinite populations and error-free replicators. It implies that selecting for a relative fitness difference of $s$ requires dissipating more than $-\ln s$ of free energy, a quantity which diverges for vanishing fitness differences, $s\to 0$. We illustrate this result using a classic model of replicators in a chemostat in Section VI.

In Section VII, we extend our results to autocatalytic sets where replication involves a cycle of cross-catalytic reactions. We first generalize our notion of fitness to autocatalytic sets and then derive generalized versions of the inequalities (1) and (2). In these generalizations, $\sigma$ refers to the Gibbs free energy dissipated by the average cross-catalytic reaction in the cycle.

Figure 1: Illustration of our thermodynamic bound (2) for three real-world molecular replicators: a self-replicating prion (baskakovFoldingPrionProtein2001, ), an RNA molecule that copies itself via a single ligation (lincolnSelfsustainedReplicationRNA2009, ), and a peptide that copies itself via “native chemical ligation” (leeSelfreplicatingPeptide1996, ; dawsonSynthesisProteinsNative1994, ).

In principle, our results can also be applied to various real-world molecular replicators, suggesting a route for experimental validation. In Fig. 1, we use published thermodynamic data to illustrate the bound (2) on three real world systems. The first is a self-replicating prion prusiner1991molecular , where $\sigma^{*}\approx 3.5$ when native and misfolded concentrations are equal baskakovFoldingPrionProtein2001 ; Note101 . The second is a self-replicating peptide that copies itself using “native chemical ligation”  leeSelfreplicatingPeptide1996 ; wangTheoreticalAnalysisDetailed2011 ; Note102 , where $\sigma^{*}\approx 5.9$ at substrates concentration $90\,\mu\mathrm{M}$ and replicator concentration $5\,\mu\mathrm{M}$ (leeSelfreplicatingPeptide1996, ). The third is a self-replicating RNA molecule that copies itself using a single RNA ligation (lincolnSelfsustainedReplicationRNA2009, ), with $\sigma^{*}\approx 5$ under in vivo concentrations (julicherMotionRNAPolymerase1998, ; erieSingleNucleotideAdditionCycle1992, ). For example, our result predicts that for the RNA replicator, selection can only discern relative fitness differences of $s\geq e^{-5}\approx 0.6\%$.

II Setup

We consider a reaction volume at constant temperature and pressure. The reaction volume contains an ideal well-mixed solution of replicators $X,X^{\prime},\dots$ and other chemical species $A_{1},A_{2},\dots$ that may serve as substrates and side products. We study this system in terms of deterministic concentrations, assuming that molecular counts are sufficiently large so that stochastic fluctuations can be ignored. We use $x$ to indicate the concentration of replicator $X$ and $\vec{a}=(a_{1},a_{2},\dots)$ to indicate the concentrations of substrates/side products $(A_{1},A_{2},\dots)$.

Each replicator $X$ undergoes a reversible autocatalytic reaction of the form

$$X+\sum_{i=1}\alpha_{i}A_{i}\rightleftharpoons X+X+\sum_{i=1}\beta_{i}A_{i},$$ (3)

where $\alpha_{i}$ and $\beta_{i}$ are stoichiometric coefficients of substrate/side products $A_{i}$. A simple example of Eq. 3 is autocatalysis from a single substrate, $X+A\rightleftharpoons X+X$, but many other schemes are also possible. The autocatalytic reaction may be elementary, or it may proceed via multiple steps. Different types of replicators will generally have different stoichiometric coefficients $\alpha_{i},\beta_{i}$ as well as other thermodynamic and kinetic parameters (discussed below).

We use $\mathcal{J}$ to indicate the net flux across the autocatalytic reaction in Eq. 3. We define the replication rate $\rho$ of replicator $X$ as the net flux per replicator,

$$\rho:=\frac{\mathcal{J}}{x}.$$ (4)

The replicators may also flow out of the volume with dilution rate $\phi\geq 0$. Accounting for both replication and dilution, the concentration of replicator $X$ changes as

$$\dot{x}=\mathcal{J}-\phi x=(\rho-\phi)x.$$ (5)

We usually leave dependence on time $t$ implicit in our notation. We will use that in steady state, any non-extinct replicator ($x>0$) must have $\rho=\phi$, meaning that dilution and replication balance.

In writing Eq. 5, we assume that the rate of non-autocatalytic formation of the replicator is negligible. We also assume that different replicators do not interact directly by consuming each other as substrates or producing each other as side products, although they may interact indirectly via shared substrates/side products $A_{i}$. For simplicity, we also ignore the spontaneous degradation of replicators; however, in Appendix A, we show that our results still hold in the presence of degradation reactions.

Some of our results are valid for closed reactors as well as open reactors that exchange matter with their environment (avanzini2022thermodynamics, ). However, our thermodynamic bound on selection (2) applies specifically to a flow reactor in steady state. Different types of flow reactors may be considered. One example is the continuous stirred-tank reactor (CSTR) where the dilution rate and inflow rates are constant, which is often used in chemical (filisettiStochasticModelEmergence2011a, ; semenovAutocatalyticBistableOscillatory2016, ) and biological experiments (dykhuizen1983selection, ; groverResourceCompetition1997, ; hoskisson2005continuous, ), and which can also arise naturally (e.g., in a pond fed by a nutrient-rich stream). To avoid confusion, we note that the CSTR is called a chemostat in the biological literature (harmandChemostat2017, ; smithTheoryChemostatDynamics1995, ), although in nonequilibrium thermodynamics, the term chemostat is usually used instead to refer to an external chemical reservoir (polettini2014irreversible, ). Another possible type of flow reactor is one where the rates of dilution and inflow can vary as a function of the chemical concentrations. An example is provided by Eigen’s quasispecies model, where the dilution rate is adjusted to maintain the total concentration of replicators constant (eigenSelforganizationMatterEvolution1971, ).

The Gibbs free energy of the replication reaction (3) is (kondepudiModernThermodynamicsHeat2015, ; raoNonequilibriumThermodynamicsChemical2016, )

$$\sigma=-\ln x+\sum_{i=1}(\alpha_{i}-\beta_{i})\ln a_{i}-\Delta G^{\circ}/RT,$$ (6)

where $R$ is the gas constant, $T$ is the temperature, and $-\Delta G^{\circ}$ is the standard Gibbs free energy,

$$-\Delta G^{\circ}=RT\Big{[}-\ln x^{\text{eq}}+\sum_{i=1}(\alpha_{i}-\beta_{i})\ln a_{i}^{\text{eq}}\Big{]}.$$ (7)

Here $x^{\mathrm{eq}}$ and $a_{i}^{\mathrm{eq}}$ are equilibrium concentrations of $X$ and $A_{i}$, as would be reached if the reactor was closed to exchange of matter and allowed to relax completely. We note that we express $\sigma$ in dimensionless units of entropy per copy. It can be converted to joules per copy as $k_{B}T\sigma$, where $k_{B}$ is Boltzmann’s constant, or to joules per mole as $RT\sigma$. In the chemistry literature, the Gibbs free energy of reaction is often written as $-\Delta G$. The Gibbs free energy of replication can be understood as a fundamental thermodynamic cost that represents lost work potential: a replication reaction that dissipates $\sigma$ of free energy can be coupled to a thermodynamically disfavored “uphill” reaction and thereby perform up to $\sigma$ of chemical work.

III Elementary and non-elementary replicators

III.1 Elementary replicators

An elementary replicator refers to the case where Eq. 3 is an elementary reaction. Then, the net flux across the reaction has the mass-action form (kondepudiModernThermodynamicsHeat2015, )

$$\mathcal{J}=rx-r^{-}x^{2},$$ (8)

where the forward and backward (pseudo) rate constants are

$$r=\kappa\prod_{i}a_{i}^{\alpha_{i}},\qquad r^{-}=\kappa e^{\Delta G^{\circ}/RT}\prod_{i}a_{i}^{\beta_{i}},$$ (9)

and $\kappa$ is a baseline rate constant. Note that $r$ and $r^{-}$ depend on the concentrations $\vec{a}$, although we always leave this dependence implicit in our notation.

The log-ratio of the forward and backward fluxes in Eq. 8 equals the Gibbs free energy of replication,

$$\sigma=\ln\frac{rx}{r^{-}x^{2}}=\ln\frac{r}{r^{-}x}\,,$$ (10)

as can be verified by comparing Eqs. 9 and 6. Eq. 10 is an instance of a general result, called the flux-force relation or local detailed balance in the literature (wachtelThermodynamicallyConsistentCoarse2018, ; raoNonequilibriumThermodynamicsChemical2016, ), which relates the Gibbs free energy of an elementary reaction with the forward and reverse fluxes (kondepudiModernThermodynamicsHeat2015, ). The flux-force relation is one of the most important results in nonequilibrium thermodynamics, since it connects the kinetic properties of a chemical reaction with its thermodynamic properties.

Actually, our results will not require that the rate constants $r$ and $r^{-}$ have the mass action form of Eq. 9, only that Eq. 8 and Eq. 10 hold. This could be used to study certain non-ideal solutions that exhibit interactions between $A_{i}$ (avanzini2021nonequilibrium, ).

III.2 Non-elementary replicators

Most real-world autocatalytic replicators cannot be treated as elementary reactions. For this reason, we also consider the case where Eq. 3 represents a reaction mechanism that proceeds via a sequence of $m$ steps,

$$\setchemfig{arrowcoeff=0.6}\schemestart X\arrow{<U=>[{\footnotesize\sum_{i}\alpha_{1,i}A_{i}}][{\footnotesize\sum_{i}\beta_{1,i}A_{i}}]}Y_{1}\arrow{<U=>[{\footnotesize\sum_{i}\alpha_{2,i}A_{i}}][{\footnotesize\sum_{i}\beta_{2,i}A_{i}}]}\dots\arrow{<U=>[\cdots][\quad\quad\cdots]}Y_{m-1}\arrow{<U=>[{\footnotesize\sum_{i}\alpha_{m,i}A_{i}}][{\footnotesize\sum_{i}\beta_{m,i}A_{i}}]}X+X\schemestop.$$ (11)

Each $Y_{k}$ is an intermediate chemical species, and each intermediate step is an elementary reversible reaction with mass-action kinetics that may involve substrates/side products $A_{i}$ with stoichiometric coefficients $\alpha_{k,i}$ and $\beta_{k,i}$. The stoichiometry of the overall reaction obeys $\alpha_{i}=\sum_{k}\alpha_{k,i}$ and $\beta_{i}=\sum_{k}\beta_{k,i}$. We assume that intermediate species are not shared between different types of replicators. For simplicity, we also assumed that degradation reactions are negligible, although our results generalize to the presence of degradation as shown Appendix A.

Figure 2: Examples of non-elementary autocatalytic replication mechanisms. Left: autocatalysis with binding, conversion, and unbinding steps. Right: templated replication of a self-complementary polymer (shown here using a dimer).

We term this kind of reaction mechanism a non-elementary replicator, although it is also called an “autocatalytic cycle” in the literature (kingAutocatalysis1978, ; virgo2016complex, ; hordijkAutocatalyticConfusionClarified2017, ). A simple example of a non-elementary replicator is a three-step mechanism with binding, conversion, and unbinding steps, see Fig. 2 (left). Another example is the step-by-step replication of a self-complementary dimer, illustrated in Fig. 2 (right), studied in numerous origin-of-life experiments (von1986self, ; patzke2007self, ; bissetteMechanismsAutocatalysis2013a, ). Yet other examples include the formose cycle (bissetteMechanismsAutocatalysis2013a, ) and the autophosphorylation of protein kinase (wangAutophosphorylationKineticsProtein2002, ; bishopStochasticBistabilityBifurcation2010, ).

Since each intermediate reaction $k\in\{1,\dots,m\}$ is elementary with mass-action kinetics, the net flux across it is

$$J_{k} =\nu_{k}y_{k-1}-\nu_{k}^{-}y_{k}.$$ (12)

Here $y_{k}$ is the concentration of intermediate species $Y_{k}$ for $k\in\{1,\dots,m-1\}$, and we use the convention $y_{0}=x,y_{m}=x^{2}$. The terms $\nu_{k}$ and $\nu_{k}^{-}$ refer to forward and backward (pseudo) rates constants of intermediate steps, typically defined in a manner analogous to Eq. 9. We note that $\nu_{k}$ and $\nu_{k}^{-}$ can depend on the concentrations $\vec{a}$, although we leave this implicit in our notation.

The production of replicator $X$ due to the autocatalytic mechanism is

$$\mathcal{J}=2J_{m}-J_{1},$$ (13)

since two copies are produced by the last step and one consumed by the first. The production of intermediate species $Y_{k}$ ($k\in\{1,\dots,m-1\}$) due to the mechanism is $J_{k}-J_{k+1}$, since one $Y_{k}$ is produced by intermediate step $k$ and one is consumed by intermediate step $k+1$. The rate of change of the concentration of the intermediate species, also accounting for dilution with rate $\phi\geq 0$, is

$$\dot{y}_{k}=J_{k}-J_{k+1}-\phi y_{k}.$$ (14)

It will be useful to express the overall flux $\mathcal{J}$ in a mass-action-like form of Eq. 8. To do so, we introduce the assumption that the relative concentrations of intermediate species $Y_{k}$ and replicator $X$ is approximately constant,

$$\frac{d}{dt}\frac{y_{k}}{x}=\frac{x\dot{y}_{k}-\dot{x}y_{k}}{x^{2}}\approx 0.$$ (15)

This assumption always holds in steady state ($\dot{y}_{k}=\dot{x}=0$). It also holds under a separation of timescales, when the relative concentrations $y_{k}/x$ quickly reach stationary values — that is, faster than those stationary values are perturbed by changes in other variables like the dilution rate and the absolute concentrations of $x$ and $\vec{a}$. This separation of timescales holds during the exponential growth phase of an initially rare replicator, when the absolute concentrations of $x$ and $y_{k}$ are very small and have a minimal effect on other variables.

We emphasize that the assumption (15) is different from the quasi-steady-state (QSS) approximation often used in biochemistry (segel1989quasi, ). QSS assumes that the absolute concentrations of intermediate species are approximately stationary, compared to the rate of change of the replicator concentration $x$. The QSS assumption can be strongly violated in growing autocatalytic replicators, especially when the intermediate species have large relative concentrations.

Plugging $\dot{x}=(\rho-\phi)x$ and $\dot{y}_{k}=J_{k}-J_{k+1}-\phi y_{k}$ into Eq. 15 and simplifying gives

$$J_{k}-J_{k+1}=\rho y_{k}.$$ (16)

Using Eq. 12, we can express the intermediate concentrations $\vec{y}=(y_{1},\dots,y_{m-1})$ in terms of a linear system

$$[(M+\rho I)\vec{y}]_{k}=\delta_{k,1}\nu_{1}x+\delta_{k,m-1}\nu_{m}^{-}x^{2},$$ (17)

where $M\in\mathbb{R}^{(m-1)\times(m-1)}$ is a matrix of intermediate rate constants,

$$M\!=\!\left[\begin{array}[]{ccccc}\nu_{1}^{-}\!\!+\!\nu_{2}&-\nu_{2}^{-}&0&0&\dots\\ -\nu_{2}&\nu_{2}^{-}\!\!+\!\nu_{3}&-\nu_{3}^{-}&0&\dots\\ 0&\dots&\dots&\dots&\nu_{m-1}^{-}\\ \dots&0&0&-\nu_{m-1}&\nu_{m-1}^{-}\!\!+\!\nu_{m}\end{array}\right]$$ (18)

For any $\rho\geq 0$, $M+\rho I$ is an “M-matrix” (shivakumar1996two, ), so it is invertible and all of entries of its inverse are nonnegative (shivakumar1996two, ; meyer1978singular, ).

We solve Eq. 17 for $\vec{y}$ and combine with Eqs. 13 and 12 to write $\mathcal{J}$ in a mass-action-like form,

$$\mathcal{J}=r(\rho)x-r^{-}(\rho)x^{2},$$ (19)

with effective rate constants given by

$$\begin{aligned} r(\rho)&:=\nu_{1}(\nu_{1}^{-}G_{11}+2\nu_{m}G_{m-1,1}-1)\\ r^{-}(\rho)&:=\nu_{m}^{-}(2-2\nu_{m}G_{m-1,m-1}-\nu_{1}^{-}G_{1,m-1})\,,\end{aligned}$$ (20)

where $G:=(M+\rho I)^{-1}$. These effective rate constants can depend on concentrations $\vec{a}$, since $\nu_{k}$ and $\nu_{k}^{-}$ depend on them, although we leave this dependence implicit in our notation.

Below we will use that the first and second derivatives of the effective rate constants obey

$$\partial_{\rho}r(\rho) \leq 0,\;\partial_{\rho}^{2}r(\rho)\geq 0$$ (21)

$$\partial_{\rho}r^{-}(\rho) \geq 0,\;\partial_{\rho}^{2}r^{-}(\rho)\leq 0$$

This follows from Eq. 20 by using matrix calculus and the fact that all entries of $G=(M+\rho I)^{-1}$ are nonnegative.

In addition, for $\rho=0$, the effective rate constants can be expressed in closed form,

$$r(0) =\Bigg{[}\sum_{k=1}^{m}\frac{\prod_{l=1}^{k-1}\nu_{l}^{-}}{\prod_{l=1}^{k}\nu_{l}}\Bigg{]}^{-1}\;\quad r^{-}(0)=r(0)\prod_{k=1}^{m}\frac{\nu_{k}^{-}}{\nu_{k}}$$ (22)

as shown in Appendix B. The $\rho=0$ regime corresponds to the limit where the replication rate $\rho$ is much slower than the rate of internal reactions, so that $\rho$ can be neglected when solving the linear system (17).

We finish by discussing the thermodynamics of non-elementary replicators. The Gibbs free energy of intermediate reaction $k$ in Eq. 11 is

$$\sigma_{k}=\ln\frac{y_{k-1}}{y_{k}}+\sum_{i=1}(\alpha_{k,i}-\beta_{k,i})\ln a_{i}-\Delta G_{k}^{\circ}/RT,$$ (23)

where as above we use $y_{0}=x,y_{m}=x^{2}$. Since intermediate reactions are elementary, the forward and backward fluxes in Eq. 12 obey the flux-force relation,

$$\sigma_{k}=\ln\frac{\nu_{k}y_{k-1}}{\nu_{k}^{-}y_{k}},$$ (24)

as can be shown explicitly when the rate constants $\nu_{k}$ and $\nu_{k}^{-}$ have mass-action kinetics similar to Eq. 9.

The Gibbs free energy of the overall replication mechanism is the sum of the Gibbs free energy of the individual steps (kondepudiModernThermodynamicsHeat2015, ),

$$\sigma=\sum_{k=1}^{m}\sigma_{k}=\ln\frac{1}{x}+\sum_{k=1}^{m}\ln\frac{\nu_{k}}{\nu_{k}^{-}}.$$ (25)

We can express $\sigma$ in a “flux-force-like” expression by combining Eqs. 25 and 22,

$$\sigma=\ln\frac{r(0)}{xr^{-}(0)}=\frac{xr(0)}{x^{2}r^{-}(0)}.$$ (26)

Comparing to Eq. 19, we see that the numerator and denominator can be interpreted as the forward and backward fluxes in the case of $\rho=0$, that is the limit when dilution is much slower than internal reactions. Ref. (wachtelThermodynamicallyConsistentCoarse2018, ) previously derived this flux-force relation for steady state conditions and reaction mechanisms such as Eq. 11, which involve a single reaction cycle.

An elementary replicator is a special case of a non-elementary replicator with a single reaction ($m=1$) and no intermediate species. In this case, the effective rate constants in Eq. 20 lose dependence on $\rho$ and reduce to Eq. 22 which in turn reduces to Eq. 9. The mass-action form of Eq. 19 reduces to Eq. 8, and the flux-force relation (26) is equivalent to Eq. 10.

IV Fitness, selection, and thermodynamics

IV.1 Fitness and selection coefficient

The concept of fitness can be defined differently in different evolutionary scenarios, and finding an appropriate definition is an important area of research in biology and ecology (metzHowShouldWe1992, ; roffDefiningFitnessEvolutionary2008, ; orrFitnessItsRole2009, ; spaakDifferentMeasuresNiche2023, ). Here, we propose a definition of fitness that is suitable for autocatalytic molecular replicators, both elementary and non-elementary.

Before proceeding, we note that one could simply define fitness as the replication rate $\rho$ at a given point in time. However, this approach runs into problems. For instance, for a reactor in steady state, all present replicators have the same replication rate, i.e., the steady-state dilution rate $\phi$, while all non-present replicators have an undefined replication rate. This makes it impossible to ask many natural questions, such as whether higher fitness replicators outcompete lower fitness ones. More generally, replication rate $\rho$ is a statistic of actual performance that does not specify how a replicator would perform in a new environment — which is usually what is desired from a fitness measure (metzHowShouldWe1992, ; roffDefiningFitnessEvolutionary2008, ; orrFitnessItsRole2009, ; spaakDifferentMeasuresNiche2023, ).

Instead, we define the fitness $f$ of replicator $X$ via the equation,

$$f=r(f)\geq 0,$$ (27)

where $r(\cdot)$ is the forward rate rate constant defined in Eq. 20. For an elementary replicator, $r(\rho)$ does not depend on $\rho$ and $f$ is simply the rate constant of the elementary reaction in Eq. 9. For a non-elementary replicator, $f$ is the nonnegative root of the algebraic expression $\alpha-r(\alpha)$. This expression is strictly increasing in $\alpha$, ranging from $0-r(0)\leq 0$ to $r(0)-r(r(0))\geq 0$, as can be deduced from Eqs. 22 and 21. There is a unique value of $0\leq f\leq r(0)$ that satisfies Eq. 27 and that can be found quickly using numerical methods such as bisection.

We note that fitness depends on the concentrations of substrates and side products $\vec{a}$, although this is left implicit in our notation. These concentrations specify the “ecological environment” of the reaction volume.

Our definition of fitness is experimentally measurable and operationally meaningful. In particular, $f$ can be understood as the initial growth rate that is achieved at small concentrations. Imagine a reactor in steady state that does not contain replicator $X$ at time $t<0$, and suppose that $X$ is introduced at a small concentration $x(0)$ at $t=0$. Suppose also that $X$ is an elementary replicator or a non-elementary replicator that obeys Eq. 15, so that relative concentrations of intermediates are approximately stationary. Assuming for simplicity that there is no dilution or degradation, the replicator’s concentration will initially grow as

$$\dot{x}=\mathcal{J}=\rho x\approx r(\rho)x,$$

as follows from Eqs. 4 and 19 and dropping second order terms in $x$. This equation $\rho=r(\rho)$ is uniquely satisfied by $\rho=f$, meaning that the concentration will grow as

$$x(t)\approx e^{tf}x(0)$$ (28)

for small $t$. Thus, $f$ has a clear operational interpretation: an initially rare mutant with fitness $f$ will increase in concentration if $f>0$ and decrease toward extinction if $f<0$. The same derivation holds for a flow reactor with dilution rate $\phi$, in which case it gives $x(t)\approx e^{t(f-\phi)}x(0)$. In biology, this type of fitness measure is called invasion fitness (liModelingMicrobialMetabolic2020a, ; arnoldiInvasionsEcologicalCommunities2022a, ; metzHowShouldWe1992, ; spaakDifferentMeasuresNiche2023, ). It has been argued to be a particularly general way to operationalize fitness (metzHowShouldWe1992, ; roffDefiningFitnessEvolutionary2008, ).

The initial growth rate is experimentally accessible as long as the initial concentration $x(0)$ is sufficiently small so that one can neglect the reverse flux (second term in Eq. 19) as well as perturbations to steady-state values of $\vec{a}$ and $\phi$ over the measurement timescale, but sufficiently large so that stochastic fluctuations can be ignored.

There is also another interpretation of fitness as the critical dilution rate, a quantity that plays an important role in chemostat studies (pirtExtensionTheoryChemostat1970, ; huEngineeringPrinciplesBiotechnology2017, ). Imagine a steady-state flow reactor with substrate concentrations $\vec{a}$ and steady-state dilution rate $\phi>0$. Suppose that the reactor contains replicator $X$ at steady-state concentration $x^{*}$. Using Eq. 19, we can express $x^{*}$ as

$$x^{*}=\begin{cases}\frac{r(\phi)-\phi}{r^{-}(\phi)}&\text{if\,}r(\phi)>\phi\\ 0&\text{otherwise}\end{cases},$$ (29)

where we used that $\rho=\phi$ in steady state. Now imagine slowly increasing the dilution rate $\phi$ while maintaining the concentrations $\vec{a}$ constant. The replicator will be pushed to extinction at the critical dilution rate $\hat{\phi}$ where $\hat{\phi}=r(\hat{\phi})$. Given Eq. 27, the critical dilution must equal the fitness $f$. The critical dilution rate is experimentally accessible, as long as the dilution rate can be manipulated while maintaining the concentrations $\vec{a}$ at constant values. Interestingly, the critical dilution rate can be defined for hypercycles and other higher-order replicators for which the replication rate vanishes at small concentrations (schloglChemicalReactionModels1972a, ; szathmarySimpleGrowthLaws1991, ).

Most simply, our definition of fitness can be considered as maximum replication rate that can be achieved by the replicator. That is, the above results can be used to show that $\rho\leq f$ always. Thus, $f$ sets an upper bound on the replication rate which is approached in the limit of low concentrations.

In addition to fitness, we also make use of the notion of the selection coefficient from evolutionary biology. The selection coefficient is a normalized measure of relative fitness difference that ranges from 0 (no difference) to 1 (maximum difference). Given two replicators $X$ and $X^{\prime}$ with fitness values $f\geq f^{\prime}$, we use a common definition of the selection coefficient as (lenski1991long, )

$$s:=1-\frac{f^{\prime}}{f}.$$ (30)

In real-world experiments, it is often difficult to measure individual concentrations of replicator and intermediate species. Often what is measured is a weighted sum of concentrations,

$$\omega=c_{X}x+\sum_{i=1}^{m-1}c_{i}y_{i},$$ (31)

given some nonnegative coefficients $c_{X}$ and $c_{i}$. The interpretations of replication rate, fitness, and critical dilution rate generalize to this situation. That is, if the replicator and all intermediate species grow at the same rate $\rho$, then their weighted sum $\omega(t)$ will also grow at the same rate. Similarly, if the replicator and all intermediate species vanish at some critical dilution rate $\hat{\phi}$, then their weighted sum will also vanish at that dilution rate.

IV.2 Thermodynamic bounds

We now derive our thermodynamic bounds. We consider a reactor containing replicator $X$ with a non-negative growth rate $\rho\geq 0$. Our first bound relates the replicator’s Gibbs free energy of replication $\sigma$, fitness $f$, and replication rate $\rho$ as

$$\sigma\geq-\ln\left(1-\frac{\rho}{f}\right),$$ (32)

which appeared as inequality (1) in the Introduction. For the special case of a flow reactor in steady state, this bound can be written in terms of the steady-state dilution rate $\phi=\rho$ as

$$\sigma^{*}\geq-\ln\left(1-\frac{\phi}{f}\right),$$ (33)

where $\sigma^{*}$ is Gibbs free energy of replication under steady state concentrations.

To derive the bound above, we combine Eqs. 4 and 19 to express the replicator’s concentration as $x=(r(\rho)-\rho)/r^{-}(\rho)$. Plugging this into the expression of $\sigma$ in Eq. 26 gives

$$\sigma=\ln\left(\frac{r(0)}{r^{-}(0)}\frac{r^{-}(\rho)}{r(\rho)-\rho}\right).$$ (34)

Inequality (32) is then equivalent to $\frac{r(0)}{r^{-}(0)}\frac{r^{-}(\rho)}{r(\rho)-\rho}\geq f/(f-\rho)$ for $0\leq\rho\leq f$, which in turn is equivalent to the nonnegativity of the function

$$h(\rho):=\frac{r(0)}{r^{-}(0)}r^{-}(\rho)(f-\rho)-f(r(\rho)-\rho)\,.$$ (35)

Taking second derivatives and using Eq. 21 shows that $\partial_{\rho}^{2}h(\rho)\leq 0$, so $h$ is concave. Inspection shows that $h(\rho)=0$ for $\rho=0$ and $\rho=f$. Therefore, $h(\rho)\geq 0$ over $0\leq\rho\leq f$, proving the inequality (32).

Our result is tight for elementary replicators where the rate constants $r(\rho)$ and $r^{-}(\rho)$ do not depend on $\rho$, so $h(\rho)=0$ for all $\rho$. More generally, our bound tends to be tighter for “effectively” elementary replicators, which have fast internal rates so that $r(\rho)$ and $r^{-}(\rho)$ depend weakly on $\rho$. For general non-elementary replicators, the bound tends to be tighter near $\rho=0$ and $\rho=f$, corresponding to equilibrium and irreversible regimes, respectively.

We also derive a thermodynamic bound on the selection coefficient, which has implications for the strength of Darwinian evolution. We consider a flow reactor in steady state with dilution rate $\phi$. Suppose that the reactor contains some replicator $X$ with fitness $f>\phi$, and that there is a second replicator $X^{\prime}$ with fitness $f^{\prime}\leq\phi$ that is pushed to extinction. Plugging into inequality (33) gives

$$\sigma^{*}\geq-\ln\left(1-\frac{f^{\prime}}{f}\right)=-\ln s,$$ (36)

which appeared as inequality (2) in the Introduction. This result is illustrated on a model of replicators in a chemostat in Section VI.

To build intuition regarding the bound (36), we may consider two extreme situations. The first is equilibrium steady state, where the replication rate $\rho$ and the Gibbs free energy of replication $\sigma$ vanish for all replicators. Here, all replicators are present in positive equilibrium concentrations which do not depend on fitness, reflecting the fact that Darwinian evolution is impossible in equilibrium (eigenSelforganizationMatterEvolution1971, ; bernstein1983darwinian, ). At the other extreme is the irreversible regime, where each replicator copies itself at its maximum possible rate $\rho=f$ and $\sigma$ diverges. Typically, steady states do not exist in this regime even in the presence of dilution, since any replicator with $f>\phi$ will grow without bound and any replicator with $f<\phi$ will decay to extinction. In the special case that the fittest replicator satisfies $f=\phi$, a non-zero steady state exists which can harbor a single type of replicator (kingAutocatalysis1978, ; lifsonCrucialStagesOrigin1997, ), since a replicator with any other fitness value will go extinct. To summarize, all replicators coexist in equilibrium ($\sigma=0$), while only the fitness replicator can exist in steady state under the irreversible regime ($\sigma\to\infty$). Intermediate values of $\sigma$ interpolate between these two extremes, permitting steady-state coexistence of some but not all replicators.

The inequalities (32) and (36) are remarkably general, being independent of most details of the chemical system. For instance, they do not depend on the number of coexisting replicators, the substrates/side products involved in replication, whether the replicators copy themselves via elementary or non-elementary reactions, whether the steady state is near or far from equilibrium, etc. They also do not depend on the dynamical mechanism that leads to a particular steady state. For example, they do not depend on whether replicators experience competitive interactions (e.g., different replicators rely on the same substrate) or not (e.g., different replicators do not share substrates but differ in their kinetic parameters).

V Example: Self-complementary dimer

Figure 3: Fitness and thermodynamic bounds illustrated on Rebek’s self-complementary dimer (tjivikuaSelfreplicatingSystem1990a, ) with 4 reactions, as in Fig. 2 (right). (a) Fitness recovers the initial replication rate given a small starting concentration. We show time-dependent concentrations of replicator $x$, bound dimer $y_{XX}$, and weighted sum $\omega$ from Eq. 37. (b) Fitness recovers the critical dilution rate in a flow reactor. (c) The bound (32) relates replication rate $\rho$, fitness $f$, and Gibbs free energy of replication $\sigma$. It is shown on the same system as in (a). (d) The bound (33) relates fitness, steady-state Gibbs free energy $\sigma^{*}$, and dilution rate $\phi$ in a steady-state flow reactor. It is shown on the same system as in (b).

We now illustrate our results using a concrete example. We consider a non-elementary replicator that copies itself via the mechanism arrow coeff=0.6 $X+A$ <=> $XA$ <V=>[$B$][] $XAB$ <=> $XX$ <=> $X+X$

which is also shown in Fig. 2 (right). $X$ is self-complementary dimer while $A$ and $B$ are substrates. The reaction $XAB\rightleftharpoons XX$ is a ligation that produces the bound dimer $XX$. This type of system was studied in many early experiments on self-replicating chemical systems (von1986self, ; zielinskiAutocatalyticSynthesisTetranucleotide1987a, ; tjivikuaSelfreplicatingSystem1990a, ; rotelloSigmoidalGrowthSelfreplicating1991, ; von_kiedrowski_parabolic_1991, ).

We parameterize the rate constants of the two binding reactions, $X+A\rightleftharpoons XA$ and $XA+B\rightleftharpoons XAB$, as

$$\nu_{1}=\kappa a \qquad\nu_{1}^{-}=\kappa e^{\Delta G_{1}^{\circ}/RT}$$

$$\nu_{2}=\kappa b \qquad\nu_{2}^{-}=\kappa e^{\Delta G_{2}^{\circ}/RT}$$

where $a$ and $b$ are concentrations of $A$ and $B$. The rate constants of the ligation step are parameterized as

$$\nu_{3}=1\qquad \nu_{3}^{-}=e^{\Delta G_{3}^{\circ}/RT}$$

and of the dimerization step $XX\rightleftharpoons X+X$ as

$$\nu_{4}=\kappa e^{-\Delta G_{4}^{\circ}/RT}\qquad\nu_{4}^{-}=\kappa$$

Where possible, we use parameter values from Rebek’s system, one of the first molecules that exhibited self-replication in the lab (tjivikuaSelfreplicatingSystem1990a, ; vonkiedrowskiMinimalReplicatorTheory1993a, ). The binding and unbinding reactions are assumed to be essentially in equilibrium, so we use a fast rate constant $\kappa=10^{9}$. The standard Gibbs free energies for the binding reactions are $-\Delta G_{1}^{\circ}/RT=-\Delta G_{2}^{\circ}/RT=\ln 60$ (favoring binding), and for dimerization it is $-\Delta G_{4}^{\circ}/RT=-\ln 630$ (favoring the bound dimer) (tjivikuaSelfreplicatingSystem1990a, ). The ligation step is assumed to be effectively irreversible, so we use a large standard Gibbs free energy of $-\Delta G_{3}^{\circ}/RT=10$.

Following experimental setups (vonkiedrowskiMinimalReplicatorTheory1993a, ), we consider the weighted sum of concentrations

$$\omega=x+y_{XA}+y_{XAB}+2y_{XX},$$ (37)

which is the total concentration of replicator and intermediates, with the bound dimer counting as two copies.

We first consider a closed reactor ($\phi=0$) which starts from nonequilibrium initial concentrations of substrates $a(0)=b(0)=8.2\,\mathrm{mM}$ (tjivikuaSelfreplicatingSystem1990a, ). We choose $x(0)=.1\,\mu\mathrm{M}$ for the initial replicator concentration. For these parameter values and substrate concentrations, we can use Eq. 20 and Eq. 27 to compute the fitness as

$$f\approx 0.14$$

Fig. 3 (a) shows the time-dependent concentrations, demonstrating that fitness (dashed line) captures the initial growth rate.

We then consider the same system, but now in a flow reactor in steady state. In Fig. 3 (b), we show steady-state concentrations across different dilution rates, with substrate concentrations maintained at steady-state values $a=b=8.2\,\mathrm{mM}$. As shown, fitness captures the critical dilution rate at which the replicator goes extinct.

Next, we illustrate our thermodynamic bounds on the same system. Fig. 3 (c) shows the Gibbs free energy of replication versus the bound (32) for the system considered in Fig. 3(a), where the replicator grows from a small initial concentration in a closed reactor. The bound is tightest in the regime of low concentrations, achieving an efficiency of $-\ln(1-\rho/f)/\sigma\approx 0.5$ around $x\approx.1\,\mu\mathrm{M}$. In Fig. 3 (d), we show the Gibbs free energy of replication versus our bound (33) for the steady-state flow reactor. The bound is tightest near the critical dilution rate, where it achieves an efficiency of $-\ln(1-\phi/f)/\sigma^{*}\approx 0.65$.

It should be noted that Rebek’s system, like many other self-complementary dimers, is sometimes described as a “parabolic replicator” with growth obeying a square-root law, $\dot{\omega}\propto\omega^{1/2}$. Square-root growth results because the autocatalytic mechanism becomes concentrated in the thermodynamically-favored dimerized form $XX$. In general, however, the type of growth varies with replicator concentration (vonkiedrowskiMinimalReplicatorTheory1993a, ). At low concentrations, the concentration of the dimerized form is small and growth is first order, as seen in the semi-logarithmic plot Fig. 3 (a). Square-root growth only appears at larger concentrations, once the concentration of the dimerized form (orange line in the figure) is sufficiently large.

VI Example: Darwinian evolution in a chemostat

We now illustrate our thermodynamic bound on selection using a classic model of autocatalytic replicators in a chemostat (continuous stirred-tank flow reactor) (schusterDynamicsEvolutionaryOptimization1985, ).

The reactor may contain up to $N$ replicator types, indicated as $X_{1},\dots,X_{N}$, where each $X_{i}$ undergoes an autocatalytic reaction

$$X_{i}+A\rightleftharpoons X_{i}+X_{i}$$ (38)

from a shared substrate $A$. The substrate $A$ is supplied at a constant rate $\phi$ and feed concentration $\gamma$, and all species flow out with constant dilution rate $\phi$.

For simplicity, we suppose that all replicators copy themselves via elementary autocatalytic reactions. The dynamics of concentrations of replicators $x_{i}$ and substrate $a$ are

$$\begin{aligned} \dot{x}_{i}&=k_{i}x_{i}(a-e^{\Delta G_{i}^{\circ}}x_{i})-\phi x_{i}\\ \dot{a}&=\phi(\gamma-a)-\sum_{i}k_{i}x_{i}(a-e^{\Delta G_{i}^{\circ}}x_{i}),\end{aligned}$$ (39)

where $k_{i}$ is a rate constant and $-\Delta G_{i}^{\circ}$ is the standard Gibbs free energy of the reaction (38). As usual, we leave dependence on time of $x(t)$ and $a(t)$ implicit.

Although the replicators do not interact directly, they experience an effective interaction due to competition for the shared substrate $A$. This system is related to models of resource competition studied in mathematical ecology and evolutionary biology (dykhuizen1983selection, ; smithTheoryChemostatDynamics1995, ; groverResourceCompetition1997, ; harmandChemostat2017, ). Moreover, this system can be mapped onto a competitive Lotka-Volterra system with an effective interaction (see Appendix C).

This type of dynamical system was considered by Schuster and Sigmund (schusterDynamicsEvolutionaryOptimization1985, ) (see also (feistelPhysicsSelfOrganization2011, )). They showed that for any strictly positive initial conditions, there is a unique steady state which governs the long-term behavior. This steady state is given by a set of coupled equations,

$$a^{*}=\gamma-\sum_{i}x_{i}^{*},\quad x_{i}^{*}=\max\{0,e^{-\Delta G_{i}^{\circ}}(a^{*}-\phi/k_{i})\}.$$ (40)

In Appendix C, we show how to solve the coupled equations in Eq. 40 by evaluating at most $N$ closed-form expressions.

The strength of selection grows with increasing dilution rate $\phi$ and/or decreasing substrate feed concentration $\gamma$, causing the replicators to be driven to extinction one-by-one in order of increasing $k_{i}$. In Appendix C, we show that replicator $X_{i}$ becomes extinct once

$$\frac{\gamma}{\phi}\leq k_{i}^{-1}+\sum_{j:k_{j}\geq k_{i}}e^{-\Delta G_{j}^{\circ}}(k_{i}^{-1}-k_{j}^{-1}).$$ (41)

Figure 4: Steady state behavior of a system of 4 elementary replicators, for varying values of the dilution rate $\phi$. Top: steady-state concentrations of the four replicators. As $\phi$ increases, the replicators are driven to extinction one-by-one (dashed vertical lines). Bottom: As predicted by the bound (36), replicator $X_{i}$ are pushed to extinction once the Gibbs free energy of the fittest replicator (blue line) crosses the selection coefficient $-\ln s_{i}=-\ln(1-f_{i}/f_{1})$.

We now consider a concrete example of 4 replicators with rate constants $(k_{1},k_{2},k_{3},k_{4})=(4,3,2,1)$ and $-\Delta G_{i}^{\circ}$ values $(1,2,3,2.5)$. Using Eq. 40, we calculate steady-state concentrations $x_{i}^{*}$ of the 4 replicators at different values of the dilution rate $\phi$, while substrate feed concentration is set to $\gamma=1$. The steady-state concentrations are shown in Fig. 4 (top). As the dilution rate increases, the replicators go extinct one-by-one in order of increasing $k_{i}$. The critical values of $\phi$ at which each replicator goes extinct, as specified by Eq. 41, are indicated with dotted vertical lines.

In Fig. 4 (bottom), we show the steady-state Gibbs free energy of replication for each replicator,

$$\sigma_{i}^{*}=\ln(a^{*}/x_{i}^{*})-\Delta G_{i}^{\circ}.$$ (42)

The values of $\sigma_{i}^{*}$ grow with increasing $\phi$, diverging to infinity as each replicator approaches extinction. We compare $\sigma_{1}$, the Gibbs free energy of replication of the fittest replicator $X_{1}$, to the selection coefficient between $X_{1}$ and $X_{i}$, $s_{i}=1-f_{i}/f_{1}$. Each replicator’s fitness is $f_{i}=k_{i}a^{*}$, so $s_{i}=1-k_{i}/k_{1}$. As predicted by our bound (36), replicator $X_{i}$ becomes extinct once $\sigma_{1}^{*}$ crosses $-\ln s_{i}$.

We note that fitness values do not determine relative concentrations in steady state. For instance, near equilibrium (small dilution rates), steady-state concentrations are determined by the standard Gibbs free energies $-\Delta G_{i}^{\circ}$ rather than fitness values. This can be seen in Fig. 4 (top): replicator $X_{3}$ has the largest steady-state concentration at small $\phi$ values, since it has the largest value of $-\Delta G_{i}^{\circ}$.

We can also consider the overall rate of Gibbs free energy dissipation due to replication. In steady state, it is

$$\dot{\Sigma} =\sum_{i}\mathcal{J}_{i}\sigma_{i}^{*}=\phi\sum_{i}x_{i}^{*}[\ln(a^{*}/x_{i}^{*})-\Delta G_{i}^{\circ}].$$ (43)

Here we used Eq. 42 and that the steady-state flux across the autocatalytic reaction of replicator $X_{i}$ is $\mathcal{J}_{i}=\phi x_{i}^{*}$.

Figure 5: Black line shows $\dot{\Sigma}$, the total rate of Gibbs free energy dissipation, Eq. 43, for the 4-replicator model as a function of the dilution rate. Shaded regions indicate contributions from different replicator populations, with colors as in Fig. 4. At the four extinction events (dotted lines), the entropy production rate is continuous but not differentiable, corresponding to second-order nonequilibrium phase transitions.

In Fig. 5, we plot $\dot{\Sigma}$ for the 4-replicator system analyzed above. We also plot the contribution from each type of replicator, $\phi x_{i}^{*}[\ln(a^{*}/x_{i}^{*})-\Delta G_{i}^{\circ}]$, as a shaded region. As before, we vary the dilution rate $\phi$ while holding fixed the substrate feed concentration at $\gamma=1$. Extinction events are marked using vertical dotted lines. Recall from Fig. 4 that the Gibbs free energy of replication $\sigma^{*}_{i}$ diverges when replicator $X_{i}$ approaches extinction. However, the concentration $x_{i}^{*}$ vanishes sufficiently fast so that the product $x_{i}^{*}\sigma_{i}^{*}\to 0$ at extinction. As we show in Appendix C, $\dot{\Sigma}$ is finite and continuous at the extinction events. Thus, under a common classification scheme (schloglChemicalReactionModels1972a, ; mcneilNonequilibriumPhaseTransitions1974, ; zhangCriticalBehaviorEntropy2016, ; tomeEntropyProductionNonequilibrium2012a, ; nguyenExponentialVolumeDependence2020, ), extinction events are second-order nonequilibrium phase transitions.

VII Cross-catalytic cycles

We finish by briefly discussing how our results generalize to certain types of autocatalytic sets (kauffman1986autocatalytic, ). For simplicity, we restrict our attention to autocatalytic sets with a uniform and cyclic organization. A general treatment the thermodynamics of cross-catalytic cycles with arbitrary topologies, kinetics, and thermodynamic parameters is an important direction for future work.

We consider an autocatalytic set that contains $n$ species $(Z_{1},\dots,Z_{n})$ and $n$ reactions, where each species $Z_{j-1}$ catalyzes the formation $Z_{j}$:

$$Z_{j-1}+\sum_{i}\alpha_{i}^{(j)}A_{i}\rightleftharpoons Z_{j-1}+Z_{j}+\sum_{i}\beta_{i}^{(j)}A_{i}.$$ (44)

The indexes are taken as $\mathrm{mod\,}n$, so $Z_{0}=Z_{n}$, and $\alpha_{i}^{(j)},\beta_{i}^{(j)}$ indicate stoichiometric coefficients of substrates/side products participating in each reaction. Each catalytic reaction in the cycle may be elementary, or it may be a non-elementary mechanism as in Eq. 11 but with initial reactant $X$ replaced by $Z_{j-1}$ and final products $X+X$ replaced by $Z_{j-1}+Z_{j}$.

Figure 6: Examples of cross-catalytic cycles. Left: a 3-element cycle. Right: templated replication of complementary polymers (shown here using dimers).

We term this kind of autocatalytic set a cross-catalytic cycle. A schematic illustration of a 3-member cross-catalytic cycle is shown in Fig. 6 (left). Cross-catalytic cycles have attracted much attention in work on the origin-of-life, both theoretical (eigenSelforganizationMatterEvolution1971, ; bagleyEvolutionMetabolism1992a, ; hordijk2019history, ) and experimental (vaidyaSpontaneousNetworkFormation2012, ). An important example of a two-member cross-catalytic cycle is the templated replication of complementary polymers, illustrated in Fig. 6 (right), which has been investigated in numerous experiments (sievers1994self, ; lincolnSelfsustainedReplicationRNA2009, ; bissetteMechanismsAutocatalysis2013a, ). In biology, a cross-catalytic cycle called the “Hinshelwood cycle” has been proposed as a general model of bacterial growth (hinshelwood136ChemicalKinetics1952, ; iyer2014universality, ).

For simplicity, we assume that the cycle is uniform, in the sense that each cross-catalytic reaction has the same kinetic, stoichiometric, and thermodynamic properties. Although this assumption seems restrictive, it suffices for studying many autocatalytic sets of fundamental interest, such as complementary pairs with two kinetically and thermodynamically similar cross-catalytic reactions, as in Fig. 6 (right).

For an autocatalytic set with a uniform cyclic organization, the cycle members $Z_{j}$ approach equal concentrations $z_{j}\approx z_{k}$ after an initial transient. In this regime, we can effectively treat each cycle member $Z_{j}$ as an independent replicator with autocatalytic flux

$$\mathcal{J}_{j}=r(\rho)z_{j}-r^{-}(\rho)z_{j}^{2}.$$ (45)

Here $\rho=\mathcal{J}_{j}/z_{j}$ is the growth rate and $r(\rho),r^{-}(\rho)$ are effective forward/backward rate constants, neither of which depend on $j$ due to the assumption of uniformity. When cross-catalytic reactions are elementary, Eq. 45 follows immediately from mass-action kinetics and the assumption $z_{j}\approx z_{j-1}$. For non-elementary cross-catalytic reactions with $m$ elementary steps, we consider the different steps that produce or consume $Z_{j}$: the $m$-th step of cross-catalytic reaction $j$ produces one $Z_{j}$, the first step of cross-catalytic reaction $j+1$ consumes one $Z_{j}$, and the $m$-th step of cross-catalytic reaction $j+1$ produces one $Z_{j}$. The overall production of $Z_{j}$ due to the cross-catalytic cycle is

$$\mathcal{J}_{j}=J_{m}^{(j)}-J_{1}^{(j+1)}+J_{m}^{(j+1)},$$ (46)

where $J_{m}^{(j)}$ is the flux across the $k$-th intermediate step in cross-catalytic reaction $j$, defined similarly to Eq. 12. Assuming uniformity of kinetics and concentrations, the intermediate fluxes are the same for all $j$, $J_{i}^{(j)}=J_{i}^{(j+1)}\equiv J_{i}$, therefore $\mathcal{J}_{j}=2J_{m}-J_{1}$ which recovers Eq. 13. We can then derive Eq. 45 using the approach we developed for treating non-elementary replicators in Section III.2. The fitness $f$ of each cycle member $Z_{j}$ is defined via the relation $f=r(f)$, since it can be shown from Eq. 45 that this determines both initial growth rate and critical dilution rate.

The Gibbs free energy of each cross-catalytic reaction $j$ is $\sigma_{(j)}=\ln[r(0)/z_{j}r^{-}(0)]$, which can be shown using a similar derivation as Eq. 26. Combining these results and using the derivation from Section IV.2 gives a bound on the Gibbs free energy of each cross-catalytic reaction $j$ as

$$\sigma_{(j)}\geq-\ln\left(1-\frac{\rho}{f}\right),$$ (47)

thereby recovering the analogue of inequality Eq. 32. In steady state with dilution rate $\phi$, this reduces to $\sigma_{(j)}^{*}\geq-\ln(1-\phi/f)$, where $\sigma_{(j)}^{*}$ is the Gibbs free energy of replication under steady-state concentrations.

We can also derive a thermodynamic bound on the strength of selection for cross-catalytic cycles, analogous to Eq. 36. We consider a flow reactor in steady state with dilution rate $\phi$ that contains a cross-catalytic cycle with fitness $f>\phi$. Suppose there is another replicator $X^{\prime}$, which may be either a cross-catalytic cycle or an individual autocatalytic reaction, that has fitness $f^{\prime}\leq\phi$ and is therefore extinction in steady state. Plugging $f^{\prime}\leq\phi$ into the inequality above gives

$$\sigma_{(j)}^{*}\geq-\ln\left(1-\frac{f^{\prime}}{f}\right)=-\ln s.$$ (48)

This result bounds the average Gibbs free energy of a cross-catalytic reaction in terms of the selection coefficient, thereby generalizing the bound (36).

We emphasize that these bounds are all stated in terms of the Gibbs free energy of a single reaction in the cross-catalytic cycle. For example, for a self-replicating complementary polymer as in Fig. 6 (right), these inequalities bound the Gibbs free energy dissipated when making a single complementary copy, that is half of the overall cycle. The Gibbs free energy dissipated in carrying out all $n$ reactions, as necessary to complete the cycle, is bounded by $-n\ln(1-\rho/f)$.

VIII Discussion and future work

In this paper, we uncovered a general relationship between dissipated Gibbs free energy, fitness, and selection strength in molecular replicators. This relationship was derived from the principle of local detailed balance, which plays a central role in nonequilibrium thermodynamics. Our results complement other recent work on fundamental stoichiometric and thermodynamic constraints on autocatalytic replication (despons2023structural, ; gagrani2023polyhedral, ; srinivas2023characterizing, ; srinivas2023thermodynamics, ; blokhuisUniversalMotifsDiversity2020, ; pinero_nonequilibrium_2018, ).

Our approach can also be related to competitive exclusion and coexistence theory from microbial ecology (smithTheoryChemostatDynamics1995, ; groverResourceCompetition1997, ; ellnerExpandedModernCoexistence2019, ). Importantly, however, our underlying question differs from the one typically asked in ecology: ecologists usually aim to explain how diversity may be maintained despite competition (ellnerExpandedModernCoexistence2019, ), while we aim to explain how selection may occur (decrease of diversity) despite the drive toward thermodynamic equilibrium. In any case, our work contributes to the study of autocatalytic systems from the perspective of theoretical ecology (leeAutocatalyticNetworksTransition1997, ; vaidyaSpontaneousNetworkFormation2012, ; szilagyi2017ecology, ; peng2020ecological, ) and evolutionary (vasasEvolutionGenes2012a, ; takeuchiEvolutionaryDynamicsRNAlike2012, ; mizuuchiEvolutionaryTransitionSingle2022, ; ametaDarwinianPropertiesTheir2021, ; ameta2021self, ).

Our first inequality (32) may be compared to a bound derived by England (Eq. 10, england2013statistical, ). However, the two bounds apply to different setups, involve different thermodynamic and dynamic quantities, and are not formally equivalent. Our results concern reversible replicators with concentration-dependent fluxes which reach steady state. England’s result concerns irreversible replicators with fixed production and degradation rates which grow exponentially without bound. Our setup may be better suited for studying thermodynamics of Darwinian evolution, since exponentially growing replicators always coexist and do not exhibit selection as traditionally understood (eigenSelforganizationMatterEvolution1971, ; bernstein1983darwinian, ). A deeper comparison of the two approaches remains for future work.

The Gibbs free energy of replication $\sigma$ plays a central role in our results. We emphasize that this quantity depends both on the intrinsic chemical properties of the replicator and on its concentration $x$, which can change over time due to reactions, dilution or inflows, etc. For a replicator in steady state, we may express these two contributions by rewriting Eq. 26 as $\sigma^{*}=\ln[r(0)/r^{-}(0)]-\ln x^{*}$. The intrinsic contribution $\ln[r(0)/r^{-}(0)]$ depends on the replicator’s standard Gibbs free energy $-\Delta G^{\circ}$ and stoichiometric coefficients $\alpha_{k,i},\beta_{k,i}$, as well as concentrations of substrates/side products $\vec{a}$ in the reaction volume. This contribution can be increased by using richer substrates (changing stoichiometric coefficients $\alpha_{k,i},\beta_{k,i}$), by enriching existing substrates (increasing $-\Delta G^{\circ}$), as well as in other ways. Such changes to the intrinsic contribution may increase or decrease $\sigma^{*}$, depending on how they affect the steady-state concentrations of the replicator and substrates/side products. They may also increase, decrease, or leave the fitness $f$ unchanged, depending on how they affect the effective rate constant of the replication mechanism. In summary, our bounds concern the overall Gibbs free energy of replication, which combines both intrinsic and concentration-dependent contributions. For this reason, they do not directly imply that there is a universal advantage for exploiting richer substrates or otherwise increasing the intrinsic free energy dissipation.

We mention some limitations of our results and directions for future work.

Like many other results derived using local detailed balance, our bound is most meaningful for molecular replicators that are not “too irreversible”, meaning that $\sigma$ is not too large. For instance, the bound (32) can be rearranged as $\rho\leq f(1-e^{-\sigma})$, which reduces to the trivial result $\rho\leq f$ once $\sigma$ grows beyond 20 or so (approximately the dissipation produced by ATP hydrolysis). Our bound on the selection coefficient, which appeared as $s\geq e^{-\sigma^{*}}$ in (2), also becomes weak in this regime. Interestingly, however, our bound (47) for cross-catalytic cycles refers to the Gibbs free energy of a single reaction in the cycle. This result may be applicable to larger-scale systems that involve large cycles consisting of many reversible cross-catalytic reactions.

We only considered deterministic concentrations, which is justified when molecular counts are large and stochastic fluctuations can be ignored. However, fluctuations cannot be ignored in small systems, nor near extinction events when concentrations are small (mcneilNonequilibriumPhaseTransitions1974, ; nitzanFluctuationsTransitionsChemical1974, ). Future work may extend our analysis to the stochastic regime. It is also interesting to investigate the assumption of ideal solution, which may be violated at sufficiently high concentrations.

We did not consider the effect of mutations. In general, mutations weaken the strength of selection (eigenSelforganizationMatterEvolution1971, ), therefore we expect that mutations can only increase the thermodynamic costs of Darwinian evolution. Future work may verify this prediction and seek stronger bounds on the thermodynamic cost of Darwinian evolution for replicators with mutations. The introduction of mutations leads to other important questions concerning the thermodynamics of evolution, such as the thermodynamic costs of finding new high-fitness replicators rather than merely selecting among existing replicators. This might be considered the thermodynamics of “the arrival of the fittest”, rather than of “the survival of the fittest” (fontanaArrivalFittestTheory1994, ; wagnerArrivalFittestHow2015, ).

Finally, our analysis of autocatalytic sets was restricted to the case where reactions are organized in a single uniform cycle. Future work may consider autocatalytic sets with more general topologies, kinetics, and thermodynamics (jain1998autocatalytic, ; blokhuisUniversalMotifsDiversity2020, ). Similarly, our analysis of multi-step reaction mechanisms was restricted to linear sequences of reactions such as Eq. 11, which may be generalized in future work.

APPENDIX A Degradation reactions

Here we show that our thermodynamic bounds can be generalized to account for degradation reactions.

Let us consider the case where the replicator $X$ undergoes degradation. In the absence of degradation, the replicator’s concentration changes as

$$\dot{x}=\mathcal{J}-\phi x=r(\rho)x-r^{-}(\rho)x^{2}-\phi x,$$ (49)

where $\mathcal{J}=r(\rho)x-r^{-}(\rho)x^{2}$ is the net flux across the autocatalytic reaction. Now suppose that the replicator undergoes a degradation reaction,

$$X+\sum_{i}\alpha_{i}^{\prime}A_{i}\to\sum_{i}\beta_{i}^{\prime}A_{i}\,,$$ (50)

where $\alpha_{i}^{\prime}$ and $\beta_{i}^{\prime}$ are stoichiometric coefficients. We assume that the degradation is irreversible and that the flux across the degradation reaction is first order with rate constant $\eta\geq 0$, so $\eta x$. As usual, the rate constant $\eta$ may depend on the concentrations of substrates/side products $\vec{a}$, though we leave this implicit.

In the presence of degradation, the replicator’s concentration changes as

$$\dot{x}=\mathcal{J}-\phi x=r(\rho)x-r^{-}(\rho)x^{2}-\eta x-\phi x.$$ (51)

We assume that $\eta<r(\rho)$, since otherwise positive growth is impossible. We then define a modified forward rate constant,

$$\tilde{r}(\rho):=r(\rho)-\eta\geq 0,$$ (52)

so that Eq. 51 becomes

$$\dot{x}=\tilde{r}(\rho)x-r^{-}(\rho)x^{2}-\phi x,$$ (53)

thereby recovering the form of Eq. 49. The flux-force relation (26) becomes an inequality,

$$\sigma=\ln\frac{r(0)}{xr^{-}(0)}\geq\ln\frac{\tilde{r}(0)}{xr^{-}(0)},$$ (54)

since $\tilde{r}(0)\leq r(0)$. We can then define fitness $\tilde{f}$ in the same way as Eq. 27, except using $\tilde{r}(\rho)$ rather than $r(\rho)$. The derivation of our bounds in Section IV.2 proceeds in the same manner as before using these new definitions. Note that the derivative relations in Eq. 21 still apply to $\tilde{r}(\rho)$, since $\eta$ does not depend on $\rho$.

Let us now consider the case of a non-elementary replicator $X$ where some intermediate species $Y_{k}$ also undergo irreversible degradation with rate constants $\pi_{k}\in[0,\nu_{k})$. To account for this, the matrix $\bar{M}$ in the linear system (17) turns into

$$\bar{M}_{ij}:=M_{ij}+\delta_{ij}\pi_{i}.$$ (55)

Matrix calculus shows that

$$\frac{\partial\bar{M}_{jk}^{-1}}{\partial\bar{M}_{ii}}=-\bar{M}_{ji}^{-1}\bar{M}_{ik}^{-1}.$$

For any set of degradation rates $\pi_{i}$, $\bar{M}$ is an M-matrix so it is invertible and all of entries of its inverse are nonnegative (shivakumar1996two, ; meyer1978singular, ). Therefore, every entry of $\bar{M}^{-1}$ decreases with increasing $\pi_{k}$. This means that the effective rate constants $\bar{r}(\rho)$ and $\bar{r}^{-}(\rho)$, defined as in Eq. 20 using $\bar{M}$, obey

$$\bar{r}(\rho)\leq r(\rho)\qquad\bar{r}^{-}(\rho)\geq r^{-}(\rho).$$

To account for the degradation of the replicator itself with rate $\eta$, we can further modify the effective forward rate constant as

$$\tilde{\bar{r}}(\rho):=\bar{r}(\rho)-\eta,$$

similarly to Eq. 52. As above, we have a flux-force inequality,

$$\sigma=\ln\frac{r(0)}{xr^{-}(0)} \geq\ln\frac{\bar{r}(\rho)}{x\bar{r}^{-}(\rho)}\geq\ln\frac{\tilde{\bar{r}}(\rho)}{x\bar{r}^{-}(\rho)}.$$

The derivation of our bounds then proceeds as before, after defining fitness $\tilde{\bar{f}}$ using the forward rate constant $\tilde{\bar{r}}(\rho)$. Note that the derivative relations in Eq. 21 still apply to $\tilde{\bar{r}}(\rho)$, since $\bar{M}$ is an M-matrix and $\eta$ does not depend on $\rho$.

In general, our thermodynamic bounds will tend to be less tight in the presence of degradation, since the flux-force equality (26) turns into an inequality. At the same time, the Gibbs free energy of replication may increase due to degradation, for instance because the steady-state concentration of replicator decreases due to degradation. Some of our bounds may also be larger, for instance the right side of the steady-state bound (33) will be larger, since fitness will decrease due to degradation.

APPENDIX B Non-elementary replicators

Here we derive closed-form expressions of the effective rate constants in the $\phi=0$ regime, as appear in Eq. 22. As discussed at the end of Section III.2, this regime applies when the internal reactions are much faster than the dilution rate.

Plugging $\phi=0$ into Eq. 16 implies that $J_{k}=J_{k+1}$, so all intermediate fluxes are equal:

$$J_{k}=J\qquad\text{for }k\in{1..m}.$$ (56)

They also all equal to the rate of production of the replicator, $\mathcal{J}=2J_{m}-J_{1}=J$. We can then rearrange Eq. 12 to write the steady-state concentrations of $Y_{k}$ as

$$y_{k}=(\nu_{k}/\nu_{k}^{-})y_{k-1}-\mathcal{J}/\nu_{k}^{-}\qquad k\in{1..m},$$

where we use the convention $y_{0}=x$, $y_{m}=x^{2}$. This is a first-order linear recurrence relation for $y_{k}$. Using an existing result (Thm. 7.1, charalambidesEnumerativeCombinatorics2002, ), this recurrence can be solved for $y_{m}=x^{2}$ starting from the initial condition $y_{0}=x$ to give

$$x^{2} =\left(x-\mathcal{J}\sum_{k=1}^{m}\frac{\prod_{l=1}^{k-1}\nu_{l}^{-}}{\prod_{l=1}^{k}\nu_{l}}\right)\prod_{k=1}^{m}\frac{\nu_{k}}{\nu_{k}^{-}}$$ (57)

where we use the convention $\prod_{l=1}^{k-1}\nu_{l}^{-}=1$ for $k=1$. Plugging in $r(0)$ and $r^{-}(0)$ from Eq. 22 gives

$$x^{2}=\left(x-\mathcal{J}/r(0)\right)\frac{r(0)}{r^{-}(0)}.$$

This can be rearranged as $\mathcal{J}=r(0)x-r^{-}(0)x^{2}$, which recovers the mass-action-like form of Eq. 19.

APPENDIX C Chemostat model

C.1 Steady state

Here we analyze the steady-state behavior of the dynamical system described by Eq. 39. To begin, let $\omega:=a+\sum_{i}x_{i}$ indicate the total concentration of substrate and replicators at time $t$. The first line of Eq. 39 means that $k_{i}x_{i}[a-e^{\Delta G_{i}^{\circ}}x_{i}]=\dot{x}_{i}+\phi x_{i}$. Plugging this into the second line of Eq. 39 and rearranging gives

$$\dot{\omega}=\phi(\gamma-\omega).$$

Thus, $\omega$ converges exponentially fast to the steady-state value $\omega^{*}=\gamma$.

We consider the long-term dynamics of the system restricted to the invariant subspace $\omega=\gamma$. Within this subspace, we can rewrite the first line of Eq. 39 as

$$\dot{x}_{i}=k_{i}x_{i}[\gamma-\sum_{j}x_{j}-e^{\Delta G_{i}^{\circ}}x_{i}]-\phi x_{i}.$$ (58)

Using an appropriate Lyapunov function, Schuster and Sigmund demonstrated that the dynamics in Eq. 58 converge to the steady state in Eq. 40 for any strictly positive initial condition $(x_{1}(0),\dots,x_{N}(0))\in\mathbb{R}_{+}^{N}$ (schusterDynamicsEvolutionaryOptimization1985, ). A similar global convergence result can also be derived from the theory of Lotka-Volterra dynamics (takeuchi_existence_1980, ). Specifically, Eq. 58 can be written as a competitive Lotka-Volterra system,

$$\dot{x}=b_{i}x_{i}+\sum_{j}R_{ij}x_{i}x_{j},$$ (59)

where $b_{i}:=k_{i}\gamma-\phi$ and $R_{ij}=-k_{i}(1+\delta_{ij}e^{\Delta G_{i}^{\circ}})$. The matrix $R$ can be expressed as $R=-K(\vec{1}\vec{1}^{T}+D)$, where $K_{ij}=\delta_{ij}k_{i}$ and $D=\delta_{ij}e^{\Delta G_{i}^{\circ}}$ are diagonal matrices and $\vec{1}$ is a vector of all 1s. Note that $\vec{1}\vec{1}^{T}+D$ is positive definite, since $\vec{1}\vec{1}^{T}$ is positive semidefinite and $D$ is positive definite. It is known that for this type of Lotka-Volterra system, any strictly positive initial condition converges to a unique globally attracting fixed point (takeuchi_existence_1980, ), which is the steady state specified by Eq. 40.

The steady state in Eq. 40 is expressed as a set of coupled equations, which can be solved in the following manner. Assume without loss of generality that the rate constants $k_{i}$ are arranged in decreasing order, $k_{1}\geq k_{2}\geq\dots\geq k_{N}$. Given Eq. 40, it must then be that $x_{i}^{*}=0$ implies $x_{j}^{*}=0$ whenever $j>i$. Suppose for the moment that the top $i\in\{0..N\}$ replicators have non-zero steady-state concentrations,

$$x^{*}_{j}=\begin{cases}e^{-\Delta G_{j}^{\circ}}(a^{*}-\phi/k_{j})&j\leq i\\ 0&j>i\end{cases}$$ (60)

Eq. 40 then gives $a^{*}=\gamma-\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}(a^{*}-\phi/k_{j})$, so

$$a^{*} =\frac{\gamma+\phi\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}k_{j}^{-1}}{1+\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}}.$$ (61)

Eqs. 60 and 61 are the solution to Eq. 40 if for all $j\in\{0..N\}$,

$$x_{j}^{*}=\max\{0,e^{-\Delta G_{j}^{\circ}}(a^{*}-\phi/k_{j})\}.$$ (62)

Given Eq. 60, Eq. 62 is satisfied once

$$a^{*}-\phi/k_{i}\geq 0\geq a^{*}-\phi/k_{i+1},$$ (63)

Therefore, to solve Eq. 40, it suffices to evaluate Eqs. 61 and 60 for $i=0,1,2,\dots$, stopping once the bounds (63) are satisfied.

C.2 Derivation of inequality (41) (condition for extinction)

Given Eq. 40, replicator $X_{i}$ is extinct once

$$a^{*}\leq\phi/k_{i}.$$ (64)

If this inequality holds, then any lower fitness replicator $X_{j}$ ($k_{j}\leq k_{i}$) must also be extinct, since then $a\leq\phi/k_{j}$. Now, combine the equations in Eq. 40 to write

$$a^{*} =\gamma-\sum_{j:x_{j}>0}e^{-\Delta G_{j}^{\circ}}(a^{*}-\phi/k_{j})$$

$$\leq\gamma-\sum_{j:k_{j}>k_{i}}e^{-\Delta G_{j}^{\circ}}(a^{*}-\phi/k_{j}),$$ (65)

where the inequality in the second line reflects that it may be that $a^{*}\leq\phi/k_{j}$ even for higher fitness replicators ($k_{j}>k_{i}$). Rearranging inequality (65) gives

$$a\leq\frac{\gamma+\phi\sum_{j:k_{j}\geq k_{i}}e^{-\Delta G_{j}^{\circ}}k_{j}^{-1}}{1+\sum_{j:k_{j}\geq k_{i}}e^{-\Delta G_{j}^{\circ}}}.$$ (66)

Given inequality (66), the bound (64) must be satisfied when

$$\frac{\gamma+\phi\sum_{j:k_{j}\geq k_{i}}e^{-\Delta G_{j}^{\circ}}k_{j}^{-1}}{1+\sum_{j:k_{j}\geq k_{i}}e^{-\Delta G_{j}^{\circ}}}\leq\phi/k_{i}.$$

Rearranging gives the inequality (41).

C.3 Nonequilibrium phase transitions at extinctions

Here we show that the overall rate of Gibbs free energy dissipation, Eq. 43, is continuous but not differentiable at extinction events, meaning that extinctions are second-order nonequilibrium phase transitions.

Suppose that the rate constants are strictly ordered as

$$k_{1}>k_{2}>\dots>k_{N}.$$ (67)

From Eq. 41, the critical dilution rate for replicator $X_{i}$ is

$$\hat{\phi}_{(i)} =\frac{\gamma}{k_{i}^{-1}+\sum_{j=1}^{i-1}e^{-\Delta G_{j}^{\circ}}(k_{i}^{-1}-k_{j}^{-1})}$$ (68)

$$=\frac{\gamma}{k_{i}^{-1}+\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}(k_{i}^{-1}-k_{j}^{-1})}$$ (69)

where in the second line we used $k_{i}^{-1}-k_{j}^{-1}=0$ for $j=i$. We treat $\phi$ as the control parameter, while holding $\gamma$ fixed.

We show that the steady-state substrate concentration $a^{*}$ is continuous as a function of the dilution rate at $\hat{\phi}_{(i)}$. When $\phi<\hat{\phi}_{(i)}$, replicators $X_{1},\dots,X_{i}$ are not extinct. We consider the limit from below,

$$\lim_{\phi\nearrow\hat{\phi}_{(i)}}a^{*} =\frac{\gamma+\hat{\phi}_{(i)}\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}k_{j}^{-1}}{1+\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}}=\hat{\phi}_{(i)}k_{i}^{-1}.$$

Here we first used Eq. 61 and then plugged in Eq. 68 and simplified using a bit of tedious algebra. When $\phi>\hat{\phi}_{(i)}$, replicators $X_{1},\dots,X_{i-1}$ are not extinct, but replicator $X_{i}$ is extinct. We consider the limit from above,

$$\lim_{\phi\searrow\hat{\phi}_{(i)}}a^{*} =\frac{\gamma+\phi\sum_{j=1}^{i-1}e^{-\Delta G_{j}^{\circ}}k_{j}^{-1}}{1+\sum_{j=1}^{i-1}e^{-\Delta G_{j}^{\circ}}}=\hat{\phi}_{(i)}k_{i}^{-1},$$

where we first used Eq. 61 and then plugged in Eq. 69 and simplified. The two limits match, so $a^{*}$ is continuous at $\hat{\phi}_{(i)}$. This also implies that steady-state replicator concentrations $x_{j}^{*}=\max\{0,e^{-\Delta G_{j}^{\circ}}(a^{*}-\phi/k_{j})\}$ from Eq. 40 are continuous at $\hat{\phi}_{(i)}$.

Next, consider the rate of Gibbs free energy dissipation at the critical point,

$$\lim_{\phi\to\hat{\phi}_{(i)}}\dot{\Sigma}=\\ \lim_{\phi\to\hat{\phi}_{(i)}}\phi\Bigg{[}\sum_{j=1}^{i-1}x_{j}^{*}\Big{(}\ln\frac{a^{*}}{x_{j}^{*}}-\Delta G_{j}^{\circ}\Big{)}+x_{i}^{*}(\ln a^{*}-\Delta G_{i}^{\circ})\Bigg{]},$$

where we used $\lim_{\phi\to\hat{\phi}_{(i)}}x_{i}^{*}\ln x_{i}^{*}=\lim_{\alpha\to 0}\alpha\ln\alpha=0.$ Given Eq. 67, $x_{j}^{*}>0$ for $j\in\{1..i-1\}$, so all the terms on the right side are finite and continuous at $\phi=\hat{\phi}_{(i)}$. Thus, $\dot{\Sigma}$ is a continuous function of $\phi$ at $\hat{\phi}_{(i)}$.

Next, we show that $\dot{\Sigma}$ is not differentiable with respect to $\phi$ at $\hat{\phi}_{(i)}$ because it has an infinite left derivative at this point (it is not hard to show that the right derivative is finite). The left derivative is

$$\partial_{\phi}^{-}\dot{\Sigma} =\sum_{j=1}^{i}\left[\partial_{\phi}^{-}x_{j}^{*}(\ln a^{*}-\ln x_{j}^{*}-1-\Delta G_{j}^{\circ})+\frac{x_{j}^{*}}{a^{*}}\partial_{\phi}^{-}a^{*}\right].$$

All terms in this expression are finite except for $-(\partial_{\phi}^{-}x_{i}^{*})\ln x_{i}^{*}$. Plugging in Eq. 40 gives

$$\partial_{\phi}^{-}x_{i}^{*} =e^{-\Delta G_{i}^{\circ}}(\partial_{\phi}^{-}a^{*}-k_{i}^{-1})$$

$$=e^{-\Delta G_{i}^{\circ}}\left[\frac{\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}k_{j}^{-1}}{1+\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}}\hat{\phi}_{(i)}k_{i}^{-1}-k_{i}^{-1}\right]$$

$$=e^{-\Delta G_{i}^{\circ}}\left[\frac{\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}(k_{j}^{-1}-k_{i}^{-1})-k_{i}^{-1}}{1+\sum_{j=1}^{i}e^{-\Delta G_{j}^{\circ}}}\right]<0,$$

where in the second line we used Eq. 61, and in the last line we used that $k_{j}>k_{i}$ and $k_{i}>0$. Hence $\partial_{\phi}^{-}x_{i}^{*}$ is a strictly negative constant, meaning that $\phi\nearrow\hat{\phi}_{(i)}$, and $x_{i}^{*}\to 0$, $-(\partial_{\phi}^{-}x_{i}^{*})\ln x_{i}^{*}\to-\infty$. This shows that the left derivative of $\dot{\Sigma}$ is negative infinite at $\phi=\hat{\phi}_{(i)}$.