#  Robustness of compositional heredity to the growth and division dynamics of
prebiotic compartments

Yoshiya J. Matsubara<EMAIL_ADDRESS>Sandeep Ameta  Simons Centre for
the Study of Living Machines, National Centre for Biological Sciences, Tata
Institute of Fundamental Research, Bengaluru, India  Shashi Thutupalli  Simons
Centre for the Study of Living Machines, National Centre for Biological
Sciences, Tata Institute of Fundamental Research, Bengaluru, India
International Centre for Theoretical Sciences, Tata Institute of Fundamental
Research, Bengaluru, India  Philippe Nghe  Laboratory of Biophysics and
Evolution, UMR CNRS-ESPCI Chemistry Biology Innovation, École Supérieure de
Physique et de Chimie Industrielles de la Ville de Paris (ESPCI Paris), PSL
University, France  Sandeep Krishna  Simons Centre for the Study of Living
Machines, National Centre for Biological Sciences, Tata Institute of
Fundamental Research, Bengaluru, India

###### Abstract

An important transition after the origin of life was the first emergence of a
Darwinian population, self-reproducing entities exhibiting differential
reproduction, phenotypic variation, and inheritance of phenotypic traits. The
simplest system we can imagine to have these properties would consist of a
compartmentalized autocatalytic reaction system that exhibits two growth
states with different chemical compositions. Identifying the chemical
composition as the phenotype, this accounts for two of the properties.
However, it is not clear what are the necessary conditions for such a chemical
system to exhibit inheritance of the compositional states upon growth and
division of the compartment. We show that for a general class of autocatalytic
chemical systems subject to serial dilution, the inheritance of compositional
information only occurs when the time interval between dilutions is below a
critical threshold that depends on the efficiency of the catalytic reactions.
Further, we show that these thresholds provide rigorous bounds on the
properties required for the inheritance of the chemical compositional state
for general growth and division cycles. Our result suggests that a serial
dilution experiment, which is much easier to set up in a laboratory, can be
used to test whether a given autocatalytic chemical system can exhibit
heredity. Lastly, we apply our results to a realistic autocatalytic system
based on the _Azoarcus_ ribozyme and suggest a protocol to experimentally test
whether this system can exhibit heredity.

bistability, autocatalytic sets, compartmentalization, growth and division,
Darwinian evolution, heredity, serial dilution, continuous stirred-tank
reactor

##  I Introduction

The origin of life on Earth involved the transition from a primordial “soup”
of chemicals to a population of self-reproducing individuals evolving under
natural selection. There have been a number of proposals for the nature of the
earliest self-reproducing entities, ranging from RNA [1, 2] to clay [3] to X.
Even within broadly supported frameworks like the RNA world [4], there are
many possibilities – some have advocated for collective autocatalytic sets [5,
6, 7, 8] while others search for the simplest RNA ribozyme that can copy
itself [9]. There are also multiple proposals for a compartment that separates
a self-reproducing individual from its environment and other individuals [10],
ranging from lipid membranes [11] to coacervates [12] to spatial separation on
surfaces [13, 14], or within hydrodynamic flows [15]. Regardless of the
details, one could describe the emergent population as consisting of
autocatalytic chemical entities confined within compartments that grow and
divide to produce new offspring individuals. We use “autocatalytic chemical
entities” and “compartments” very broadly to encompass all the possibilities
described above. For such a system to be considered a population evolving
under natural selection [16], it must have certain properties. Following
Godfrey-Smith [17], who analyzes and builds on formulations by Lewontin [18],
Endler [19], Ridley [20], and others, we expect the individuals to exhibit:

  1. 1.

Phenotypic variation

  2. 2.

Differential reproduction

  3. 3.

Inheritance of phenotypic traits (“heredity”)

We ask: if one is provided an autocatalytic chemical system that could be
confined to compartments that grow and divide as the system consumes food
molecules to self-reproduce from the environment, how can one test whether it
has all the above properties? One of the simplest dynamical (chemical) systems
that have the properties enumerated by Godfrey-Smith consists of a bistable
chemical system, which exhibits two growth states with different chemical
compositions and, in general, different growth rates. We define a growth state
as one where the concentrations of the chemicals comprising the system grow
without bound (often exponentially) as the system consumes food molecules, but
where the chemical composition – the relative concentrations of the chemical
components – reaches a steady state. Identifying the chemical composition to
be the individual’s phenotype accounts for phenotypic variation, and assuming
the chemical composition affects the growth and division process accounts for
differential reproduction 111This could occur via a variety of mechanisms. For
example, osmotic pressure due to the difference in the composition between the
inside and outside of the compartment may induce its growth [11]. The ACS
could also produce the precursors of the compartment (e.g., lipid molecules)
[69, 70]. Another possibility is that the ACS energetically drives the growth,
and shape instability triggers the division of compartments [12].. In the
language of dynamical systems, the third property of inheritance translates to
the stability of the two growth states (hence our term ‘bistable chemical
system’) (see Fig. 1). Here, by stability, we mean that the system does not
spontaneously transition from one growth state to the other _when it divides
into two offspring entities_ 222In dynamical systems, stability often refers
to stability against random stochasticity, e.g., due to thermal noise. It is
true that for an autocatalytic chemical system to exhibit an inheritance of
its phenotype (its chemical composition), it must also exhibit a certain
amount of stability against noise. Some small enough probability of
transitions to different states due to noise can be subsumed under phenotypic
variation (indeed, this may be the only source of variation available), but
too much will destroy the property of the heredity of states. Later we provide
some results from stochastic simulations, but we largely assume that if the
system, in the absence of noise, is stable upon division, then it satisfies
the third property of inheritance of phenotypic traits..

The inheritance of the compositional information has been previously debated
in a number of models of autocatalytic networks [23, 24, 25, 26, 27]. One
suggestion has been that the network must contain multiple ‘viable
autocatalytic cores’ [25, 28] in order to exhibit heredity. Another suggestion
is that the chemical system could show multiple compositional states as the
response to environmental changes (e.g., change of the food set molecules)
[29, 30, 31]. However, it is unclear whether such mechanisms could be stable
enough to stochastic noise or environmental fluctuations, let alone to the
growth and division dynamics of compartments, which would be needed for the
inheritance of information across generations.

Multistability in chemical systems has been extensively studied in the context
of epigenetic memory in gene regulatory networks [32], signaling pathways
[33], metabolic networks [34], chiral symmetry breaking [35, 36], or enzymatic
cascades [37, 38, 39, 40, 41], etc. Generally speaking, positive feedback in
the network structure is necessary though not sufficient for multistability
[42], and more detailed conditions have also been suggested in specific
contexts [43] 333The necessary conditions for bistable chemical reactions are:
(i) positive feedback (e.g., autocatalysis), (ii) filtering noise, and (iii)
preventing explosion (e.g., conservation law of the components) [43]. In
addition, ‘nonlinearity (or ‘ultrasensitivity’) is required in the positive
feedback [33]. Later, we discuss the minimum autocatalytic chemical system
that satisfies the above condition.. However, these conditions for
bi/multistability have been investigated mainly in chemostat or continuously-
stirred-tank-reactor (CSTR) scenarios where there is a constant influx and
outflux [45]. It has not systematically been investigated under the conditions
where these reaction systems are enclosed within a compartment that
dynamically grows and divides.

Thus, our question narrows to: if one is provided an autocatalytic chemical
system capable of exhibiting two distinct growth states, what are the
necessary conditions to ensure that these states will be stable upon growth
and division of the compartments that enclose these entities? We study this
question theoretically using deterministic dynamical models of a class of
autocatalytic chemical systems that exhibit two exponential growth states. We
analyze the stability of these states when subjected to the additional growth
and division process. We first show that the dynamics of a general growth and
division process can range from a chemostat-like (CSTR) condition, where the
rate equations are subject to influx and outflux of molecules at a constant
rate, to sudden periodic reductions of the chemical concentrations. A
formulation that models growth and division as a serial dilution protocol [46]
captures these extremes and everything in between. We find the inheritance of
the compositional state under serial dilution only when the time interval
between dilutions is below a critical threshold that depends on the efficiency
of the catalytic reactions. A similar critical threshold is found for the
average dilution rate. These thresholds provide bounds on the properties of
general growth and division cycles which, if violated, will destroy the
heredity of the chemical compositional state. Moreover, because serial
dilution is a protocol much more easily implemented in the laboratory than the
growth and division of compartments, this provides a relatively easy
experimental test for whether a given autocatalytic chemical system is capable
of exhibiting heredity. We illustrate our results with models of a realistic
instance of such an autocatalytic chemical system, based on a ribozyme from
the _Azoarcus_ bacterium engineered to self-reproduce from smaller their
fragments [47, 48, 30].

Figure 1: A sketch of the schematic figures for explaining compartments with
and without heredity. Red, blue, and grey circles represent variations of
reproducing compartments with different compositions within them under the
same environment. (Top) reproduction of compartments with heredity: daughters
have the same composition as their parents. (Bottom) Heredity can break down
in two ways: either, the chemical composition converges to the same single
grey state irrespective of the initial state, or spontaneous transitions
between the two growth states, red and blue, are so rapid that all information
about the initial state is lost very quickly.

##  II Results

###  II.1 Mapping general growth and division scenarios to the serial dilution
protocol

First, we discuss how compartments’ growth and division (GD) affect the
chemical reaction systems encapsulated within them. As compartments grow and
divide, the chemical concentrations within them change due to the chemical
reactions occurring, but also get diluted due to increases in the compartment
volume. If the compartment volume grows between division as $V(t)$, then the
dilution rate is $\phi(t)\equiv\frac{dV}{dt}(t)/V(t)$. An influx of substrates
from outside the compartment may also occur, increasing the amount of those
components. Thus, the chemical rate equations must include terms for such
influx and dilution. Initially, we consider three assumptions on such a GD
process, some of which will be relaxed in later sections:

  1. 1.

Compartments divide into $m$ equal-sized smaller ones periodically, at time
intervals of $\varDelta t$ (see Fig. 2A).

  2. 2.

We assume a well-mixed condition inside a compartment. This implies that the
_concentrations_ do not change at divisions because the chemical components
are partitioned proportionally to the volumes of the daughter compartments.

  3. 3.

For simplicity, we assume the influx rate of substrates $\sigma(t)$ is
proportional to the dilution rate $\phi(t)$ (this is not an important
assumption; if $\sigma(t)$ is an arbitrary function with periodicity
$\varDelta t$ our results do not qualitatively change).

Typical dilution protocols that are used in laboratories, such as serial
dilution (SD) or the continuous stirred-tank reactor (CSTR), are special cases
of the compartment GD scenario: If $\phi(t)$ (and $\sigma(t)$) have sharp
spikes at times $t=n\varDelta t\quad(n=1,2,3,\dots)$, it corresponds to SD, in
which after each time interval $\varDelta t$ the chemical compositions are
diluted by $m=e^{\varDelta t\bar{\phi}}$ fold, and $s^{tot}(1-\frac{1}{m})$
substrate is added at the beginning of the next cycle (where $\bar{\phi}$ and
$\bar{\sigma}$ are the average dilution and influx rates over one division
cycle; see Methods and Models for details). Similarly, a CSTR corresponds to
an extreme case of a GD process where influx and dilution rates are constants:
$\phi(t)=\bar{\phi}$ and $\sigma(t)=\bar{\sigma}$. This also corresponds to an
SD protocol with an infinitesimal short interval of cycles, $\varDelta
t\rightarrow 0$ (see Methods and Models for details).

As the general case of GD protocols lies between the SD and CSTR protocols
(see Fig. 2B), and CSTR is a special case of the SD, in the following
sections, we will investigate the reaction dynamics of competing autocatalytic
reaction sets (ACSs) under the SD protocol. Later we will return to the GD
protocols and show that the condition for the bistability of the ACSs is
bounded by the conditions obtained for the SD and the CSTR protocols.

Figure 2: (A) Schematic diagram of dilution protocols. Growth and Division
(GD) cycles of a compartment with volume $V(t)$: After a period $\varDelta t$
during which the compartment grows, it divides into $m$ smaller compartments
each with volume $m=e^{\bar{\phi}\varDelta t}$ fold of the initial volume.
Continuous Stirred-Tank Reactor (CSTR) protocol: a substrate is supplied, and
compositions are diluted with the same constant rate $\bar{\phi}$. Serial
Dilution (SD) protocol: each interval $\varDelta t$, compositions are diluted
with the factor $m$. (B) The time course of the total dilution inside the
compartment undergoes (i.e., $\exp(\int_{0}^{t}\phi(t^{\prime})dt^{\prime})$).
The black solid and dotted lines correspond to the serial dilution (SD) and
continuous stirred-tank reactor (CSTR) protocols, respectively. The blue line
corresponds to a case where the compartment grows at rate
$\frac{dV}{dt}\propto V^{\alpha}$ ($\alpha=4$) and splits into $m$ equal-sized
daughters when it divides (see Material and Methods).

###  II.2 Nonlineaerlities are essential for two competing autocatalytic
entities to exhibit heredity under serial dilution

First, we test whether an autocatalytic chemical reaction system can exhibit
bistability under the SD protocol. We consider a simple class of autocatalytic
reaction systems consisting of two identical (but distinguishable)
autocatalytic entities, ${\rm{X}}_{1}$ and ${\rm{X}}_{2}$, which consume the
same substrate S (see Fig. 3A). This system can exhibit two distinct growth
states. The rate equations for this class of systems are:

$\frac{dx_{i}}{dt}=sr(x_{i})x_{i},$ (1)

where $i=1,2$, and $s$ is the concentration of the substrate S. $r(x_{i})$ is
the reproduction rate of ${\rm{X}}_{i}$ as a function of $x_{i}$ 444Here, we
consider the reproduction rate functions are symmetric between ${\rm{X}}_{1}$
and ${\rm{X}}_{2}$ for the sake of simplicity. Later we also consider the
asymmetric case in the Sec. II E.. We assume $r(x)$ is a differentiable and
non-negative function for $x\geq 0$, but otherwise place no restrictions on
its form. Later we will discuss realistic possibilities for $r(x)$.

We consider this reaction system under the SD protocol with cycle interval
$\varDelta t$ and the dilution factor $m$ ($=e^{\bar{\phi}\varDelta t}$) (see
Methods and Models). The assumptions on the protocol from the previous section
imply that the total concentration of the components $s^{tot}=s+x_{1}+x_{2}$
remains constant at $s^{tot}=\frac{\bar{\sigma}}{\bar{\phi}}$.

In the long timescale, after sufficiently many SD cycles, the trajectory of
the chemical composition reaches a stationary periodic orbit (Fig. 3B). If the
system does _not_ exhibit inheritance, then it must settle into the same
stable trajectory for every initial condition. As we assumed two ACSs which
are completely symmetric, this trajectory must be one where the concentrations
are equal (i.e., $x_{1}=x_{2}$). In contrast, if there is the inheritance of
the compositional state, then across different initial conditions, the system
must exhibit bistability, i.e., two stable trajectories. Again, due to
symmetry, in each of these two trajectories, one of the components,
${\rm{X}}_{1}$ or ${\rm{X}}_{2}$, must be dominant. Therefore, a sufficient
condition for bistability under the SD protocol can be obtained by showing the
_instability_ of the symmetrical ($x_{1}=x_{2}$) trajectory.

Introducing the notation $\chi=x_{1}+x_{2}$ and $\delta=x_{1}-x_{2}$, for
trajectories close to the symmetrical one we can assume that $\delta\ll\chi$.
We can then derive the following relation (see details of the derivation in
Appendix Sec. A.1):

$\frac{\delta(t)}{\chi(t)}=\frac{r(\frac{\chi(t)}{2})}{r(\frac{\chi(0)}{2})}\frac{\delta(0)}{\chi(0)}.$
(2)

If $\delta(t)/\chi(t)$ at the end of a cycle, $\delta(\varDelta
t)/\chi(\varDelta t)$, is larger than that at the beginning,
$\delta(0)/\chi(0)$, the trajectory is unstable, otherwise it is stable.
Therefore, the sufficient condition for bistability under the SD protocol is

$r\big{(}\frac{\chi(\varDelta t)}{2}\big{)}>r\big{(}\frac{\chi(0)}{2}\big{)}.$
(3)

That is, the stability of the compositional trajectory is determined by
whether the reproduction rate at the end of a cycle
$r\big{(}\frac{\chi(\varDelta t)}{2}\big{)}$ is larger than that at the
beginning $r\big{(}\frac{\chi(0)}{2}\big{)}$ or not 555Note that if $r(x)x$ is
monotonic, this is also a necessary condition. However, in general, this is
only a sufficient and not necessary condition. For example, there are cases
where both symmetric and asymmetric trajectories are stable if $r(x)x$ is non-
monotonic. For example, $r(x)x=\epsilon+\kappa
x(x^{2}-\frac{3}{2}(\alpha+\beta)+3\alpha\beta)$..

For example, if $r(x)x$ is linear, (e.g., $r(x)x={\epsilon+\kappa x}$, as is
the case for the competitive Azoarcus ribozymes discussed later in Sec.
II.6.1) only the growth state with $\delta=0$ (i.e., the symmetrical
trajectory) is always stable. Thus, for the system to show bistability and
heredity, $r(x)x$ must be a nonlinear function of $x$.

Figure 3: (A) Schematics of the competing autocatalytic entities, $\rm{X}_{1}$
and $\rm{X}_{2}$; they are converted from a substrate S (solid arrows)
catalyzed by itself (dashed arrows). (B) The red and blue curves represent the
nullclines for the composition of the entities $\rm{X}_{1}$ and $\rm{X}_{2}$
at just before the dilution, $x_{1}(-0)$ and $x_{2}(-0)$, respectively, in a
case with $r(x)x=\epsilon+\kappa x^{2}$. The cross points represent the
stable/unstable fixed points of concentrations of $\rm{X}_{1}$ and
$\rm{X}_{2}$ at just before the dilution, $x_{1}^{*}(-0)$ and $x_{2}^{*}(-0)$.
See Methods and Models for the detail of the drawing of the nullclines. We set
$\varDelta t=1$, $2$, $\kappa=8$, $\epsilon=0.5$, and $\bar{\phi}=1$.

###  II.3 Heredity of compositional state requires serial dilution interval to
be below a critical threshold

Next, we analyze how bistability varies across the parameter space for the
serial dilution (SD) protocol, consisting of the cycle interval $\varDelta t$
and the dilution factor $m$ ($=e^{\bar{\phi}\varDelta t})$. Here, we assume
nonlinear reproduction rate functions that can show the bistability (i.e.,
satisfying the condition Eq. 3 under the appropriate SD protocol). In this
subsection, we use $r(x)x=\epsilon+\kappa x^{2}$ as an example of such a
function. This function would be realized in a catalytic system consisting of
two chemical reactions that generate each autocatalytic entity: a
“spontaneous” or “background” chemical reaction where substrate S is converted
to species X at a rate $\epsilon$, and a catalyzed reaction where S is
converted to X in the presence of another molecule of X with the efficiency
$\kappa$ (The square dependence of $x$ is realized by dimeric catalyst [38],
or multi-step reactions such as those in the modified _Azoarcus_ system
described later in Sec. II.6.2.). However, the results in this subsection are
true for all nonlinear functions of $r(x)$ which are convex, i.e.,
$\frac{d^{2}r}{dx^{2}}>0$.

Fig. 4A shows the bifurcation diagram of the concentrations just before a
dilution in the stationary trajectory, as we vary the cycle interval
$\varDelta t$ in the SD protocol while keeping the dilution rate $\bar{\phi}$
fixed (note that the dilution factor $m$($=e^{\bar{\phi}\varDelta t}$) is
_not_ fixed). The bifurcation occurs at $\varDelta t=\varDelta t_{c}^{sd}$. If
$\varDelta t$ is more than this critical value, the system is no longer
bistable, i.e., it does not exhibit heredity of the compositional state.
Similarly, if we fix the interval $\varDelta t$ and vary the dilution rate
$\bar{\phi}$, the same bifurcation at which the bistability disappears is
observed at $\bar{\phi}={\phi}_{c}^{sd}$ (Fig. 4C). The phase diagram of the
parameters in the protocols, $\varDelta t$, and $\bar{\phi}$ is drawn in Fig.
4D.

Figure 4:  (A) The bifurcation diagram with a varying interval of each
dilution cycle $\varDelta t$ in a case with $r(x)x=\epsilon+\kappa x^{2}$. The
pitchfork bifurcation occurs at $\varDelta t=\varDelta t_{c}^{sd}$ when the
upper and lower fixed points merge into one fixed point. (B) The dotted line
represents the critical point $\varDelta t_{c}^{sd}$, which divides the
regions where the system has bistability or not. The solid line represents the
theoretical line for $\varDelta t_{c}^{sd}$ determined by the relation:
$\varDelta
t_{c}^{sd}\sim\frac{1}{\bar{\phi}}\log\left(\frac{\kappa}{4\epsilon}(s^{tot})^{2}\right).$
(C) The bifurcation diagram with a varying long-term dilution rate
$\bar{\phi}$. (D) The dotted line represents the critical point $\varDelta
t_{c}^{sd}$, which divides the regions where the system has bistability or
not. Each colored dotted line represents the difference in the background
reaction rate $\epsilon$. The theoretical lines (solid) are the same as in
(B). We set the default parameters as $\epsilon=0.5$, $\kappa=8$ and
$\bar{\phi}=1$.

The critical value $\varDelta t_{c}^{sd}$ depends on the reproduction rate
function $r(x)$ and its kinetic parameters (Fig. 4B). Using Eq. 3, the
critical $\varDelta t$ at which the system loses bistability, $\varDelta
t_{c}^{sd}$, in a case with $r(x)x=\epsilon+\kappa x^{2}$ is derived as
$\varDelta
t_{c}^{sd}\sim\frac{1}{\bar{\phi}}\log\left(\frac{\kappa}{4\epsilon}(s^{tot})^{2}\right)$
(Fig. 4B). Intuitively, we can realize this form as the condition that the
background reaction dominates the catalyzed reaction at the start of each
cycle (just after each dilution), i.e., $\epsilon>\kappa(x^{*})^{2}$, and
$x^{*}$ is roughly $x^{*}\sim\frac{1}{2}\frac{s^{tot}}{m}$ if all of the
substrate $\rm{S}$ added at the beginning of a cycle is converted to the
$\rm{X}_{1}$ or $\rm{X}_{2}$ by the end of the cycle.

The critical value for $\bar{\phi}^{sd}$ is determined in a similar way.
Especially, in the CSTR limit (i.e., $\varDelta t\rightarrow 0$), the
condition Eq. 3 becomes $\frac{dr}{dx}(\frac{\chi^{*}}{2})<0$, where
$\chi^{*}$ is such that $\chi^{*}r(\frac{\chi^{*}}{2})-\bar{\phi}=0$ (see
Appendix Sec A.2). Then, the critical dilution rate under CSTR,
$\phi_{c}^{cstr}$ is determined as $\bar{\phi}$ such that breaks the above
condition. In a case with $r(x)x=\epsilon+\kappa x^{2}$, $\phi_{c}^{cstr}$ is
derived as $\phi_{c}^{cstr}=2s^{tot}\sqrt{\epsilon\kappa}-4\epsilon$.

###  II.4 Critical interval in a general growth and division process is
bounded by that in the serial dilution protocol

Now, we return to the general GD protocols with the dilution rate
$\phi(t)={\frac{dV}{dt}}/{V}$ discussed in Sec. A. Under such GD protocols,
parametrized by cycle interval $\varDelta t$ and long-term dilution rate
$\bar{\phi}$, similar to the SD protocol, we discuss the condition for the
same ACSs to exhibit bistability.

Our main result is that the parameter region exhibiting bistability for the
general GD protocol is bounded by that of SD and CSTR: if $r(x)$ is a convex
function of $x$, i.e., $\frac{d^{2}r}{dx^{2}}>0$, all of the parameter regions
($\varDelta t$ and $\bar{\phi}$) where there is bistability under SD are
included within the bistable region under the general dilution protocols,
which in turn is included within the region exhibiting bistability under the
CSTR protocol (see Appendix Sec. A.3 for the proof).

For example, we can consider the compartment growth dynamics obeying,
$\frac{dV}{dt}=\bar{\phi}_{\alpha}V^{\alpha},$ where $\alpha$ is the order of
the growth, and $\bar{\phi}_{\alpha}$ is a constant depending on $\alpha$.
Here, to compare sensibly across the different protocols, we assume the growth
rate of the compartment volume on long timescales in each case is the same,
i.e., SD with the dilution factor $m=e^{\bar{\phi}\varDelta t}$ or the CSTR
with the dilution rate $\bar{\phi}$ (see Methods and Models for the details).
As Fig. 5 shows, the critical interval $\varDelta t_{c}$ that is the upper
limit for a system with bistability in the general case is bounded from below
by the critical $\varDelta t^{sd}_{c}$ for the SD protocol: $\varDelta
t_{c}^{sd}\leq\varDelta t_{c}<\varDelta t_{c}^{cstr},$ under the fixed
$\bar{\phi}$, where $\varDelta t_{c}^{cstr}$ is infinite or otherwise zero
($\varDelta t_{c}^{cstr}=0$ means that there is no bistability under any
$\varDelta t$ values). On the other hand, if $\varDelta t$ is fixed,
$\phi_{c}^{sd}\leq\phi_{c}\leq\phi_{c}^{cstr},$ where $\phi_{c}^{sd}$,
$\phi_{c}^{cstr}$ and $\phi_{c}$ are the threshold for $\phi$ under SD, CSTR
and general protocols, respectively.

Figure 5: Phase diagram for the system with/without heredity using the
dilution rate $\bar{\phi}$ and dilution interval $\varDelta t$ as parameters,
in a case with $r(x)x=\epsilon+\kappa x^{2}$. The solid line represents the
boundary between with and without bistability under SD. The dashed lines
represent that under the GD with different growth laws (the exponents $\alpha$
of the growth law for compartment; $\alpha=2,4,6$). The vertical thin line
represents that under GD with $\alpha=1$, corresponding to the CSTR protocol.
We set the parameters as $\kappa=8,\epsilon=0.5$.

###  II.5 Robustness of results

Two key results we have are: (i) for a serial dilution protocol, there is a
critical dilution interval and dilution rate for the inheritance to be
possible, (ii) these critical values provide bounds for the properties of a
general growth and division process in order to enable inheritance. We have
checked the robustness of these results to a number of variations and
extensions of our models, namely:

1\. Autocatalytic sets with reversible chemical reactions (Appendix Sec. B.1)
Making all reactions reversible is more chemically realistic and allows
convergence to thermal equilibrium in the absence of dilution protocols.
Interestingly, we find that, unlike the irreversible reaction case, the region
of bistability is bounded in the reversible case for the parameter
$\bar{\phi}$. That is, for the SD protocol with fixed $\varDelta t$, there is
both an upper and a lower critical $\bar{\phi}$ (see Fig. S3 in Appendix).
Importantly, even with reversible reactions, the parameter space for the
general GD protocol contains that for SD protocols, as we found for
irreversible reactions.

2\. ACSs with asymmetric kinetic rates (Appendix Sec. B.2)
Our results are similar even when the two competing ACSs have different
kinetic rate constants (catalytic efficiency), although the bifurcation where
the bistability disappears is discontinuous (see Fig. S4). The boundary in the
parameter space for ACSs that have bistability under general GD protocols is
also bounded by the boundaries for the SD and CSTR protocols (see Appendix).

3\. Reproduction rates that depend on both entities (Appendix Sec. B.3)
The rate equation Eq. 1 and the condition Eq. 3 can be extended to a more
general reproduction rate function $r(x_{1},x_{2})$, which also depends on the
concentration of the other entity; for example, this includes mutual
inhibition between the entities. Our results are unchanged by this
generalization.

4\. Stochasticity in reaction systems or protocols (Appendix Sec. B.4)
If the number of molecules is small (e.g., the reaction dynamics occur inside
sufficiently small compartments), stochastic fluctuations are not negligible.
Close to the deterministic bifurcation transitions, these fluctuations cause
random transitions between states, which destroys the information to be
inherited. Thus, the parameter space where the system has heredity is narrower
than the deterministic case. However, the existence of a critical threshold in
the dilution interval or the dilution rate and the bistable parameter region
for general GD cycles being bounded by that of SD and CSTR, are robust to the
addition of stochastic noise in the chemical reaction system. The results are
also robust to stochasticity in the dilution protocols; for example, if the
division time $\varDelta t$ is randomly distributed (see Appendix Sec. B.4.2).

5\. Differential reproduction of different chemical compositional states
(Appendix Sec. B.5)
Finally, we investigated situations where the dilution rates depend on the
chemical composition, and the system has different growth rates for the two
growth states. As noted in the introduction, this is, in fact, a crucial
property necessary for a population of compartmentalized chemical reaction
systems to undergo Darwinian evolution. We consider cases where $\phi$ depends
both symmetrically or asymmetrically on $x_{1}$ and $x_{2}$. We have also
checked the case where the division interval, not the rate, depends on the
chemical composition. In all cases, our results continue to hold.

###  II.6 Autocatalytic sets based on the _Azoarcus_ ribozyme

In this section, we apply the framework we have developed to an experimentally
realized ACS based on the _Azoarcus_ ribozyme [47, 48, 30]. We examine a few
simplifications and variants of this ribozyme system, along with systems
consisting of two _Azoarcus_ ribozymes competing for the same nutrients. We
show that some of these variants can exhibit two (exponential) growth states,
and some are not. In particular, we find a modified version of the _Azoarcus_
ribozyme, which incorporates additional catabolic steps [51], can exhibit
bistability under competition for shared resources. Applying our results on
critical thresholds to this system, we can suggest a serial dilution protocol
that can test whether this modified _Azoarcus_ system can exhibit the
inheritance of its phenotypic state.

The Azoarcus ribozyme $\mathbf{WXYZ}$ can be assembled from two fragments by
the reaction,

${}_{\mathrm{{M}}}\mathbf{WXY}_{\mathrm{{N}}}+\mathbf{Z}\rightarrow\
{}_{\mathrm{{M}}}\mathbf{WXY}_{\mathrm{{N}}}\mathbf{Z}.$ (4)

${}_{\mathrm{{M}}}\mathbf{WXY}_{\mathrm{{N}}}$ and $\mathbf{Z}$ are the
fragments, and ${\mathrm{{M}}}$ and ${\mathrm{{N}}}$ are bases at the ends of
$\mathbf{WXY}$ 666The bases M and N are the middle nucleotide of the 3nt
recognition element in WXY, called the internal guide sequence (IGS) and the
tag sequence [48].. $\mathrm{{M}}$ and $\mathrm{{N}}$ can be arbitrary bases
$\mathrm{\\{A,U,C,G\\}}$, thus there are 16 different types of this engineered
Azoarcus ribozyme. This reaction is catalyzed specifically by a ribozyme if
$\mathrm{{M}}$ in the fragment is a complementary base to $\mathrm{{N}}$ in
the ribozyme. Certain types of the ribozyme can catalyze the formation
reaction of themselves: e.g.,
${}_{\mathrm{{M}}}\mathbf{WXY}_{\mathrm{{N}}}\mathbf{Z}$ such that
$\mathrm{{M}}=\mathrm{C}$ and $\mathrm{{N}}=\mathrm{G}$ or
$\mathrm{{M}}=\mathrm{U}$ and $\mathrm{{N}}=\mathrm{A}$. Besides being
catalyzed by the corresponding ribozymes, the reaction is also (weakly)
catalyzed by a non-covalent complex between the corresponding fragments, or
non-specifically by non-corresponding ribozymes; we call these ‘background
reactions’ [48].

####  II.6.1 Absence of bistability in competing ACSs based on the original
engineered Azoarcus ribozyme

We imagine a particular case of the system in Fig. 6A, where $\rm{X_{1}}$ and
$\rm{X_{2}}$ are two distinct types of self-catalyzing Azoarcus ribozyme, made
from two distinct $\mathbf{WXY}$ fragments and a common $\mathbf{Z}$ fragment.
For simplicity, we assume the two sets of reactions occur with symmetric
kinetic rates. We assume the $\mathbf{WXY}$ fragments are abundant, whereas
$\mathbf{Z}$ is not and limits the reaction rates. Thus, $\mathbf{Z}$ acts as
the common substrate S.

Then, the rate equations for the concentrations of $\rm{X_{1}}$ and
$\rm{X_{2}}$, $x_{1}$ and $x_{2}$ are described as Eq. 1 with the linear
reproduction function $r(x)x=\epsilon+\kappa x$, where $\epsilon$ is the rate
constant of the background reaction, and $\kappa$ is the catalytic strength of
$\rm{X}_{1}$ and $\rm{X}_{2}$. 777Here, the background reaction rate is
approximated as a constant, although this reaction is due to catalyzed
reaction by the non-covalent ribozymes or non-corresponding ribozymes.
Assuming that the WXY fragments are abundant, the concentration of non-
covalent ribozymes WXY:Z is $\sim s$ approximately. We assume their catalytic
activity is non-specific, whose reaction rate is $\epsilon$. Next, we assume
corresponding and non-corresponding ribozymes catalyze at efficiency
$\tilde{\kappa}^{\prime}$ and $\tilde{\kappa}^{\prime}$, respectively. Then,
the rate of the reaction $\rm{S}\rightarrow\rm{X_{1}}$ is
$s\epsilon+x_{2}\tilde{\kappa}^{\prime}+x_{1}\tilde{\kappa}$. If we assume
that $\epsilon=\tilde{\kappa}^{\prime}$, the rate is given as $\epsilon+\kappa
x_{1}$, where $\kappa=\tilde{\kappa}_{1}-\epsilon$. As mentioned earlier, for
simplicity, we assume the catalytic efficiency of $\rm{X_{1}}$ and
$\rm{X_{2}}$ are equal.

As already discussed in Sec. B, competing ACSs with such linear reproduction
rate functions cannot be bistable – for all initial conditions, the system
eventually reaches the state with an equal amount of the two ribozymes. More
precisely, in Sec. B, only the local stability of the symmetric state is
shown. In this case, however, we can further show the _global_ stability of
this symmetric state 888Here, the time derivative of difference of
concentration between $\rm{X_{1}}$ and $\rm{X_{2}}$ is
$\frac{d}{dt}\left(\frac{\delta}{\chi}\right)=\frac{1}{\chi^{2}}(\frac{d\delta}{dt}\chi-\delta\frac{d\chi}{dt})=-2\epsilon
s\frac{\delta}{\chi},$ where $\delta=x_{1}-x_{2}$ and $\chi=x_{1}+x_{2}$.
Therefore, $\delta/\chi$ decays exponentially into zero in the characteristic
relaxation time scale $\tau=\frac{1}{2s\epsilon}$. For example, if
$s\epsilon\sim 0.01\min^{-1}$, the half time $\tau$ is estimated as $\sim
50\min$.. Therefore, the system has no heredity, which is consistent with
previous experiments [30].

####  II.6.2 Azoarcus system coupled with metabolism exhibits bistability

As discussed in Sec. B, a reproduction rate function $r(x)$ with a higher
order of catalysis is necessary for bistability. In the Azoarcus system, this
has previously been realized by engineering a variant where the system is
coupled to catabolism and anabolism reactions:

$\begin{split}{}_{\mathrm{{M}}}\mathbf{WXY}_{\mathrm{{N}}}\mathbf{\mathchar
45\relax
mod}\rightarrow{}_{\mathrm{{M}}}\mathbf{WXY}_{\mathrm{{N}}}+\mathbf{\mathchar
45\relax mod},\\\
{}_{\mathrm{{M}}}\mathbf{WXY}_{\mathrm{{N}}}+\mathbf{Z}\rightarrow\
{}_{\mathrm{{M}}}\mathbf{WXY}_{\mathrm{{N}}}\mathbf{Z},\end{split}$ (5)

where $\mathbf{\mathchar 45\relax mod}$ represents an extra sequence joined to
fragments $\mathbf{WXY}$. Here, the first reaction represents a catabolism
reaction that processes the modified fragment to a substrate that can
participate in ribozyme synthesis. The second one represents an anabolic
reaction that joins the fragments to form the ribozyme. We denote $\mathbf{Z}$
as $\rm{S}$, $\mathbf{WXYZ}$ as $\rm{X}_{1}$ and $\rm{X}_{2}$ and
$\mathbf{WXY}$ as $\rm{X^{\prime}_{1}}$ and $\rm{X^{\prime}_{2}}$. Two
ribozymes $\rm{X_{1}}$ and $\rm{X_{2}}$ synthesize themselves from the shared
substrate $\rm{S}$, and there is the intermediate state $\rm{X^{\prime}_{1}}$
and $\rm{X^{\prime}_{2}}$ during the synthesis (see Fig. 6B).

Figure 6: Schematic diagram of the model chemical reaction systems: (A) non-
modified two competing autocatalytic ribozymes and (B) the ones coupled with
anabolism and catabolism reactions. Solid arrows represent reactions and
dashed represent catalysis. Species enclosed by dashed circles are supplied
from outside. (C) Time course of the difference of the concentration of the
ribozyme species, $\delta=x_{1}-x_{2}$, starting from the initial condition
where $\delta=1$ or $=-1$. The gray circle represents the time course for the
model (A) under the SD protocol with the interval $\varDelta t=1.5$ (Eq. 4
with $r(x)x=\epsilon+\kappa x$; we considered the reversible reaction). The
black circle and square represent time courses for the model (B) under the SD
with $\varDelta t=1.5$ and $=3.0$, respectively. We set
$\epsilon=0.25,\kappa=2.5$, and $b=0.1$. (D) The parameter region for the
system with the bistability. In the dark-grey area, the system has bistability
under SD protocol. The thin vertical lines (the boundaries for the light-grey
area) represent the lower and upper boundary for $\bar{\phi}$ for the
bistability under the CSTR. We set $\epsilon=0.25,\kappa=2.5$, and $b=0.1$.

The concentrations of chemical species obey the rate equations

$\displaystyle\frac{dx^{\prime}_{i}}{dt}$ $\displaystyle=(\epsilon+\kappa
x_{i})(1-(1+b)sx^{\prime}_{i}+bx_{i}),$ (6) $\displaystyle\frac{dx_{i}}{dt}$
$\displaystyle=(\epsilon+\kappa x_{i})(sx^{\prime}_{i}-bx_{i}),$

where $i=1$ or $2$, $\epsilon$ is the spontaneous reaction rate, $\kappa$ is
the catalytic efficiency of the ribozymes, and $b$ ($\ll 1$) is the relative
rate of backward reaction compared with the forward one 999The ACSs based on
the _Azoarcus_ system show slow backward reactions since they are based on the
recombination reactions of nucleotides [51].. We again consider the ACSs
system under the SD protocol with the interval $\varDelta t$ and an $m$-fold
dilution factor. Also, the total concentration of S, ${\rm X_{1}}$ and ${\rm
X_{2}}$ is kept as a constant $s^{tot}=s+x_{1}+x_{2}$.

In contrast with the original system (Fig. 6A), this system indeed exhibits
bistability for appropriate values of the parameters $\varDelta t$ and
$\bar{\phi}$ (Fig. S12), i.e., it reaches either the $\rm{X}_{1}$-dominated
state or the $\rm{X}_{2}$-dominated state, depending on the initial
composition of the species (Fig. 6C).

####  II.6.3 Bounds on the critical $\varDelta t$ and $\bar{\phi}$ for
observing inheritance of compositional state

As in the previous models, there is the region for the kinetic parameters,
$\varDelta t$ and $\bar{\phi}$, where the system exhibits bistability (see
Fig. 6D). For this modified Azoarcus system, the concrete value that would
exhibit heredity under SD is predicted that the dilution interval lies within
50-125 min, and the dilution factor per cycle lies between 2.5-11-fold (We
assumed $\kappa=0.1\rm{\mu M}^{-2}\min^{-1}$, $\epsilon=0.01\rm{\mu
M}^{-1}\min^{-1}$, $\bar{\phi}=0.04\min^{-1}$ and $b=0.1$; see Fig. S13)

Further, we checked numerically (see Fig. 6D and Appendix Sec. C) that under
the alternative dilution protocols (various functions of $\phi(t)$) upper
bound of $\varDelta t$ for the bistability is bounded by $\varDelta
t_{c}^{sd}$. These results for the _Azoarcus_ based system are robust even if
the reaction network is modified: we have tested alternative structures,
asymmetric kinetic rates, and stochastic reactions or compartment dynamics
(see Appendix Sec. D).

##  III Discussion

We have studied mathematical models of a very general class of chemical
reaction systems consisting of two ACSs competing for a shared resource. When
enclosed within growing and dividing compartments, such a system is one of the
simplest capable of exhibiting the heredity of its compositional state. In
other words, it can exhibit two distinct compositional states that are stable
to the growth and division of the compartment. This, along with differential
reproduction rates and some sources of variation of the compositional state,
is a key property for a chemical system to form a Darwinian population.

Our main result is that such bistable chemical systems exhibit the heredity of
their two compositional states when the growth and division cycle parameters
lie within certain critical bounds. These are bounded by the critical values
for these parameters in a serial dilution protocol. Thus, a serial dilution
experiment, which is much easier to set up in a laboratory than a general
growth and division scenario, can be used to test whether a given
autocatalytic chemical system can exhibit heredity of its compositional
states.

We used these results to analyze a chemical system based on the _Azoarcus_
ribozyme, which has been realized experimentally, and predicted that it would
indeed exhibit heredity under serial dilution for a specific range of values
for the dilution interval and dilution factor (see Results section II.6.3).
The _Azoarcus_ system is a particularly intriguing experimental system for
studying heredity because it can go beyond our present theoretical analysis in
two ways:
(i) It can be engineered to exhibit any number of compositional states. For
example, by choosing bases M and N to be AU, UA, CG, and GC, one can have four
ACSs competing for the same resource. Such a system can encode more than 1-bit
of information 101010Moreover, the _Azoarcus_ ribozyme can form cross-
catalytic networks [47, 30] (e.g., choosing bases M and N to be CC and GG, or
AA and UU); in such cases, each network is a unit of self-reproduction (called
an ‘autocatalytic core’) and could compete with other units..
(ii) Each compositional state can be composed of a large number of chemical
species. The _Azoarcus_ ribozyme can catalyze not only the formation of itself
but also the production of diverse RNA sequences [57]. Here, self-reproducing
ribozymes (ACSs) correspond to the ‘autocatalytic cores,’ and the sequences
produced by the ribozyme correspond to their ‘peripheries’ [58, 59, 25].
Species in the periphery can nevertheless play important roles, for instance,
in the differential reproduction of the compositional states.
Of course, extending our present theoretical analysis to multistable systems
which can encode more than 1-bit of information and to more complex
autocatalytic systems consisting of cores and peripheries is relatively
straightforward, and easy to be implemented experimentally. The _Azoarcus_
system can guide such theoretical extensions, but we expect our core results
to remain the same for such more complex reaction systems.

The condition for bi/multistability of autocatalytic systems has been
previously discussed in the context of self-reproduction and the origins of
life. For example, Giri-Jain [40] found a class of ACSs that exhibit
bistability under CSTR conditions, with one state growing and the other a non-
growing state [41]. Remarkably, their network model also needs two steps to
show bistability (if mass-action kinetics is assumed), like the _Azoarcus_
system coupled to metabolic reactions. Note that the bistability in their
model needs high catalytic efficiency of self-catalysts (e.g., $10^{4}$ order
of magnitude) 111111This type of bistability: low and high catalyst
concentration states, also could appear in our model if the background
reaction rate is small enough (nearly zero). However, such bistability is
easily destroyed under the SD protocol unless the interval is very short.. It
is noteworthy that our model exhibits bistability even with relatively low
efficiency of catalysts (or equivalently, with a high rate of background
reactions), which is more plausible for a prebiotic scenario.

Some previous studies have shown how heredity of composition can arise even
without bistability due to the difference in reproduction rate and competition
between compartments. For example, see Segre-Lancet [23], Kaneko-Yomo [61], or
Doulcier et al. [62], for the models with ‘multilevel selection’ [63]
(including ‘transient compartmentalization’ [64]). However, such mechanisms
require fine-tuning for protocols, such as the size of compartments, and it is
controversial whether the compositional states can be sustained enough for
selection pressures to act on them [24]. In contrast, we demonstrated that
compositional information in our models is robustly inherited in a single
lineage of compartments across a wide range of values for kinetic constants or
protocol parameters. Note that compartments such as coacervates are known for
their high permeability, and they exchange their contents (including large
ribozymes) at faster rates. Even in such a case, the inheritance of
compositional information could be robust for the kind of bistable systems we
have studied. [65].

##  IV Methods and Models

###  IV.1 Simulation details

All simulations of chemical reaction dynamics involved solving ordinary
differential equations (and a discrete mapping). This was done using the
adaptive Runge-Kutta method (Dormand–Prince method) coded by the C++ language
[66].

_The serial dilution protocol:_ The composition of entities
$x_{1}(t),x_{2}(t),...,x_{M}(t)$ under the serial dilution (SD) protocol
evolve according to a rate equation,

$\frac{d\bm{x}}{dt}=\bm{f}(\bm{x}),$ (7)

where $\bm{x}(t)$ is a time-dependent compositional vector (e.g., in a case
with a model in Sec. II.F,
$\bm{x}(t)=\\{s(t),x^{\prime}_{1}(t),x^{\prime}_{2}(t),x_{1}(t),x_{2}(t)\\}$).
At each interval $t=n\varDelta t$ ($n=1,2,\dots$), the all entities are
diluted (and the substrate S is added), that is, the composition change by a
discrete mapping:

$\bm{x}(n\varDelta t+0)=\frac{1}{m}\bm{x}(n\varDelta
t-0)+s^{tot}(1-\frac{1}{m})\hat{\bm{1}}_{s},$ (8)

where $\hat{\bm{1}}_{s}$ is a unit vector for the substrate S (e.g., in a case
with a model in Sec. II.F, $\hat{\bm{1}}_{s}=\\{1,0,0,0,0\\}$), $n=0,1,\dots$,
and $n\varDelta t-0$ and $n\varDelta t+0$ represent the time right before and
after the dilution at $t=n\varDelta t$ ($n=0,1,\dots$). If the reaction
dynamics in Eq. 14 does not change the total sum of composition,
$x^{tot}=\sum_{i}x_{i}$, the repeats of the dilution cycle, i.e.,the mapping
in Eq. 15 results in the steady state with $x^{tot}=s^{tot}$. In the main
text, we fix $x^{tot}$ as $s^{tot}$, and the concentration of free substrate
is $s=s^{tot}-\sum x_{i}$ in the rate equations.

If we set $\varDelta t$ enough small, i.e., the system is diluted repeatedly
at quite a short interval, and set $m$ and $s^{tot}$ as
$m=e^{\bar{\phi}\varDelta t}$ and $s^{tot}=\frac{\bar{\sigma}}{\bar{\phi}}$,
the dynamics and the steady-state of the species are the same as that in the
continuous stirred-tank reactor (CSTR) with a dilution rate $\bar{\phi}$ and a
substrate supply rate $\bar{\sigma}$ [45].

_The dilution by the growth of the compartment:_ Formally, the rate equation
for the chemical composition $\bm{x}$ under the general dilution protocols is
expressed as,

$\frac{d\bm{x}}{dt}=\sigma(t)\hat{\bm{1}}_{s}+\bm{f}(\bm{x})-\phi(t)\bm{x},$
(9)

where $\bm{f}(\bm{x})$ is an arbitrary reaction dynamics, and $\sigma(t)$ and
$\phi(t)$ is the time-dependent supply rate of S and the dilution rate. If
$\sigma(t)=\phi(t)$ and the dynamics $\bm{f}(\bm{x})$ conserves the total
concentration $x^{tot}$, then $x^{tot}$ is constant at the steady state.

In the case of the compartment growth scenario, $\phi(t)$ is determined as
$\phi(t)={\frac{dV}{dt}}/{V}$, where $V$ is the volume of the compartment. For
example, we can consider the arbitrary order of compartment growth.

$\frac{dV}{dt}=\bar{\phi}_{\alpha}V^{\alpha},$ (10)

where $\alpha$ is the order of the growth, and $\bar{\phi}_{\alpha}$ is a
constant depending on $\alpha$. If $\alpha=1$, the growth is exponential. For
example, if we consider the volume growth is proportional to the surface area,
i.e., $\frac{dV}{dt}=r_{\alpha}S$, and if the vesicle is spherical, the
surface area $S$ is $S=V^{\frac{2}{3}}$, then $\alpha$ = 2/3 [67, 68]. The
volume $V$ is solved as
$V(t)=\left(\bar{\phi}_{\alpha}(1-\alpha)t+V_{0}^{1-\alpha}\right)^{\frac{1}{1-\alpha}}$.
Here, we assume the growth speed of compartment’s volume in the long time
scale is the same as the exponential growth with rate $\bar{\phi}$ (i.e., CSTR
with dilution rate $\bar{\phi}$); that is, $V(\varDelta
t)=V_{0}e^{\bar{\phi}\varDelta t}$. Then, $\bar{\phi}_{\alpha}$ should be
$\bar{\phi}_{\alpha}=\frac{V_{0}^{1-\alpha}}{(1-\alpha)\varDelta
t}(e^{\bar{\phi}(1-\alpha)\varDelta t}-1)$. Therefore, the dilution rate is

$\phi(t)=\frac{e^{\bar{\phi}(1-\alpha)\varDelta t}-1}{(1-\alpha)\varDelta
t\left(1+\frac{t}{\varDelta t}(e^{\bar{\phi}(1-\alpha)\varDelta
t}-1)\right)}.$ (11)

Note that $\phi(t)$ approaches $\phi(t)=\bar{\phi}$ as $\alpha\rightarrow 1$
(i.e., the same as the condition under the CSTR). On the other hand, $\phi(t)$
approaches $\phi(t)=\infty$ if $t=n\varDelta t\quad(n=1,2,\dots)$ or
$\phi(t)=0$ if otherwise as $\alpha\rightarrow\infty$ (i.e., the serial
dilution condition).

_Nullclines under SD protocol:_ The map between the chemical composition at
the beginning (end) of one cycle to that at the beginning (end) of the next
cycle, $P:\bm{x}(n\varDelta t+0)\mapsto\bm{x}\big{(}(n+1)\varDelta
t+0\big{)},n=0,1,..$, is interpreted as the Poincaré map obtained using the
Poincaré section: $t=n\varDelta t$. The trajectory $\bm{x}(t)$ is stable if
and only if the corresponding fixed point in the Poincaré map is stable. Then,
the stability of the trajectories can be determined from the intersections of
the nullclines of the discrete map (see Fig. 3B for a depiction of these
nullclines). The nullcline for $\rm X_{1}$, $x_{1}=f(x_{2})$, is obtained by
the fixing the concentration of $\rm X_{2}$ at the beginning of every cycle,
$x_{2}(n\varDelta t+0)=x_{2}$, and calculating the stationary concentration of
$\rm X_{1}$ at the beginning of cycles $x_{1}=x_{1}^{*}(+0)$ by repeating the
map enough times.

###### Acknowledgements.

We thank Sanjay Jain and Angad Yuvraj for the fruitful discussions. We
acknowledge support from the Indo-French Centre for the Promotion of Advanced
Research under project no. 5904-3, the Department of Atomic Energy (India)
under project no. RTI4006, the Simons Foundation (Grant No. 287975), and
computational facilities at NCBS.

## References

  * Woese [1967] C. R. Woese, The genetic code :the molecular basis for genetic expression., Proc. Natl. Acad. Sci. U.S.A (1967).
  * Crick [1968] F. H. Crick, The origin of the genetic code, J. Mol. Biol. 38, 367 (1968).
  * Cairns-Smith [1966] A. G. Cairns-Smith, The origin of life and the nature of the primitive gene, Journal of Theoretical Biology 10, 53 (1966).
  * Gilbert [1986] W. Gilbert, Origin of life: The rna world, nature 319, 618 (1986).
  * Kauffman [1986] S. A. Kauffman, Autocatalytic sets of proteins, Journal of Theoretical Biology 119, 1 (1986).
  * Jain and Krishna [1998] S. Jain and S. Krishna, Autocatalytic sets and the growth of complexity in an evolutionary model, Physical Review Letters 81, 5684 (1998).
  * Blokhuis  _et al._ [2020] A. Blokhuis, D. Lacoste, and P. Nghe, Universal motifs and the diversity of autocatalytic systems, Proceedings of the National Academy of Sciences 117, 25230 (2020).
  * Ameta  _et al._ [2021a] S. Ameta, Y. J. Matsubara, N. Chakraborty, S. Krishna, and S. Thutupalli, Self-reproduction and darwinian evolution in autocatalytic chemical reaction systems, Life 11, 308 (2021a).
  * Adamski  _et al._ [2020] P. Adamski, M. Eleveld, A. Sood, Á. Kun, A. Szilágyi, T. Czárán, E. Szathmáry, and S. Otto, From self-replication to replicator systems en route to de novo life, Nature Reviews Chemistry 4, 386 (2020).
  * Mizuuchi and Ichihashi [2021] R. Mizuuchi and N. Ichihashi, Primitive compartmentalization for the sustainable replication of genetic molecules, Life 11, 191 (2021).
  * Chen  _et al._ [2004] I. A. Chen, R. W. Roberts, and J. W. Szostak, The emergence of competition between model protocells, Science 305, 1474 (2004).
  * Zwicker  _et al._ [2017] D. Zwicker, R. Seyboldt, C. A. Weber, A. A. Hyman, and F. Jülicher, Growth and division of active droplets provides a model for protocells, Nature Physics 13, 408 (2017).
  * Boerlijst and Hogeweg [1991] M. C. Boerlijst and P. Hogeweg, Spiral wave structure in pre-biotic evolution: hypercycles stable against parasites, Physica D: Nonlinear Phenomena 48, 17 (1991).
  * Szabó  _et al._ [2002] P. Szabó, I. Scheuring, T. Czárán, and E. Szathmáry, In silico simulations reveal that replicators with limited dispersal evolve towards higher efficiency and fidelity, Nature 420, 340 (2002).
  * Krieger  _et al._ [2020] M. S. Krieger, S. Sinai, and M. A. Nowak, Turbulent coherent structures and early life below the kolmogorov scale, Nature communications 11, 1 (2020).
  * Charlat  _et al._ [2021] S. Charlat, A. Ariew, P. Bourrat, M. Ferreira Ruiz, T. Heams, P. Huneman, S. Krishna, M. Lachmann, N. Lartillot, L. Le Sergeant d’Hendecourt, _et al._ , Natural selection beyond life? a workshop report, Life 11, 1051 (2021).
  * Godfrey-Smith [2007] P. Godfrey-Smith, Conditions for evolution by natural selection, The Journal of Philosophy 104, 489 (2007).
  * Lewontin [1970] R. C. Lewontin, The units of selection, Annual Review of Ecology and Systematics , 1 (1970).
  * Endler [2020] J. A. Endler,  _Natural Selection in the Wild.(MPB-21), Volume 21_ (Princeton University Press, 2020).
  * Ridley [2003] M. Ridley,  _Evolution, 3rd Edition_ (Wiley-Blackwell, 2003).
  * Note [1] This could occur via a variety of mechanisms. For example, osmotic pressure due to the difference in the composition between the inside and outside of the compartment may induce its growth [11]. The ACS could also produce the precursors of the compartment (e.g., lipid molecules) [69, 70]. Another possibility is that the ACS energetically drives the growth, and shape instability triggers the division of compartments [12].
  * Note [2] In dynamical systems, stability often refers to stability against random stochasticity, e.g., due to thermal noise. It is true that for an autocatalytic chemical system to exhibit an inheritance of its phenotype (its chemical composition), it must also exhibit a certain amount of stability against noise. Some small enough probability of transitions to different states due to noise can be subsumed under phenotypic variation (indeed, this may be the only source of variation available), but too much will destroy the property of the heredity of states. Later we provide some results from stochastic simulations, but we largely assume that if the system, in the absence of noise, is stable upon division, then it satisfies the third property of inheritance of phenotypic traits.
  * Segré  _et al._ [2000] D. Segré, D. Ben-Eli, and D. Lancet, Compositional genomes: prebiotic information transfer in mutually catalytic noncovalent assemblies, Proceedings of the National Academy of Sciences 97, 4112 (2000).
  * Vasas  _et al._ [2010] V. Vasas, E. Szathmáry, and M. Santos, Lack of evolvability in self-sustaining autocatalytic networks constraints metabolism-first scenarios for the origin of life, Proceedings of the National Academy of Sciences 107, 1470 (2010).
  * Vasas  _et al._ [2012] V. Vasas, C. Fernando, M. Santos, S. Kauffman, and E. Szathmáry, Evolution before genes, Biology direct 7, 1 (2012).
  * Hordijk  _et al._ [2012] W. Hordijk, M. Steel, and S. Kauffman, The structure of autocatalytic sets: Evolvability, enablement, and emergence, Acta Biotheoretica 60, 379 (2012).
  * Guttenberg  _et al._ [2015] N. Guttenberg, M. Laneuville, M. Ilardo, and N. Aubert-Kato, Transferable measurements of heredity in models of the origins of life, PloS one 10, e0140663 (2015).
  * Lancet  _et al._ [2018] D. Lancet, R. Zidovetzki, and O. Markovitch, Systems protobiology: origin of life in lipid catalytic networks, Journal of The Royal Society Interface 15, 20180159 (2018).
  * Colomer  _et al._ [2020] I. Colomer, A. Borissov, and S. P. Fletcher, Selection from a pool of self-assembling lipid replicators, Nature Communications 11, 1 (2020).
  * Ameta  _et al._ [2021b] S. Ameta, S. Arsène, S. Foulon, B. Saudemont, B. E. Clifton, A. D. Griffiths, and P. Nghe, Darwinian properties and their trade-offs in autocatalytic rna reaction networks, Nature Communications 12, 1 (2021b).
  * Peng  _et al._ [2022] Z. Peng, J. Linderoth, and D. A. Baum, The hierarchical organization of autocatalytic reaction networks and its relevance to the origin of life, PLOS Computational Biology 18, e1010498 (2022).
  * Gardner  _et al._ [2000] T. S. Gardner, C. R. Cantor, and J. J. Collins, Construction of a genetic toggle switch in escherichia coli, Nature 403, 339 (2000).
  * Ferrell Jr and Xiong [2001] J. E. Ferrell Jr and W. Xiong, Bistability in cell signaling: How to make continuous processes discontinuous, and reversible processes irreversible, Chaos: An Interdisciplinary Journal of Nonlinear Science 11, 227 (2001).
  * Peng  _et al._ [2020] Z. Peng, A. M. Plum, P. Gagrani, and D. A. Baum, An ecological framework for the analysis of prebiotic chemical reaction networks, Journal of Theoretical Biology 507, 110451 (2020).
  * Frank [1953] F. C. Frank, On spontaneous asymmetric synthesis, Biochimica et biophysica acta 11, 459 (1953).
  * Laurent  _et al._ [2021] G. Laurent, D. Lacoste, and P. Gaspard, Emergence of homochirality in large molecular systems, Proceedings of the National Academy of Sciences 118, e2012741118 (2021).
  * Maity  _et al._ [2019] I. Maity, N. Wagner, R. Mukherjee, D. Dev, E. Peacock-Lopez, R. Cohen-Luria, and G. Ashkenasy, A chemically fueled non-enzymatic bistable network, Nature Communications 10, 1 (2019).
  * Wagner  _et al._ [2020] N. Wagner, R. Mukherjee, I. Maity, S. Kraun, and G. Ashkenasy, Programming multistationarity in chemical replication networks, ChemSystemsChem 2, e1900048 (2020).
  * Schlögl [1972] F. Schlögl, Chemical reaction models for non-equilibrium phase transitions, Zeitschrift für physik 253, 147 (1972).
  * Giri and Jain [2012] V. Giri and S. Jain, The origin of large molecules in primordial autocatalytic reaction networks, PloS ONE 7, e29546 (2012).
  * Matsubara and Kaneko [2016] Y. J. Matsubara and K. Kaneko, Optimal size for emergence of self-replicating polymer system, Physical Review E 93, 032503 (2016).
  * Wilhelm [2009] T. Wilhelm, The smallest chemical reaction system with bistability, BMC systems biology 3, 1 (2009).
  * Craciun  _et al._ [2006] G. Craciun, Y. Tang, and M. Feinberg, Understanding bistability in complex enzyme-driven reaction networks, Proceedings of the National Academy of Sciences 103, 8697 (2006).
  * Note [3] The necessary conditions for bistable chemical reactions are: (i) positive feedback (e.g., autocatalysis), (ii) filtering noise, and (iii) preventing explosion (e.g., conservation law of the components) [43]. In addition, ‘nonlinearity (or ‘ultrasensitivity’) is required in the positive feedback [33]. Later, we discuss the minimum autocatalytic chemical system that satisfies the above condition.
  * Blokhuis  _et al._ [2018] A. Blokhuis, D. Lacoste, and P. Gaspard, Reaction kinetics in open reactors and serial transfers between closed reactors, The Journal of Chemical Physics 148, 144902 (2018).
  * Erez  _et al._ [2020] A. Erez, J. G. Lopez, B. G. Weiner, Y. Meir, and N. S. Wingreen, Nutrient levels and trade-offs control diversity in a serial dilution ecosystem, Elife 9, e57790 (2020).
  * Vaidya  _et al._ [2012] N. Vaidya, M. L. Manapat, I. A. Chen, R. Xulvi-Brunet, E. J. Hayden, and N. Lehman, Spontaneous network formation among cooperative rna replicators, Nature 491, 72 (2012).
  * Yeates  _et al._ [2016] J. A. Yeates, C. Hilbe, M. Zwick, M. A. Nowak, and N. Lehman, Dynamics of prebiotic rna reproduction illuminated by chemical game theory, Proceedings of the National Academy of Sciences 113, 5030 (2016).
  * Note [4] Here, we consider the reproduction rate functions are symmetric between ${\rm{X}}_{1}$ and ${\rm{X}}_{2}$ for the sake of simplicity. Later we also consider the asymmetric case in the Sec. II E.
  * Note [5] Note that if $r(x)x$ is monotonic, this is also a necessary condition. However, in general, this is only a sufficient and not necessary condition. For example, there are cases where both symmetric and asymmetric trajectories are stable if $r(x)x$ is non-monotonic. For example, $r(x)x=\epsilon+\kappa x(x^{2}-\frac{3}{2}(\alpha+\beta)+3\alpha\beta)$.
  * Arsène  _et al._ [2018] S. Arsène, S. Ameta, N. Lehman, A. D. Griffiths, and P. Nghe, Coupled catabolism and anabolism in autocatalytic rna sets, Nucleic Acids Research 46, 9660 (2018).
  * Note [6] The bases M and N are the middle nucleotide of the 3nt recognition element in WXY, called the internal guide sequence (IGS) and the tag sequence [48].
  * Note [7] Here, the background reaction rate is approximated as a constant, although this reaction is due to catalyzed reaction by the non-covalent ribozymes or non-corresponding ribozymes. Assuming that the WXY fragments are abundant, the concentration of non-covalent ribozymes WXY:Z is $\sim s$ approximately. We assume their catalytic activity is non-specific, whose reaction rate is $\epsilon$. Next, we assume corresponding and non-corresponding ribozymes catalyze at efficiency $\tilde{\kappa}^{\prime}$ and $\tilde{\kappa}^{\prime}$, respectively. Then, the rate of the reaction $\rm{S}\rightarrow\rm{X_{1}}$ is $s\epsilon+x_{2}\tilde{\kappa}^{\prime}+x_{1}\tilde{\kappa}$. If we assume that $\epsilon=\tilde{\kappa}^{\prime}$, the rate is given as $\epsilon+\kappa x_{1}$, where $\kappa=\tilde{\kappa}_{1}-\epsilon$.
  * Note [8] Here, the time derivative of difference of concentration between $\rm{X_{1}}$ and $\rm{X_{2}}$ is $\frac{d}{dt}\left(\frac{\delta}{\chi}\right)=\frac{1}{\chi^{2}}(\frac{d\delta}{dt}\chi-\delta\frac{d\chi}{dt})=-2\epsilon s\frac{\delta}{\chi},$ where $\delta=x_{1}-x_{2}$ and $\chi=x_{1}+x_{2}$. Therefore, $\delta/\chi$ decays exponentially into zero in the characteristic relaxation time scale $\tau=\frac{1}{2s\epsilon}$. For example, if $s\epsilon\sim 0.01\min^{-1}$, the half time $\tau$ is estimated as $\sim 50\min$.
  * Note [9] The ACSs based on the _Azoarcus_ system show slow backward reactions since they are based on the recombination reactions of nucleotides [51].
  * Note [10] Moreover, the _Azoarcus_ ribozyme can form cross-catalytic networks [47, 30] (e.g., choosing bases M and N to be CC and GG, or AA and UU); in such cases, each network is a unit of self-reproduction (called an ‘autocatalytic core’) and could compete with other units.
  * Jeancolas  _et al._ [2021] C. Jeancolas, Y. J. Matsubara, M. Vybornyi, C. N. Lambert, A. Blokhuis, T. Alline, A. D. Griffiths, S. Ameta, S. Krishna, and P. Nghe, Rna diversification by a self-reproducing ribozyme revealed by deep sequencing and kinetic modelling, Chemical Communications 57, 7517 (2021).
  * Jain and Krishna [2002] S. Jain and S. Krishna, Crashes, recoveries, and “core shifts” in a model of evolving networks, Physical Review E 65, 026103 (2002).
  * Kaneko [2005] K. Kaneko, On recursive production and evolvabilty of cells: catalytic reaction network approach, Advanced in Chemical Physics 130, 543 (2005).
  * Note [11] This type of bistability: low and high catalyst concentration states, also could appear in our model if the background reaction rate is small enough (nearly zero). However, such bistability is easily destroyed under the SD protocol unless the interval is very short.
  * Kaneko and Yomo [2002] K. Kaneko and T. Yomo, On a kinetic origin of heredity: minority control in a replicating system with mutually catalytic molecules, Journal of Theoretical Biology 214, 563 (2002).
  * Doulcier  _et al._ [2020] G. Doulcier, A. Lambert, S. De Monte, and P. B. Rainey, Eco-evolutionary dynamics of nested darwinian populations and the emergence of community-level heredity, Elife 9, e53433 (2020).
  * Takeuchi  _et al._ [2020] N. Takeuchi, N. Mitarai, and K. Kaneko, Scaling laws of multilevel selection: a striking difference between continuous-trait and binary-trait models, arXiv (2020), arXiv:2005.04421.
  * Matsumura  _et al._ [2016] S. Matsumura, Á. Kun, M. Ryckelynck, F. Coldren, A. Szilágyi, F. Jossinet, C. Rick, P. Nghe, E. Szathmáry, and A. D. Griffiths, Transient compartmentalization of rna replicators prevents extinction due to parasites, Science 354, 1293 (2016).
  * Ameta  _et al._ [2022] S. Ameta, M. Kumar, N. Chakraborty, Y. J. Matsubara, S. Prashanth, D. Gandavadi, and S. Thutupalli, Multispecies autocatalytic rna reaction networks in coacervates, bioRxiv https://doi.org/10.1101/2022.11.01.514660 (2022).
  * Press  _et al._ [2007] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical recipes 3rd edition: The art of scientific computing_ (Cambridge university press, 2007).
  * Shirt-Ediss  _et al._ [2015] B. Shirt-Ediss, R. V. Solé, and K. Ruiz-Mirazo, Emergent chemical behavior in variable-volume protocells, Life 5, 181 (2015).
  * Ruiz-Herrero  _et al._ [2019] T. Ruiz-Herrero, T. G. Fai, and L. Mahadevan, Dynamics of growth and form in prebiotic vesicles, Physical review letters 123, 038102 (2019).
  * Luisi [1985] P. L. Luisi, Enzymes hosted in reverse micelles in hydrocarbon solution, Angewandte Chemie International Edition 24, 439 (1985).
  * Luisi and Varela [1989] P. L. Luisi and F. J. Varela, Self-replicating micelles— a chemical version of a minimal autopoietic system, Origins of Life and Evolution of the Biosphere 19, 633 (1989).
  * Gardiner  _et al._ [1985] C. W. Gardiner _et al._ ,  _Handbook of stochastic methods_ , Vol. 3 (springer Berlin, 1985).
  * Gillespie [2000] D. T. Gillespie, The chemical langevin equation, The Journal of Chemical Physics 113, 297 (2000).
  * Note [12] The dilution factor $m$ could also be randomly distributed, or depending on the random variable interval at each cycle $\varDelta t$. In these cases, the results do not qualitatively differ.

##  Appendix A Inheritable variety in general autocatalytic systems

###  A.1 The necessary condition under the serial dilution protocol

We consider the competing autocatalytic entities, $\rm{X}_{1}$ and
$\rm{X}_{1}$, under the SD protocol. The rate equation for the concentration
for the entities, $x_{1}$ and $x_{2}$, during one cycle until the dilution,
are Eq. 1 in the main text,

$\frac{d{x}_{i}}{dt}=s(\\{x_{j}\\},t)r(x_{i})x_{i},$ (S1)

where $s(\\{x_{j}\\},t)$ is the concentration of the substrate S that is
consumed in the replication reactions and $r(x_{i})$ is the reproduction rate
of $x_{i}$, respectively. Here we assume $s$ is symmetrical under exchange of
$x_{i}$s, i.e., $s(x_{1},x_{2})=s(x_{2},x_{1})$ and $s$, $x_{i}$s satisfy some
conservation law (e.g., $s+x_{1}+x_{2}=s^{tot}$). Also, we assume $s$ and $r$
are differentiable and non-negative functions for $x_{1}\geq 0$ and $x_{2}\geq
0$.

Here, we define $\chi=x_{1}+x_{2}$ and $\delta=x_{1}-x_{2}$, respectively.
Then, the time derivative of them are

$\frac{d\chi}{dt}=sr(\frac{\chi}{2})\chi+\mathcal{O}(\delta^{2}),\quad\frac{d\delta}{dt}=s\left(r(\frac{\chi}{2})+\frac{1}{2}\chi\frac{dr}{dx}(\frac{\chi}{2})\right)\delta+\mathcal{O}(\delta^{2}),$
(S2)

where we used the expansion
$r\left(\frac{1}{2}(\chi\pm\delta)\right)=r(\frac{\chi}{2})\pm\frac{dr}{dx}(\frac{\chi}{2})\frac{\delta}{2}+\mathcal{O}\left((\frac{\delta}{2})^{2}\right)$,
and assumed $\delta$ is small compared with $\chi$, i.e., the concentrations
of two catalysts are nearly equal, $x_{1}\sim x_{2}$. Then,

$\frac{d}{dt}\bigl{(}\frac{\delta}{\chi}\Bigr{)}=\frac{1}{\chi^{2}}(\frac{d\delta}{dt}\chi-\frac{d\chi}{dt}\delta)=s\left(\frac{1}{2}\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}\chi\right)\frac{\delta}{\chi}.$
(S3)

The integration of $\frac{d}{dt}(\frac{\delta}{\chi})/(\frac{\delta}{\chi})$
lead to

$\displaystyle\log\left|\frac{\delta(t)}{\chi(t)}\right|$
$\displaystyle=\int_{0}^{t}\frac{1}{2}s\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}\chi
dt$ $\displaystyle+$ $\displaystyle\log\left|\frac{\delta(0)}{\chi(0)}\right|$
(S4)
$\displaystyle=\int_{\frac{\chi(0)}{2}}^{\frac{\chi(t)}{2}}\frac{\frac{dr}{dt}(x)}{r(x)}dx$
$\displaystyle+$ $\displaystyle\log\left|\frac{\delta(0)}{\chi(0)}\right|,$

where we used $\frac{dt}{d\chi}=1/(sr(\frac{\chi}{2})\chi)$. Therefore,

$\frac{\delta(t)}{\chi(t)}=\frac{r(\frac{\chi(t)}{2})}{r(\frac{\chi(0)}{2})}\frac{\delta(0)}{\chi(0)},$
(S5)

which is Eq. 2 in the main text.

Now, we consider the serial dilution protocol, i.e., the amounts of $x_{1}$
and $x_{2}$ are multiplied by $m^{-1}$ at the end of a cycle $t=\varDelta t$.
In the stationary trajectory, $\chi(t)$ should satisfy the condition
$\chi(\varDelta t)=m\chi(0)$. If the condition

$\frac{\delta(\varDelta t)}{\chi(\varDelta t)}>\frac{\delta(0)}{\chi(0)}$ (S6)

is met, the difference between $x_{1}$ and $x_{2}$, $\delta/\chi$ is magnified
during a cycle. Therefore, the stationary trajectory with the equal
concentration of $\rm{X}_{1}$ and $\rm{X}_{2}$ (i.e., $\delta=0$) is unstable;
if otherwise, the stationary trajectory is stable.

Thus, surprisingly, whether the trajectory is stable or not is determined by
only whether the replication rate at the end of a cycle
$r(\frac{\chi(\varDelta t)}{2})$ is larger than that at the beginning
$r(\frac{\chi(0)}{2})$ or not. Roughly, the replication rates at the beginning
and the end are interpreted as the background and catalyzed reaction rates.
For example, if $r(x)=\frac{\epsilon+\kappa x}{x}$ (the system based on the
Azoarcus ribozyme) the stationary trajectory with $\delta=0$ is always stable.
If $r(x)=\frac{(\epsilon+\kappa x)^{2}}{x}$, $\varDelta t_{c}$ is calculated
as $\varDelta t_{c}=\log\left|\frac{\epsilon+\kappa s^{tot}}{\epsilon+\kappa
s^{tot}m^{-1}}\right|^{2}\sim 2\log\left(1+\frac{\kappa
s^{tot}}{\epsilon}\right)$.

###  A.2 The necessary condition for the heredity under the CSTR

Next, we consider the competing autocatalytic systems under the CSTR
condition, where the constant dilution rate is $\bar{\phi}$. Similar as in the
previous section, the rate equations for $x_{1}$ and $x_{2}$ are

$\frac{dx_{i}}{dt}=s(\\{x_{j}\\},t)r(x_{i})x_{i}-\bar{\phi}x_{i},$ (S7)

where the notations are the same as in the previous section.

Here, the time derivative of $\chi$ and $\delta$ ($\chi=x_{1}+x_{2}$ and
$\delta=x_{1}-x_{2}$) are

$\frac{d\chi}{dt}=sr\Bigl{(}\frac{\chi}{2}\Bigr{)}\chi-\bar{\phi}\chi+\mathcal{O}(\delta^{2}),\quad\frac{d\delta}{dt}=s\left(r\Bigl{(}\frac{\chi}{2}\Bigr{)}+\frac{1}{2}\chi\frac{dr}{dt}\Bigl{(}\frac{\chi}{2}\Bigr{)}\right)\delta-\bar{\phi}\delta+\mathcal{O}(\delta^{2}),$
(S8)

where we used
$r\left(\frac{1}{2}(\chi\pm\delta)\right)=r(\frac{\chi}{2})\pm\frac{dr}{dx}(\frac{\chi}{2})\frac{\delta}{2}+O\left((\frac{\delta}{2})^{2}\right)$,
and assumed $\delta$ is small. The same calculation in the previous section
leads to

$\frac{d}{dt}\Bigl{(}\frac{\delta}{\chi}\Bigr{)}=\frac{1}{\chi^{2}}(\frac{d\delta}{dt}\chi-\frac{d\chi}{dt}\delta)=s\left(\frac{1}{2}\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}\chi\right)\frac{\delta}{\chi}+\mathcal{O}((\delta/\chi)^{2}).$
(S9)

Thus, given the steady state concentration $\chi^{*}$ such that
$s^{*}r(\frac{\chi^{*}}{2})=\bar{\phi}$, the state with $\delta=0$ is stable
if $\frac{dr}{dx}(\frac{\chi^{*}}{2})$ is negative. For example, if
$r(x)=\frac{\epsilon+\kappa x}{x}$ (the _Azoarcus_ based system with only one
step), $\frac{dr}{dx}(x)=-\epsilon/x^{2}$ is negative for all $x$, therefore
the state with $\delta=0$ is always stable. While in a case with
$r(x)=\frac{\epsilon+\kappa x^{2}}{x}$,
$\frac{dr}{dx}(x)=-\epsilon/x^{2}+\kappa$. Therefore, if
$\bar{\phi}>\phi_{c}^{cstr}$, where
$\phi_{c}^{cstr}=s^{tot}\sqrt{\epsilon\kappa}-4\epsilon$, the state with
$\delta=0$ is stable.

Note that, for the stability of the steady state $\chi^{*}$, the condition
$-r(\frac{\chi^{*}}{2})+\frac{1}{2}(s^{tot}-\chi^{*})\frac{dr}{dx}(\frac{\chi^{*}}{2})<0$
should be satisfied. The $\chi$-direction is always stable if
$\delta$-direction is stable, i.e., $\frac{dr}{dx}(\frac{\chi^{*}}{2})<0$.

###  A.3 Under the general dilution protocol (compartment growth and division)

Lastly, we discuss the heredity of the system under the general dilution
scenarios:

$\frac{dx_{i}}{dt}=s(\\{x_{j}\\},t)r(x_{i})x_{i}-\phi(t)x_{i},$ (S10)

where $\phi(t)$ is the time-dependent dilution rate due to the growth of the
compartment volume $V$, $\phi(t)=\frac{dV}{dt}/V$. Here we restrict $\phi(t)$
to periodic functions $\phi(t+\varDelta t)=\phi(t)$ such that
$\int_{0}^{\varDelta t}\phi(t)dt=\bar{\phi}\varDelta t$. Also
$\sigma(t)=s^{tot}\phi(t)$, thus the total mass of the substrate is kept as a
constant, $s+\sum_{i}x_{i}=s^{tot}$. We further assume that $r(x)$ is a convex
function, i.e., $\frac{d^{2}r}{dt^{2}}(x)>0$.

Under the above setup, the critical interval $\varDelta t_{c}$, which divides
the region with and without heredity, is bounded by both that under the serial
dilution and the CSTR:

  1. 1.

The period $\varDelta t_{c}$ for the protocol under any $\phi(t)$ is bounded
from the above by that under the CSTR $\varDelta t_{c}^{cstr}$, $\varDelta
t_{c}\leq\varDelta t_{c}^{cstr}$. Here, $\varDelta t_{c}^{cstr}$ is either
$\infty$ or $0$.

  2. 2.

$\varDelta t_{c}$ is bounded from the bottom by that under the serial dilution
$\varDelta t_{c}^{sd}$, $\varDelta t_{c}^{sd}\leq\varDelta t_{c}$.

To prove this (non-rigorously), we use the dilution rate function $\phi(t)$

$\phi(t)\equiv\begin{cases}\bar{\phi}_{1}\quad(0\leq t\leq\tau),\\\
\bar{\phi}_{2}\quad(\tau\leq t\leq\varDelta t).\\\ \end{cases}$ (S11)

Figure S1: Schematics for the proof in Sec. A.3. Left figures represent the
time-dependent dilution rate in one dilution cycle, $\phi(t)$, and right
figures represent $\chi(t)$. The top figures are for the perturbation from the
CSTR condition (The dashed lines represent the corresponding case under the
CSTR condition). The bottom figures are for the perturbation from the SD
condition (The dashed lines represent the corresponding SD).

##### The perturbation from the CSTR condition:

Firstly, we consider the perturbation from the CSTR protocol fixing
$\int_{0}^{\varDelta t}\phi(t)=\bar{\phi}\varDelta t$, i.e.,
$\bar{\phi}_{1}=\bar{\phi}+\frac{\xi}{\tau}$ and
$\bar{\phi}_{2}=\bar{\phi}-\frac{\xi}{\varDelta t-\tau}$ (see Fig. S1). Here,
we define the steady-state concentration under the CSTR with the constant
dilution rate $\bar{\phi}$, $\chi^{*}$, is such that
$(s^{tot}-\chi^{*})r(\chi^{*})\chi^{*}-\bar{\phi}=0$. Here we consider a small
perturbation for $\bar{\phi}$, $\bar{\phi}\rightarrow\bar{\phi}+\xi_{i}$,
where $\xi_{1}=\frac{\xi}{\tau}$ and $\xi_{2}=-\frac{\xi}{\varDelta t-\tau}$,
and the deviation for $\chi^{*}$ due to it,
$\chi^{*}\rightarrow\chi^{*}+\eta_{i}$, which should satisfies

$\displaystyle(s^{tot}-\chi^{*}-\eta_{i})r\Bigl{(}\frac{\chi^{*}+\eta_{i}}{2}\Bigr{)}-(\bar{\phi}+\xi_{i})=0,$
(S12)
$\displaystyle(s^{tot}-\chi^{*}-\eta_{i})\left(r\Bigl{(}\frac{\chi^{*}}{2}\Bigr{)}+\frac{\eta_{i}}{2}\frac{dr}{dx}\Bigl{(}\frac{\chi^{*}}{2}\Bigr{)}+\frac{\eta_{i}^{2}}{8}\frac{d^{2}r}{dx^{2}}\Bigl{(}\frac{\chi^{*}}{2}\Bigr{)}+\mathcal{O}(\eta_{i}^{3})\right)-(\bar{\phi}+\xi_{i})=0,$
$\displaystyle\eta_{i}\left(-r\Bigl{(}\frac{\chi^{*}}{2}\Bigr{)}+(s^{tot}-\chi^{*})\frac{1}{2}\frac{dr}{dx}\Bigl{(}\frac{\chi^{*}}{2}\Bigr{)}\right)+\eta_{i}^{2}\left(-\frac{1}{2}\frac{dr}{dx}\Bigl{(}\frac{\chi^{*}}{2}\Bigr{)}+(s^{tot}-\chi^{*})\frac{1}{8}\frac{d^{2}r}{dx^{2}}\Bigl{(}\frac{\chi^{*}}{2}\Bigr{)}\right)+\mathcal{O}(\eta_{i}^{3})-\xi_{i}=0,$

where we used the expansion
$r\bigl{(}\frac{\chi^{*}+\eta_{i}}{2}\bigr{)}=r\bigl{(}\frac{\chi^{*}}{2}\bigr{)}+\frac{\eta_{i}}{2}\frac{dr}{dx}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}+\frac{\eta_{i}^{2}}{8}\frac{d^{2}r}{dx^{2}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}+\mathcal{O}(\eta_{i}^{3})$,
and $(s^{tot}-\chi^{*})r(\chi^{*})-\bar{\phi}=0$. Then, $\eta_{i}$ is
determined as

$\eta_{i}=-\tilde{\xi}_{i}-\frac{1}{2}\mathcal{R}\tilde{\xi}_{i}^{2}+\mathcal{O}(\xi^{3}),$
(S13)

where we define

$\displaystyle\tilde{\xi}_{i}=\frac{\xi_{i}}{r(\frac{\chi^{*}}{2})-\frac{1}{2}(s^{tot}-\chi^{*})\frac{dr}{dx}(\frac{\chi^{*}}{2})},\quad\mathcal{R}=\frac{\frac{dr}{dx}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}-\frac{1}{4}(s^{tot}-\chi^{*})\frac{d^{2}r}{dx^{2}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}}{r\bigl{(}\frac{\chi^{*}}{2}\bigr{)}-\frac{1}{2}(s^{tot}-\chi^{*})\frac{dr}{dx}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}}.$
(S14)

Note that, from the condition for the stability of the steady state
$\chi^{*}$,
$\frac{d}{dt}\left((s^{tot}-x)r(\frac{x}{2})-\bar{\phi}\right)|_{x=\chi^{*}}<0$,
$\frac{1}{2}(s^{tot}-\chi^{*})\frac{dr}{dx}(\frac{\chi^{*}}{2})-r(\frac{\chi^{*}}{2})<0$.

The same as the previous sections, the integration of
$\frac{d}{dt}(\frac{\delta}{\chi})/(\frac{\delta}{\chi})$ from $t=0$ to
$=\varDelta t$ lead to

$\displaystyle\log\left|\frac{\delta(\varDelta t)}{\chi(\varDelta
t)}\right|-\log\left|\frac{\delta(0)}{\chi(0)}\right|=\frac{1}{2}\int_{0}^{\varDelta
t}(s^{tot}-\chi)\chi\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}dt,$ (S15)
$\displaystyle\approx\tau(s^{tot}-\chi^{*}-\eta_{1})\frac{\chi^{*}+\eta_{1}}{2}\frac{dr}{dx}\Bigl{(}\frac{\chi^{*}+\eta_{1}}{2}\Bigr{)}+(\varDelta
t-\tau)(s^{tot}-\chi^{*}-\eta_{2})\frac{\chi^{*}+\eta_{2}}{2}\frac{dr}{dx}\Bigl{(}\frac{\chi^{*}+\eta_{2}}{2}\Bigr{)},$

where we assumed that when $\phi$ is changed (at $t=0$ and $t=\tau$), the
relaxation time to the steady concentration is much smaller than $\tau$ and
$\varDelta t-\tau$ (see Fig. S1), whose contribution is of the order of
$\xi^{3}$. Then, substituting Eq. S13 into the above,

$\displaystyle=$ $\displaystyle\varDelta
t(s^{tot}-\chi^{*})\frac{\chi^{*}}{2}\frac{dr}{dx}\Bigl{(}\frac{\chi^{*}}{2}\Bigr{)}$
(S16) $\displaystyle+(\tau\eta_{1}+(\varDelta
t-\tau)\eta_{2})\left((\frac{1}{2}s^{tot}-\chi^{*})\frac{dr}{dx}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}+(s^{tot}-\chi^{*})\frac{\chi^{*}}{4}\frac{d^{2}r}{dx^{2}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}\right)$
$\displaystyle+(\tau\eta_{1}^{2}+(\varDelta
t-\tau)\eta_{2}^{2})\left((s^{tot}-2\chi^{*})\frac{1}{4}\frac{d^{2}r}{dx^{2}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}-\frac{1}{2}\frac{dr}{dx}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}+(s^{tot}-\chi^{*})\frac{\chi^{*}}{16}\frac{d^{3}r}{dx^{3}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}\right)+\mathcal{O}(\xi^{3}),$
$\displaystyle=$ $\displaystyle\varDelta
t(s^{tot}-\chi^{*})\frac{\chi^{*}}{2}\frac{dr}{dx}\Bigl{(}\frac{\chi^{*}}{2}\Bigr{)}-\frac{1}{2}\xi^{2}\Bigl{(}\frac{1}{\tau}+\frac{1}{\varDelta
t-\tau}\Bigr{)}\Biggl{[}\mathcal{R}\left((\frac{1}{2}s^{tot}-\chi^{*})\frac{dr}{dx}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}+(s^{tot}-\chi^{*})\frac{\chi^{*}}{4}\frac{d^{2}r}{dx^{2}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}\right)$
$\displaystyle+(\chi^{*}-\frac{1}{2}s^{tot})\frac{d^{2}r}{dx^{2}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}+\frac{dr}{dx}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}-(s^{tot}-\chi^{*})\frac{\chi^{*}}{8}\frac{d^{3}r}{dx^{3}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}\Biggr{]}+\mathcal{O}(\xi^{3}),$

where the first term is the deviation under the CSTR with the dilution rate
$\bar{\phi}$, $\log\left[\frac{\delta^{cstr}(\varDelta
t)}{\chi^{cstr}(\varDelta t)}/\frac{\delta^{cstr}(0)}{\chi^{cstr}(0)}\right]$.
When $\frac{dr}{dx}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}=0$,

$=-\frac{1}{2}\xi^{2}\Bigl{(}\frac{1}{\tau}+\frac{1}{\varDelta
t-\tau}\Bigr{)}\left((\chi^{*}-\frac{1}{2}s^{tot})\frac{d^{2}r}{dx^{2}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}-\frac{1}{16}(s^{tot}-\chi^{*})^{2}{\chi^{*}}\frac{\bigl{(}\frac{d^{2}r}{dx^{2}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}\bigr{)}^{2}}{r\bigl{(}\frac{\chi^{*}}{2}\bigr{)}}-\frac{1}{8}(s^{tot}-\chi^{*})\chi^{*}\frac{d^{3}r}{dx^{3}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}\right)+\mathcal{O}(\xi^{3}).$
(S17)

Here, if we approximate that $\chi^{*}$ is nearly saturated, i.e.,
$\chi^{*}\sim s^{tot}$,

$\sim-\frac{1}{4}\xi^{2}\Bigl{(}\frac{1}{\tau}+\frac{1}{\varDelta
t-\tau}\Bigr{)}s^{tot}\frac{d^{2}r}{dx^{2}}\bigl{(}\frac{\chi^{*}}{2}\bigr{)}+\mathcal{O}(\xi^{3}).$
(S18)

Therefore, under the $\phi(t)$ with any choice of $\xi$ and $\tau$,
$\frac{\delta(\varDelta t)}{\chi(\varDelta
t)}/\frac{\delta(0)}{\chi(0)}<\frac{\delta^{cstr}(\varDelta
t)}{\chi^{cstr}(\varDelta t)}/\frac{\delta^{cstr}(0)}{\chi^{cstr}(0)}$. Thus,
if the stationary trajectory with $\delta=0$ is stable under the CSTR, i.e.,
$\frac{\delta^{cstr}(\varDelta t)}{\chi^{cstr}(\varDelta
t)}/\frac{\delta^{cstr}(0)}{\chi^{cstr}(0)}<1$, it is also true for that under
the protocol with the dilution rate $\phi(t)$, therefore, $\varDelta
t_{c}\leq\varDelta t_{c}^{cstr}$.

Further, we can divide the region for $\bar{\phi}_{1}$ or $\bar{\phi}_{2}$,
and add more steps for the function $\phi(t)$. The same analysis as the above
reveals that in each addition $\frac{\delta(\varDelta t)}{\chi(\varDelta
t)}/\frac{\delta(0)}{\chi(0)}$ declines monotonically, thus the upper limit
for heredity $\varDelta t_{c}$ also declines monotonically.

##### The perturbation from the serial dilution protocol:

Secondly, we also consider the perturbation from the case with the serial
dilution; here we consider $\bar{\phi}_{1}$ and $\bar{\phi}_{2}$ as
$\bar{\phi}_{1}=0$ and $\bar{\phi}_{2}=\bar{\phi}\frac{\varDelta t}{\varDelta
t-\tau}$ (see Fig. S1). Note that when $\varDelta t-\tau\rightarrow 0$, the
dilution protocol becomes the same as the serial dilution.

The deviation of $\frac{\delta}{\chi}$ in one cycle under the dilution
protocol with $\phi(t)$,

$\displaystyle\log\left|\frac{\delta(\varDelta t)}{\chi(\varDelta
t)}\right|-\log\left|\frac{\delta(0)}{\chi(0)}\right|=\frac{1}{2}\int_{0}^{\varDelta
t}(s^{tot}-\chi)\chi\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}dt,$ (S19)
$\displaystyle=\int_{\frac{\chi(0)}{2}}^{\frac{\chi(\tau)}{2}}\frac{\frac{dr}{dt}}{r}dx+\frac{1}{2}\int_{\tau}^{\varDelta
t}(s^{tot}-\chi)\chi\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}dt$

We assume $\chi^{sd}(\varDelta t)\sim\chi(\tau)$, because in the both cases
$\chi(t)$ is saturated (i.e., S is exhausted) if $\varDelta t$ is enough large
(see Fig. S1). Then, $\chi^{sd}(0)=\chi(\tau)\exp(-\bar{\phi}\varDelta t)$.
From the rate equation
$\frac{d\chi}{dt}=\left((s^{tot}-\chi)r\bigl{(}\frac{\chi}{2}\bigr{)}-\bar{\phi}\right)\chi$

$\displaystyle\chi(\varDelta t)=\chi(0)=\chi(\tau)\exp\left(-\varDelta
t\bar{\phi}+\int_{\tau}^{\varDelta
t}(s^{tot}-\chi)r\Bigl{(}\frac{\chi}{2}\Bigr{)}dt\right),$ (S20)
$\displaystyle=\chi^{sd}(0)\exp\left(\int_{\tau}^{\varDelta
t}(s^{tot}-\chi)r\Bigl{(}\frac{\chi}{2}\Bigr{)}dt\right).$

On the other hand, the deviation of $\frac{\delta^{sd}}{\chi^{sd}}$ in one
cycle under SD is derived as

$\log\left|\frac{\delta^{sd}(\varDelta t)}{\chi^{sd}(\varDelta
t)}\right|-\log\left|\frac{\delta^{sd}(0)}{\chi^{sd}(0)}\right|=\int_{\frac{\chi^{sd}(0)}{2}}^{\frac{\chi(0)}{2}}\frac{\frac{dr}{dx}}{r}dx+\int_{\frac{\chi(0)}{2}}^{\frac{\chi(\tau)}{2}}\frac{\frac{dr}{dx}}{r}dx,$
(S21)

Here, the first term is,

$\int_{\frac{\chi^{sd}(0)}{2}}^{\frac{\chi(0)}{2}}\frac{\frac{dr}{dx}}{r}dx=\log\left|\frac{r\Bigl{(}\frac{\chi(0)}{2}\Bigr{)}}{r\Bigl{(}\frac{\chi^{sd}(0)}{2}\Bigr{)}}\right|=\frac{\chi^{sd}(0)}{2}\frac{\frac{dr}{dx}\Bigl{(}\frac{\chi^{sd}(0)}{2}\Bigr{)}}{r\Bigl{(}\frac{\chi^{sd}(0)}{2}\Bigr{)}}\int_{\tau}^{\varDelta
t}(s^{tot}-\chi)r\Bigl{(}\frac{\chi}{2}\Bigr{)}dt+\mathcal{O}((\varDelta
t-\tau)^{2}),$ (S22)

where in the approximation, we used the assumption that $\varDelta t-\tau$ is
small.

From the comparison between Eq. S21 Eq. S19,

$\displaystyle\log\left|\frac{\delta(\varDelta t)}{\chi(\varDelta
t)}\right|-\log\left|\frac{\delta(0)}{\chi(0)}\right|-\left(\log\left|\frac{\delta^{sd}(\varDelta
t)}{\chi^{sd}(\varDelta
t)}\right|-\log\left|\frac{\delta^{sd}(0)}{\chi^{sd}(0)}\right|\right)$ (S23)
$\displaystyle=\frac{1}{2}\int_{\tau}^{\varDelta
t}(s^{tot}-\chi(t^{\prime}))\left[\chi(t^{\prime})\frac{dr}{dx}\Bigl{(}\frac{\chi(t^{\prime})}{2}\Bigr{)}-\frac{r\Bigl{(}\frac{\chi(t^{\prime})}{2}\Bigr{)}\chi^{sd}(0)\frac{dr}{dx}\Bigl{(}\frac{\chi^{sd}(0)}{2}\Bigr{)}}{r\Bigl{(}\frac{\chi^{sd}(0)}{2}\Bigr{)}}\right]dt^{\prime}.$

Thus, if
$\frac{\chi(t)\frac{dr}{dx}\bigl{(}\frac{\chi(t)}{2}\bigr{)}}{r\bigl{(}\frac{\chi(t)}{2}\bigr{)}}>\frac{\chi^{sd}(0)\frac{dr}{dx}\bigl{(}\frac{\chi^{sd}(0)}{2}\bigr{)}}{r\bigl{(}\frac{\chi^{sd}(0)}{2}\bigr{)}}$
is satisfied, the above is always positive. As $\chi(t)>\chi^{sd}(0)$, it
reveals that $\frac{\delta(\varDelta t)}{\chi(\varDelta
t)}/\frac{\delta(0)}{\chi(0)}>\frac{\delta^{sd}(\varDelta
t)}{\chi^{sd}(\varDelta t)}/\frac{\delta^{sd}(0)}{\chi^{sd}(0)}$
$\frac{\delta(\varDelta t)}{\chi(\varDelta t)}$ for the serial dilution is
always less than that for the protocol with the dilution rate function
$\phi(t)$, if $\frac{x\frac{dr}{dt}(x)}{r(x)}$ is an increasing function of
$x$. Therefore, the critical period of time for $\phi(t)$, $\varDelta t_{c}$
is always larger than that for the serial dilution, $\varDelta t^{sd}_{c}$,
i.e., $\varDelta t^{sd}_{c}\leq\varDelta t_{c}$.

##  Appendix B robustness of results

Figure S2: Schematics of the alternative models for the competing
autocatalytic entities, $\rm{X}_{1}$ and $\rm{X}_{2}$; they are converted from
a substrate S (solid arrows) catalyzed by itself (dashed arrows). (A) non-
modified model (the same as Fig. 3A in the main text) (B) with reversible
chemical reactions (Sec. B.1) (C) with asymmetric kinetics rates (Sec. B.2)
(D) with more general reproduction rate function (dashed arrow with bar head
represents the inhibition; Sec. B.3)

###  B.1 Autocatalytic sets with reversible chemical reactions

In the main text, we mainly considered chemical reaction systems with
irreversible reactions. However, all reactions should be reversible to be
chemically consistent and converge to thermal equilibrium in the absence of
dilution protocols. Here, we expand our previous analysis to examine how
reversible reactions change the phase diagram of bistability.

If the reactions in the system in Fig. 3A is reversible, the rate equations
are

$\frac{d{x}_{i}}{dt}=(s-bx_{i})r(x_{i})x_{i},$ (S24)

where $b$ is the relative rate of the backward reaction (see Fig. S2).

Then, $\frac{d\chi}{dt}$ and $\frac{d\delta}{dt}$ are modified as

$\displaystyle\frac{d\chi}{dt}=$
$\displaystyle(s-\frac{b}{2}\chi)r(\frac{\chi}{2})\chi$ $\displaystyle+$
$\displaystyle\mathcal{O}(\delta^{2}),$ (S25)
$\displaystyle\frac{d\delta}{dt}=$ $\displaystyle
s\left(r(\frac{\chi}{2})+\frac{1}{2}\chi\frac{dr}{dt}(\frac{\chi}{2})\right)\delta-b\left(r(\frac{\chi}{2})+\frac{1}{4}\chi\frac{dr}{dt}(\frac{\chi}{2})\right)\chi\delta$
$\displaystyle+$ $\displaystyle\mathcal{O}(\delta^{2}).$

Then,

$\frac{d}{dt}\bigl{(}\frac{\delta}{\chi}\Bigr{)}=\left[(s-\frac{b}{2}\chi)(\frac{1}{2}\frac{dr}{dt}\Bigl{(}\frac{\chi}{2}\Bigr{)}\chi)-\frac{1}{2}b(r(\frac{\chi}{2})\chi)\right]\frac{\delta}{\chi}.$
(S26)

As a result,

$\log\left|\frac{\delta(t)}{\chi(t)}\right|=\left[\log(r(x))+\frac{b}{2+b}\log(s^{tot}-(2+b)x)\right]_{x=\frac{\chi(0)}{2}}^{x=\frac{\chi(t)}{2}}+\log\left|\frac{\delta(0)}{\chi(0)}\right|,$
(S27)

where we used $s=s^{tot}-\chi$. Thus, the reversible reactions alter Eq. 2 in
the main text as

$\frac{\delta(t)}{\chi(t)}=\frac{r(\frac{\chi(t)}{2})}{r(\frac{\chi(0)}{2})}\left(\frac{s^{tot}-(1+\frac{b}{2})\chi(t)}{s^{tot}-(1+\frac{b}{2})\chi(0)}\right)^{\frac{b}{2+b}}\frac{\delta(0)}{\chi(0)}.$
(S28)

Note that, from Eq. S25, $s^{tot}-(1+\frac{b}{2})\chi(\infty)=0$ at the
equilibrium, thus $\frac{\delta(t)}{\chi(t)}$ eventually converges to 0.

Similar to the case with irreversible reaction ($b=0$) in the main text, the
system shows bifurcation that the bistability disappears when varying the
cycle interval $\varDelta t$ or the dilution rate $\bar{\phi}$. Interestingly,
unlike the irreversible reaction case ($b=0$), the region of bistability is
bounded in the reversible case ($b>0$) for the parameter $\bar{\phi}$. Thus,
for the SD protocol with fixed $\varDelta t$, there is both an upper and a
lower critical $\bar{\phi}$ (see Fig. S3). Note that the boundary of the
region is close to the irreversible case when $\bar{\phi}$ is large and the
critical value depends on $\epsilon$ similarly to the $b=0$ case. Importantly,
even with reversible reactions, the parameter space for the general GD
protocol contains that for SD protocols, as we found for irreversible
reactions.

Figure S3: Phase space of the bistability in the case with reversible
reactions (Eq. S25). (A) The colored lines represent the boundary between with
or without bistability in different backward reaction rates $b$. The other
parameters are set as $\epsilon=0.5$, $\kappa=8$ and $\bar{\phi}=1$. (B) The
boundary between with or without bistability under various functions of the
dilution rates $\phi(t)$. We set $b=0.1$.

###  B.2 ACSs with asymmetric kinetic rates

In the previous sections, we considered symmetric competing ACSs. However, our
results are similar even if the two competing ACSs have different kinetic rate
constants (catalytic efficiency), although the bifurcation where the
bistability disappears is discontinuous.

In Eq. 1 in the main text, we assumed the reproduction rate function,
$r(x_{i})$ is the same between the two self-reproducing entities. Here, we
also discuss the rate functions of two entities, $r_{1}(x_{1})$ and
$r_{2}(x_{2})$ are different, i.e., $r_{1}(x)\neq r_{2}(x)$. The simplest
example is $r_{i}(x)x=\epsilon+\kappa_{i}x^{2}$, where
$\kappa_{1}\neq\kappa_{2}$. In this case, if the difference between the
catalytic strength of two entities $|\kappa_{1}-\kappa_{2}|$ is not so large,
the system has the bistability, as the nullclines shows (Fig. S4A). Further,
as the bifurcation diagram shows (Fig. S4B), even in this case, the critical
point for the interval $\varDelta t$ and the dilution rate $\bar{\phi}$ exist,
although the transition is discontinuous.

In such cases, the boundary in the parameter space for ACSs that have
bistability in the general GD protocols is also bounded by the boundaries for
the SD and CSTR protocols.

Figure S4:  (A) The red and blue curves represents the nullclines for the
composition of the entities $\rm{X}_{1}$ and $\rm{X}_{2}$ at just before the
dilution, $x_{1}(-0)$ and $x_{2}(-0)$, respectively, in a case with the
reproduction rate function are asymmetrical,
$r_{i}(x)x=\epsilon+\kappa_{i}x^{2}$ ($i=1,2$) and $\kappa_{1}<\kappa_{2}$).
The cross points represents the stable/unstable fixed points of concentrations
of $\rm{X}_{1}$ and $\rm{X}_{2}$ at just before the dilution, $x_{1}^{*}(-0)$
and $x_{2}^{*}(-0)$. See Methods and Models for the detail of the drawing of
the nullclines. We define $\kappa_{1}=\kappa-\frac{\iota}{2}$ and
$\kappa_{2}=\kappa+\frac{\iota}{2}$, and set $\varDelta t=1$, $\kappa=8$,
$\epsilon=0.5$, $\bar{\phi}=1$, and $\kappa_{2}-\kappa_{1}=\iota=0.5$ (left)
and $2$ (right).(B) The bifurcation diagram of the concentrations of the
entities $\rm{X}_{1}$ and $\rm{X}_{2}$ with a varying period of the dilution
cycles $\varDelta t$. In contrast with the symmetric case (i.e., $\iota=0$),
the $\rm{X}_{1}$-dominant state (and thus the bistability) disappears
discontinuously at around $\varDelta t_{c}^{sd}\sim 1.0$. (C) The dependence
of the critical period $\varDelta t_{c}$ of the cycle on $\bar{\phi}$.

###  B.3 More general reproduction rate function

Here, we consider the autocatalytic system is under the SD protocol, in which
the reproduction rate function depends on both $x_{1}$ and $x_{2}$, including
cases where $\rm{X}_{1}$ and $\rm{X}_{1}$ mutually inhibit their synthesis
(see Fig. S2D). For example, in the case of the genetic toggle switch model
[32], $r(x_{i},x_{\bar{i}})=\frac{1}{x_{i}}(\frac{1}{1+x_{\bar{i}}^{n}})$,
where $n$ is the Hill coefficient.

The rate equations are

$\frac{d{x}_{i}}{dt}=sr(x_{i},x_{-i})x_{i}.$ (S29)

where $-i$ represents 2 or 1, if $i=1$ or 2, respectively.

Here, the same as before, we define $\chi=x_{1}+x_{2},\delta=x_{1}-x_{2}$,
then

$\frac{d\chi}{dt}=sr(\frac{\chi}{2},\frac{\chi}{2})\chi+\mathcal{O}(\delta^{2}),\quad\frac{d\delta}{dt}=s\left(r(\frac{\chi}{2},\frac{\chi}{2})+\frac{1}{2}\chi\frac{\partial
r}{\partial
x_{i}}(\frac{\chi}{2},\frac{\chi}{2})-\frac{1}{2}\chi\frac{\partial
r}{\partial
x_{-i}}(\frac{\chi}{2},\frac{\chi}{2})\right)\delta+\mathcal{O}(\delta^{2}),$
(S30)

where we used the expansion

$r\left(\frac{1}{2}(\chi\pm\delta),\frac{1}{2}(\chi\mp\delta)\right)=r(\frac{\chi}{2},\frac{\chi}{2})\pm\frac{\partial
r}{\partial
x_{i}}(\frac{\chi}{2},\frac{\chi}{2})\frac{\delta}{2}\mp\frac{\partial
r}{\partial
x_{-i}}(\frac{\chi}{2},\frac{\chi}{2})\frac{\delta}{2}+\mathcal{O}(\delta^{2}),$

and assumed $\delta$ is small.

$\frac{d}{dt}\bigl{(}\frac{\delta}{\chi}\Bigr{)}=\frac{1}{\chi^{2}}(\frac{d\delta}{dt}\chi-\frac{d\chi}{dt}\delta)=\frac{1}{2}s\chi\left(\frac{\partial
r}{\partial x_{i}}(\frac{\chi}{2},\frac{\chi}{2})-\frac{\partial r}{\partial
x_{-i}}(\frac{\chi}{2},\frac{\chi}{2})\right)\frac{\delta}{\chi}.$ (S31)

The integration of $\frac{d}{dt}(\frac{\delta}{\chi})/(\frac{\delta}{\chi})$
lead to

$\displaystyle\log\left|\frac{\delta(t)}{\chi(t)}\right|$
$\displaystyle=\int_{0}^{t}\frac{1}{2}s\chi\left(\frac{\partial r}{\partial
x_{i}}(\frac{\chi}{2},\frac{\chi}{2})-\frac{\partial r}{\partial
x_{-i}}(\frac{\chi}{2},\frac{\chi}{2})\right)dt$ $\displaystyle+$
$\displaystyle\log\left|\frac{\delta(0)}{\chi(0)}\right|$ (S32)
$\displaystyle=\int_{\frac{\chi(0)}{2}}^{\frac{\chi(t)}{2}}\frac{\frac{\partial
r}{\partial x_{i}}(x,x)-\frac{\partial r}{\partial x_{-i}}(x,x)}{r(x,x)}dx$
$\displaystyle+$ $\displaystyle\log\left|\frac{\delta(0)}{\chi(0)}\right|,$

where we used $\frac{dt}{d\chi}=1/(sr(\frac{\chi}{2},\frac{\chi}{2})\chi)$.
Therefore,

$\frac{\delta(t)}{\chi(t)}=\frac{R_{i}(\frac{\chi(t)}{2})}{R_{i}(\frac{\chi(0)}{2})}\frac{R_{-i}(\frac{\chi(0)}{2})}{R_{-i}(\frac{\chi(t)}{2})}\frac{\delta(0)}{\chi(0)},$
(S33)

where $R_{i}(x)$ and $R_{-i}(x)$ are defined as
$R_{i}(x)=\exp\left(\int^{x}\frac{\frac{\partial r}{\partial
x_{i}}(x^{\prime},x^{\prime})}{r(x^{\prime},x^{\prime})}dx^{\prime}\right)$
and $R_{\bar{i}}(x)=\exp\left(\int^{x}\frac{\frac{\partial r}{\partial
x_{-i}}(x^{\prime},x^{\prime})}{r(x^{\prime},x^{\prime})}dx^{\prime}\right)$.
Therefore, the sufficient condition for the bistability of the system is

$\frac{R_{i}(\frac{\chi(t)}{2})}{R_{-i}(\frac{\chi(t)}{2})}>\frac{R_{i}(\frac{\chi(0)}{2})}{R_{-i}(\frac{\chi(0)}{2})}.$
(S34)

Here, if we assume that the separation of variables is possible for the
reproduction rate function $r(x_{i},x_{-i})$,
$r(x_{i},x_{-i})=f(x_{i})g(x_{-i})$; then, $R_{i}(x)=f(x)$ and
$R_{-i}(x)=g(x)$. Note that, if $g(x)=1$, then the condition Eq. S34 reduced
to Eq. 3 in the main text. For example, in a case with the genetic toggle
switch model, the above condition is $\frac{1+(\frac{\chi(\varDelta
t)}{2})^{n}}{\frac{\chi(\varDelta
t)}{2}}>\frac{1+(\frac{\chi(0)}{2})^{n}}{\frac{\chi(0}{2}}$. Thus this model
could show the bistability if $n>1$.

###  B.4 Robustness to the stochasticity in reaction systems or protocols

####  B.4.1 stochasticity in the reaction dynamics and transition between
states

If the number of molecules is small (e.g., the reaction dynamics occur inside
sufficiently small compartments), stochastic fluctuations are not negligible.
Close to the deterministic bifurcation transitions, these fluctuations cause
random transitions between states, which destroys the information to be
inherited. Thus, the parameter space where the system has heredity is narrower
than the deterministic case. However, our main results – the existence of a
critical threshold in the dilution interval or the dilution rate and the
bistable parameter region for general GD cycles being bounded by that of SD
and CSTR – are robust to the addition of stochastic noise in the chemical
reaction system.

Previously we discussed the deterministic model, Eq. 1. If the volume of the
compartment is small (i.e., the system size) the system, the stochastic
fluctuation in the reaction dynamics is non-negligible large [71]. The
chemical Langevin equations [71, 72], which is corresponding to Eq. 1 is,

$\frac{dx_{i}}{dt}=sr_{i}(x_{i})x_{i}+\sqrt{sr_{i}(x_{i})V^{-1}}\eta_{i}(t),$
(S35)

where $\eta_{i}$ are i.i.d. Gaussian random variables with the correlation
function $\langle\eta_{i}(t)\eta_{i}(t^{\prime})\rangle=\delta(t-t^{\prime})$.
$V$ is the volume of the compartment, and $V^{-1}$ corresponds to the
intensity of noise.

Here, we assume the serial dilution protocol (the volume $V$ is fixed during
one cycle). In contrast with the deterministic case, the transition could
occur from $\rm X_{1}$($\rm X_{2}$)-dominant state to the other one in the
presence of the stochastic noise (Fig. S5A). We numerically calculate the
transition time as in Fig. S6A. We also calculate the transition time in the
case that catalytic strength of the entities are asymmetric (see Appendix Sec.
S4 in Fig. S6B.

Figure S5: (A) The stochastic reaction dynamics under the SD. Red and Blue in
test tubes (with finite and fixed volume $V$) represent $\rm X_{1}$ and $\rm
X_{2}$-dominant states, respectively. In contrast with the deterministic case,
the fluctuation of concentration causes the transition between states. The
(averaged) transition time from $\rm X_{1}$($\rm X_{2}$)-dominant to $\rm
X_{2}$($\rm X_{2}$)-dominant state is noted as $T_{1\rightarrow
2}$($T_{2\rightarrow 1}$). If $T_{1\rightarrow 2}$ or $T_{2\rightarrow 1}$ are
less than the interval $\varDelta t$, the system has no heredity at all (see
Fig. 1 in the main text). (B) The stochastic dilution protocol. The interval
at each cycle $\varDelta t_{i}$ is randomly distributed, and the dilution
factor $m_{i}$ is $m=e^{\bar{\phi}\langle\varDelta t\rangle}$.  Figure S6: (A)
The transition time from $\rm X_{1}$ to $\rm X_{2}$, or the reverse direction,
$T_{1\rightarrow 2}$ or $T_{2\rightarrow 1}$, varying the dilution rate
$\bar{\phi}$, under SD protocol. The transition time $T_{1\rightarrow 2}$ is
calculated as the average of 10000 trials (the time until $x_{2}>x_{1}$
starting from $(x_{1},x_{2})=(s^{tot},0)$). The lines with different colors
represent the difference in the volume $V$. Solid lines represent the interval
$\varDelta t=1$, and dashed ones $\varDelta t=2$. We set
$\epsilon=0.5,\kappa=8$, in $r(x)x=\epsilon+\kappa x^{2}$. (B) The transition
time between states in the case with asymmetric catalytic strength (from $\rm
X_{1}$-dominant to $\rm X_{2}$, $T_{1\rightarrow 2}$ (blue), and the reverse
direction (green)). We set $\iota=0.25$, or $\iota=0.5$
($\kappa_{1}=\kappa-\frac{\iota}{2}$ and $\kappa_{2}=\kappa+\frac{\iota}{2}$)
in the solid or dashed lines, respectively, and $\varDelta t=1$.

####  B.4.2 stochasticity in the dilution protocols

Not only the chemical reaction dynamics but also the dilution protocol, in
general, could be stochastic. For example, in the compartment growth and
division scenario, the division events could be the stochastic process. The
results are also robust to stochasticity in the dilution protocols; for
example, if the division time $\varDelta t$ is randomly distributed.

Here, we consider the stochastic serial dilution protocol at which the
interval of each dilution event is randomly determined (Fig. S5B). For
example, we consider a case that the interval $\varDelta t$ is randomly
distributed (e.g., a exponential distribution, whose probability density
function $f(\tau)$ is $f(\tau)=\lambda\exp(-\frac{\tau}{\lambda}$) where the
mean $\lambda$ is $\lambda=\langle\varDelta t\rangle$). We consider the ratio
of each dilution $m$ is constant $m=e^{\bar{\phi}\langle\varDelta t\rangle}$
(see Fig. S7A for the time course) 121212The dilution factor $m$ could also be
randomly distributed, or depending on the random variable interval at each
cycle $\varDelta t$. In these cases, the results do not qualitatively differ..
In this case, also, there exists the threshold for the time scale
$\langle\varDelta t\rangle$ for the system with heredity, although the
critical interval $\langle\varDelta t\rangle_{c}^{sd}$ is less than that for
the deterministic case $\varDelta t_{c}^{sd}$ (see Fig. S7B).

Figure S7: (A) Time series of $\rm X_{1}$ and $\rm X_{2}$ in a case with
$r(x)x=\epsilon+\kappa x^{2}$. The interval $\varDelta t$ is randomly
distributed, while $m$ is fixed as $e^{\langle\varDelta t\rangle\phi}$. We set
$\langle\varDelta t\rangle=0.5$, $\epsilon=0.5$ and $\kappa=8$. (B) The
bifurcation diagram of the averaged value of $x_{1}$ and $x_{2}$ just before
the dilutions in the stationary trajectory.

###  B.5 Differential reproduction of different chemical compositional states

Finally, we investigate situations where the dilution rates depend on the
chemical composition, and the system has different growth rates between the
two states. As noted in the introduction, this is, in fact, a crucial property
necessary for a population of compartmentalized chemical reaction systems to
undergo Darwinian evolution (see Fig. S8).

In such a case, the dilution rate $\hat{\phi}\big{(}t,\bm{x}(t)\big{)}$ is
given by $\hat{\phi}\big{(}t,\bm{x}(t)\big{)}\equiv{\frac{dV}{dt}}/{V},$ such
that

$\int_{0}^{\varDelta
t}\hat{\phi}\big{(}t,\bm{x}^{*}_{(i)}(t)\big{)}=\bar{\phi}_{i}\varDelta t,$
(S36)

where $i=1,2$ and $\bm{x}^{*}_{(i)}(t)$ represents the steady-state
trajectories of $\bm{x}$, where $x_{i}$ is dominant. Here, differential
reproduction implies $\bar{\phi}_{1}\neq\bar{\phi}_{2}$. Recall that
$m=e^{\bar{\phi}\varDelta t}$ is the number of offspring at one generation
with interval $\varDelta t$ (main text Sec. A).

As a simple example, we consider the dilution due to the growth of
compartments which depends on their components,

$\hat{\phi}\big{(}t,\bm{x}(t)\big{)}=\phi(t)-\gamma\delta,$ (S37)

where $\gamma$ is a constant (control parameter), $\delta=x_{1}-x_{2}$, and
${\phi}(t)$ is the dilution rate independent with the component; here, we
choose Eq. 11 in Material and Methods as ${\phi}(t)$, such that
$\int_{0}^{\varDelta t}\phi(t)dt=\bar{\phi}\varDelta t$. Naturally, in this
case, the growth rate of the two states are different,
$\bar{\phi}_{1}<\bar{\phi}_{2}$, if $\gamma>0$. Notably, this dependency of
$\hat{\phi}$ on the composition does not change the critical long-term
dilution rate at which the system loses bistability, $\phi_{c}$ (Fig. S10A),
i.e., $\phi_{c}$ does not depend on $\gamma$. Generally, the bistability of
the system does not seem to be affected by the asymmetric factor of
$\hat{\phi}$. We can see this by considering the autocatalytic system is under
arbitrary protocols with the dilution rate depending on the chemical
composition,

$\frac{dx_{i}}{dt}=sr(x_{i})x_{i}-\phi(t,\chi,\delta).$ (S38)

where $\chi=x_{1}+x_{2},\delta=x_{1}-x_{2}$. The dependency of $\phi$ on
$\delta$ is interpreted as differential reproduction between the states. The
same calculation as in Appendix Sec. A leads to

$\frac{d\chi}{dt}=sr(\frac{\chi}{2})\chi-\phi(t,\chi,0)\chi+\mathcal{O}(\delta),\quad\frac{d\delta}{dt}=s\left(r(\frac{\chi}{2})+\frac{1}{2}\chi\frac{dr}{dx}(\frac{\chi}{2})\right)\delta-\phi(t,\chi,0)\delta+\mathcal{O}(\delta^{2}),$
(S39)

where $\chi\gg\delta$ is assumed. Then, the deviation of $\delta$ in one cycle
is calculated as the same as in Appendix Sec. A,

$\displaystyle\log\left|\frac{\delta(t)}{\chi(t)}\right|$
$\displaystyle=\int_{0}^{t}\frac{1}{2}s\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}\chi
dt$ $\displaystyle+$
$\displaystyle\log\left|\frac{\delta(0)}{\chi(0)}\right|,$ (S40)

which, notably, does not depend on the asymmetric part of
$\phi(t,\chi,\delta)$.

Next, we consider a case where $\phi$ depends on symmetrically $x_{1}$ and
$x_{2}$, i.e.,
$\hat{\phi}\big{(}t,x_{1},x_{2}\big{)}=\hat{\phi}(t,x_{2},x_{1})$:

$\hat{\phi}\big{(}t,\bm{x}(t)\big{)}=\phi(t)-\gamma\chi,$ (S41)

where $\chi=x_{1}+x_{2}$, while the catalytic strength of two entities,
$\kappa_{1}$ and $\kappa_{2}$, are asymmetrical, i.e.,
$\kappa_{1}\neq\kappa_{2}$, as discussed in Appendix Sec. B.2. In this case,
also the growth rate is different for different states (Fig. S10).

Figure S8: Schematic diagram for the differential reproduction of
compartments. The compartments with the $\rm{X}_{1}$\- and
$\rm{X}_{2}$-dominant states are represented as red and blue circles,
respectively. Each compartment grows and produces
$m_{1}=e^{\bar{\phi}_{1}\varDelta t}$ (or $m_{2}=e^{\bar{\phi}_{2}\varDelta
t}$) daughters with the equal volumes in the interval $\varDelta t$. In this
figure, the $\rm X_{2}$-dominant compartment grows faster and produces more
offspring, i.e., $m_{2}>m_{1}$. Figure S9: Schematic diagrams of competing
autocatalytic entities encapsulated by a growing compartment. The grey arrow
represents the growth of the compartment, which is promoted (dotted arrow) or
inhibited (bar-headed dotted arrow) by entities. (A) The entities have
symmetric catalytic strength (the same model as Fig. S2A), but affect the
compartment growth asymmetrically (see Eq. S14). (B) The entities have
asymmetric catalytic strength (the same model as Fig. S2C), but affect the
growth symmetrically (see Eq. S18). Figure S10: (A) The bifurcation diagram
with a varying cycle interval $\varDelta t$ in a case with
$r(x)x=\epsilon+\kappa x^{2}$, the asymmetric dilution rate
$\hat{\phi}(t,\bm{x}(t))=\phi(t)-\gamma\delta$, where $\phi(t)$ is given by
Eq. 11. The grey curve represents the concentration, $x_{1}$ (or $x_{2}$), at
just before dilution in the stationary trajectory when $\gamma=0$, and solid
and dashed colored curves represent the concentration when $\gamma=1$ in $\rm
X_{1}$-dominant and $\rm X_{2}$-dominant states, respectively; the red and
blue curves represents $x_{1}$ and $x_{2}$, respectively. We set the
parameters as $\kappa=8,\epsilon=0.5$ and $\alpha=2$. (B) The reproduction
rates of the compartment (i.e., the long-term cycle averaged dilution rate),
$\frac{1}{\varDelta t}\int_{T}^{T+\varDelta t}\hat{\phi}(t,\bm{x}(t))dt$ as
$T\rightarrow\infty$, at each compositional state (stationary trajectory
$\\{\bm{x}(t)\\}$), varying the interval $\varDelta t$. The dotted line
represents the critical point $\varDelta t_{c}^{sd}$ which divides the regions
where the system has bistability or not. (C) The same bifurcation diagram as
(A) in a case with $r_{i}(x)x=\epsilon+\kappa_{i}x^{2}$ and the symmetric
dilution rate $\hat{\phi}(t,\bm{x}(t))=\phi(t)+\gamma\chi$. (D) The
reproduction rates of the compartment, under the same setup as (B). We set the
parameters as $\iota=-1$, $\kappa=8,\epsilon=0.5$ and $\alpha=2$.

We also consider a situation where the division interval $\varDelta t$ depends
on the composition $\bm{x}$. For example, we decide whether $\varDelta t$ is
equal to $\varDelta t_{1}$ if $x_{1}>x_{2}$ at the beginning of cycle, or
$\varDelta t_{2}$ otherwise. Trivially, in this case also, the system has the
bistability if $\varDelta t_{1}$ and $\varDelta t_{2}$ are below the critical
value discussed in the main text Sec. C, $\varDelta t_{1},\varDelta
t_{2}<\varDelta t_{c}$.

Overall, even the case that growth rate $\phi$ or/and the division interval
$\varDelta t$ depend on the chemical composition of the system, we can
guarantee the bistability if $\phi$ at each stationary trajectory,
$\bm{x}_{(1)}^{*}(t)$ and $\bm{x}_{(2)}^{*}(t)$ satisfies
$\bar{\phi}_{1},\bar{\phi}_{2}<\phi_{c}^{sd}$ and $\varDelta t_{1},\varDelta
t_{2}<\varDelta t_{c}^{sd}$.

##  Appendix C ACS based on _Azoarcus_ ribozyme coupled with metabolism

In the case of the system based on _Azoarcus_ ribozyme coupled with
catabolism/anabolism reaction (Eq. 6 in Sec. II.F in the main text), we drew
the same figures as the model (Eq. 1) in Sec. II.B, C and D. These results are
qualitatively similar as shown in Fig. S11 and Fig. S12.

Figure S11: The red and blue curves represent the nullclines for $\rm{X}_{1}$
and $\rm{X}_{2}$, respectively, just before dilution. We set $\varDelta
t=1.5$, $3$, $\kappa=2.5$, $\epsilon=0.25$, $\bar{\phi}=1$, and $b=0.1$.
Figure S12: (A) The bifurcation diagram for the concentration at just before
dilution varying the cycle interval $\varDelta t$. The lines with different
colors represent the difference in the background reaction rate $\epsilon$.
(B) The lines represent the critical point $\varDelta t_{c}$, which divides
the regions where the system has bistability or not. We set the parameters as
$\kappa=2.5$ and $\bar{\phi}=1$. Figure S13: The parameter region for the
system based on _Azoarcus_ ribozyme (Eq. 6) with the bistability under SD
protocol. The same plot as Fig. 6D in the main text, but its horizontal axis
is the dilution factor $m$, $m=\exp(\bar{\phi}\varDelta t)$, instead of
$\bar{\phi}$. The lines with different colors represent the boundary between
with/without bistability under different backward reaction rates $b$. We set
$\epsilon=0.25,\kappa=2.5$, and $b=0.1$.

##  Appendix D Robustness of the results for the system based on _Azoarcus_
ribozyme

###  D.1 In the case with asymmetric catalytic efficiency

A similar relation also appears even when the catalytic activity of two
species are different, i.e., $\kappa_{1}\neq\kappa_{2}$ (Fig. S14A). Although
the bifurcation at which the bistability disappears is discontinuous, it is at
the similar $\varDelta t_{c}$ unless the difference between $\kappa_{1}$ and
$\kappa_{2}$ is not so large Fig. S14B.

Figure S14: (A) The concentrations of the species $\rm{X}_{1}$ and
$\rm{X}_{2}$ with a varying the period of the dilution cycles $\varDelta t$,
when the catalytic strengths are asymmetric, $\kappa_{1}\neq\kappa_{2}$. We
set $\kappa_{1}=\kappa-\frac{\iota}{2}$ and
$\kappa_{2}=\kappa+\frac{\iota}{2}$, and $\kappa_{1}-\kappa_{2}=\iota=0.1$. In
contrast with the symmetric case (i.e., $\iota=0$), the $\rm{X}_{2}$-dominant
state (and thus the bistability) disappears discontinuously at around
$\varDelta t\sim 1.71$. We set the other parameters as
$\kappa=2.5,\epsilon=0.25$, $b=0.1$ and $\bar{\phi}=1$. (B) The dependence of
the critical interval $\varDelta t_{c}^{sd}$ of the cycle on the dilution rate
$\bar{\phi}$ when the catalytic strengths are asymmetric:
$\kappa_{1}-\kappa_{2}=\iota=0,0.05,0.1$, and $0.2$. We set the other
parameters as $\kappa=2.5,\epsilon=0.25$, and $\bar{\phi}=1$.

###  D.2 The variation of two-step ACSs with or without heredity

Here, assuming the mass action kinetics, we investigate alternative models
with two-step reactions, similar to the model discussed in the main text (Fig.
S15A).

Figure S15: Schematics of the alternative models for the symmetric competing
autocatalytic entities with two reaction steps. Solid arrows represent the
reactions, and dotted arrows represent the catalysis. (A) The autocatalytic
ribozymes coupled with anabolism and catabolism reactions in the main text
Fig. 6B. (B) A similar reaction system with two catalyzed reaction steps, but
the substrate S is consumed to produce $\rm X^{\prime}_{1}$ and $\rm
X^{\prime}_{2}$ instead of $\rm X_{1}$ and $\rm X_{2}$. (C) Two competing
Lotka-Volterra model.

Firstly, if we assume the substrate $\rm{S}$ is consumed in another reaction
(see Fig. S15B), the result does not change qualitatively; Here, assuming the
CSTR condition, we consider the modified rate equations,

$\frac{dx^{\prime}_{i}}{dt}=(\epsilon+\kappa
x_{i})(s-x^{\prime}_{i})-\bar{\phi}x^{\prime}_{i},\quad\frac{dx_{i}}{dt}=(\epsilon+\kappa
x_{i})x^{\prime}_{i}-\bar{\phi}x_{i},$ (S42)

where $s^{tot}=s+x^{\prime}_{1}+x^{\prime}_{2}+x_{1}+x_{2}$.
$\rm{X}^{\prime}_{1}$ and $\rm{X}^{\prime}_{2}$ are converted from the shared
substrate S, and which are further converted into $\rm{X}_{1}$ and
$\rm{X}_{2}$, respectively.

Then, we consider the dynamics near the steady state, and if we assume
$x^{\prime}_{1}$ and $x^{\prime}_{2}$ can be adiabatically eliminated from
$\frac{dx^{\prime}_{i}}{dt}=0$,

$\frac{dx_{i}}{dt}=\frac{s(\epsilon+\kappa x_{i})^{2}}{\epsilon+\kappa
x_{i}+\bar{\phi}}-\bar{\phi}x_{i}.$ (S43)

Thus, this rate equation corresponds to Eq. 1 with the reproduction rate
function $xr(x)=\frac{(\epsilon+\kappa x)^{2}}{\epsilon+\kappa x+\bar{\phi}}$.
Further, in this model also $\rm{X}_{1}$ and $\rm{X}_{2}$ compete for the same
substrate $\rm{S}$; therefore, this model also shows the bistability, if the
condition $\frac{d}{dt}r(\chi^{*}/2)>0$ is satisfied, where
$s^{*}r(\chi^{*}/2)=\phi$. (Later, the condition for this reaction system to
be bistable under the serial dilution protocol is derived.)

Secondly, for the system show bistability, both of the reactions, from
$\rm{S}$ to $X^{\prime}_{1}$ and from $X^{\prime}_{1}$ to $\rm{X}_{1}$ have to
be catalyzed by $\rm{X}_{1}$; for example, we modify the model as the reaction
$\rm{S}$ to $X^{\prime}_{1}$ is catalyzed by $X^{\prime}_{1}$ instead of
$\rm{X}_{1}$ (Fig. S15C). Under the CSTR condition, the rate equations of the
model are

$\frac{dx^{\prime}_{i}}{dt}=(\epsilon+\kappa
x^{\prime}_{i})(s-x^{\prime}_{i})-\bar{\phi}x^{\prime}_{i},\quad\frac{dx_{i}}{dt}=(\epsilon+\kappa
x_{i})x^{\prime}_{i}-\bar{\phi}x_{i},$ (S44)

where $s^{tot}=s+x^{\prime}_{1}+x^{\prime}_{2}+x_{1}+x_{2}$. This system has
only one stable fixed point and does not show bistability. This is because,
$\rm{X}_{1}$ and $\rm{X}_{2}$ does not compete for the same resource for their
replications but $\rm{X}^{\prime}_{1}$ and $\rm{X}^{\prime}_{2}$ does. Then,
the effective reproduction rate functions for $X^{\prime}_{1}$ and
$X^{\prime}_{2}$ do not satisfy the condition for the symmetric state to be
unstable.

In conclusion, if the mass action kinetics is assumed, for the competing
autocatalytic chemical reaction networks sharing the same substrate to have
bistability (i.e., the nonlinear reproduction rate function), it is required
at least two reaction steps to produce the autocatalytic entities (catalysts),
which is catalyzed by the entities themselves.

###  D.3 derivation of $\varDelta t_{c}^{sd}$ in a case with two-step
catalyzed reactions

We further consider the two competing chemical reaction networks with two
catalyzed reaction steps in Fig. S15B, under the SD protocol:

$\frac{dx^{\prime}_{i}}{dt}=(s(\\{x^{\prime}_{j}\\},\\{x_{j}\\},t)-x^{\prime}_{i})r(x_{i}),\quad\frac{dx_{i}}{dt}=x^{\prime}_{i}r(x_{i}),$
(S45)

where $r(x)=\epsilon+\kappa x$. Here, we define

$\begin{split}\chi=x_{1}+x_{2},\quad\chi^{\prime}=x_{1}^{\prime}+x_{2}^{\prime},\quad\delta=x_{1}-x_{2},\quad\delta^{\prime}=x_{1}^{\prime}-x_{2}^{\prime}.\end{split}$
(S46)

The time derivative of the above is derived as

$\displaystyle\frac{d\delta}{dt}$
$\displaystyle=\frac{1}{2}\delta\chi^{\prime}\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}+\delta^{\prime}r\Bigl{(}\frac{\chi}{2}\Bigr{)}$
$\displaystyle+$
$\displaystyle\mathcal{O}(\delta^{2},\delta\delta^{\prime},\delta^{\prime
2}),$ (S47) $\displaystyle\frac{d\chi}{dt}$
$\displaystyle=\chi^{\prime}r\Bigl{(}\frac{\chi}{2}\Bigr{)}$ $\displaystyle+$
$\displaystyle\mathcal{O}(\delta^{2},\delta\delta^{\prime},\delta^{\prime
2}),$ $\displaystyle\frac{d\delta^{\prime}}{dt}$
$\displaystyle=s\delta\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}-\delta^{\prime}r\Bigl{(}\frac{\chi}{2}\Bigr{)}-\frac{1}{2}\chi^{\prime}\delta\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}$
$\displaystyle+$
$\displaystyle\mathcal{O}(\delta^{2},\delta\delta^{\prime},\delta^{\prime
2}),$ $\displaystyle\frac{d\chi^{\prime}}{dt}$
$\displaystyle=2sr\Bigl{(}\frac{\chi}{2}\Bigr{)}-\chi^{\prime}r\Bigl{(}\frac{\chi}{2}\Bigr{)}$
$\displaystyle+$
$\displaystyle\mathcal{O}(\delta^{2},\delta\delta^{\prime},\delta^{\prime
2}),$

where we used
$r(x_{i})=r(\frac{\chi}{2})\pm\frac{\delta}{2}\frac{dr}{dx}(\frac{\chi}{2})+\mathcal{O}(\delta^{2})$.

To determine the deviation of $\delta$ in one cycle, we integrate
$\frac{d\delta}{dt}/\delta$,

$\displaystyle\int_{0}^{t}\frac{\frac{d\delta}{dt}}{\delta}dt$
$\displaystyle=\frac{1}{2}\int_{0}^{t}\chi^{\prime}\frac{dr}{dx}\Bigl{(}\frac{\chi}{2}\Bigr{)}dt$
$\displaystyle+$
$\displaystyle\int_{0}^{t}r\Bigl{(}\frac{\chi}{2}\Bigr{)}\frac{\delta^{\prime}}{\delta}dt,$
(S48)
$\displaystyle=\int_{\frac{\chi(0)}{2}}^{\frac{\chi(t)}{2}}\frac{\frac{dr}{dx}}{r}dx$
$\displaystyle+$
$\displaystyle\int_{\frac{\chi(0)}{2}}^{\frac{\chi(t)}{2}}\frac{2\delta^{\prime}}{\chi^{\prime}\delta}dx,$

where we used the change of the variable
$\frac{dt}{d\chi}=1/(\chi^{\prime}r(\frac{\chi}{2}))$.

Here, to calculate Eq. S7, we have to estimate $\delta^{\prime}/\chi^{\prime}$
in the second integral in Eq. S7. If we assume in the second term, in the
dominant part of the integrate, $\delta^{\prime}/\chi^{\prime}$ has the
scaling relation $\delta^{\prime}/\chi^{\prime}\sim\mathcal{O}(1)$,

$\frac{\delta^{\prime}}{\chi^{\prime}}\approx\frac{\frac{d\delta^{\prime}}{dt}}{\frac{d\chi^{\prime}}{dt}}=\frac{1}{2}\frac{\frac{dr}{dx}}{r}\delta-\frac{\delta^{\prime}}{2s-\chi^{\prime}}\approx\frac{1}{2}\frac{\frac{dr}{dx}}{r}\delta,$
(S49)

where we used $\frac{d}{dt}(\frac{\delta^{\prime}}{\chi^{\prime}})=0$ in the
first approximation, and in the second approximation, we assumed
$s\gg\chi^{\prime}$. We substitute this to the above equation,

$\log\left|\frac{\delta(t)}{\delta(0)}\right|\approx
2\log\left|\frac{r(\frac{\chi(t)}{2})}{r(\frac{\chi(0)}{2})}\right|.$ (S50)

The threshold of $\varDelta t$ for the symmetric (i.e., $\delta=0$) trajectory
to be unstable is,

$\varDelta t_{c}^{sd}=\frac{2}{\bar{\phi}}\log\left(1+\frac{\kappa
s^{tot}}{2\epsilon}\right)-C_{0},$ (S51)

where $C_{0}$ is a constant value, which is determined by numerically as
$C_{0}\sim 3.6$ (see Fig. S16A).

While we use the asymptotic relation $\frac{\delta^{\prime}}{\chi^{\prime
2}}\sim\mathcal{O}(1)$ as $t$ becomes large. Then,

$\frac{\delta^{\prime}}{\chi^{\prime}}\approx\frac{1}{2}\frac{\frac{d\delta^{\prime}}{dt}}{\frac{d\chi^{\prime}}{dt}}\approx\frac{1}{4}\frac{\frac{dr}{dx}}{r}\delta,$
(S52)

where $\frac{d}{dt}(\frac{\delta^{\prime}}{\chi^{\prime 2}})=0$ in the first
approximation, and $s\gg\chi^{\prime}$ in the second as like before. Thus,

$\log\left|\frac{\delta(t)}{\delta(0)}\right|\approx\frac{3}{2}\log\left|\frac{r(\frac{\chi(t)}{2})}{r(\frac{\chi(0)}{2})}\right|,$
(S53)

and then

$\varDelta t_{c}^{sd}=\frac{3}{2\phi}\log\left(1+\frac{\kappa
s^{tot}}{2\epsilon}\right)-C_{1},$ (S54)

where $C_{1}$ is numerically determined as $C_{1}\sim 0.477$ (see Fig. S16B).

Figure S16: (A) The solid line represents the critical point $\varDelta
t_{c}$, which divides the regions where the system has bistability or not. The
solid line represents the theoretical line for $\varDelta t_{c}$ determined by
the relation in Eq. S51. (B) The interval threshold $\varDelta t_{c}$ vs the
catalytic efficiency $\kappa$, fitted by
$\frac{3}{2}\log(1+\kappa/2/\epsilon)-C_{1}$. We set
$\bar{\phi}=\bar{\sigma}=1$.