Abstract

One proposed scenario for the emergence of biochemical oscillations is that they may have provided the basic mechanism behind cellular self-replication by growth and division. However, alternative scenarios not requiring any chemical oscillation have also been proposed. Each of the various protocell models proposed to support one or another scenario comes with its own set of specific assumptions, which makes it difficult to ascertain whether chemical oscillations are required or not for cellular self-replication. This article compares these two cases within a single whole-cell model framework. This model relies upon a membrane embedding a chemical reaction network (CRN) synthesizing all the cellular constituents, including the membrane, by feeding from an external nutrient. Assuming the osmolarity is kept constant, the system dynamics are governed by a set of nonlinear differential equations coupling the chemical concentrations and the surface-area-to-volume ratio. The resulting asymptotic trajectories are used to determine the cellular shape by minimizing the membrane bending energy (within an approximate predefined family of shapes). While the stationary case can be handled quite generally, the oscillatory one is investigated using a simple oscillating CRN example, which is used to identify features that are expected to hold for any network. It is found that cellular self-replication can be reached with or without chemical oscillations, and that a requirement common to both stationary and oscillatory cases is that a minimum spontaneous curvature of the membrane is required for the cell to divide once its area and volume are both doubled. The oscillatory case can result in a greater variety of cellular shape trajectories but raises additional constraints for cellular division and self-replication: (i) the ratio of doubling time to oscillation period should be an integer, and (ii) if the oscillation amplitude is sufficiently high, then the spontaneous curvature must be below a maximum value to avoid early division before the end of the cycle. Because of these additional stringent constraints, it is likely that early protocells did not rely upon chemical oscillations. Biochemical oscillations typical of modern evolved cells may have emerged later through evolution for other reasons (e.g., metabolic advantage) and must have required additional feedback mechanisms for such a self-replicating system to be robust against even slight environmental variations (e.g., temperature fluctuations).

1 Introduction

Biochemical oscillations are ubiquitous in living cells. They have been observed in various contexts ranging from metabolism (e.g., glycolytic oscillations, cyclic AMP oscillations) to specific proteins such as cyclins (see [27] for a review). An open question is when such oscillations emerged: Were they already present in early protocells? Or did they appear later in evolution? A related question is that of the benefits and costs brought by such oscillations for a self-replicating biochemical system.

One tentative explanation for the emergence of biochemical oscillations is that they might have provided the basic mechanism behind cellular division, as investigated in [37, 38, 23]. Oscillations in the concentration of biochemical species inside the cell should also result in oscillations in the total cell osmolarity, unless its shape (characterized by the surface-area-to-volume ratio) also oscillates so as to keep the osmolarity constant and balanced with that of the extracellular medium. Such oscillations have been numerically observed by Surovstev et al. [37, 38] in a whole-cell model relying upon a membrane embedding an oscillating chemical reaction network (CRN) synthesizing all cell constituents including the membrane itself. In [37], the cell was assumed to take the shape of two overlapping spheres, the size and offset of which were uniquely determined from the concentrations and surface-area-to-volume ratio trajectories. The oscillations resulted in cellular division and self-replication, provided the oscillation period matched the doubling time. In [38], the cell shape determination was refined by minimizing the membrane bending energy, and the inclusion of a contractile ring was found necessary to result in cell constriction and division. In a similar line of thought, Martin et al. [23] have also considered several examples of oscillating CRNs embedded in a growing protocell, but volume-changing effects were taken into account in a simplified phenomenological manner and the issue of the cellular shape was not tackled.

However, cellular growth and division has also been observed in other numerical models, without requiring chemical oscillations. Macía and Solé [22] have developed a protocell model in which the cell volume grows under the influence of a single chemical reaction transforming a nutrient into a cellular constituent, and in which osmotic pressure and dynamic water flow across the membrane were taken into account. Cell growth, constriction, and division were numerically observed, but this model assumed that the single chemical reaction was catalyzed by enzymatic centers arbitrarily positioned at opposite ends of the spherical cell at the beginning of the cell cycle. Besides, lipids were assumed to be in constant supply in the inside and outside media and were not synthesized by the chemical reaction. This matches certain protocell experimental schemes [11, 1], but is different from the generic minimum working cell considered by Surovstev et al. [37, 38]. Mavelli et al. [24] and Shirt-Ediss et al. [36, 35] have also proposed protocell models relying upon simple non-oscillating CRNs, for which the membrane precursor is synthesized by the embedded reaction network or not, and in which osmotic pressure is quasi-instantaneously balanced across the membrane. But how the cellular shape varies across the cell cycle, and how this may lead to cellular division or not, were not considered. From a different perspective, other authors [9, 10, 16] have focused on vesicle growth and possible self-replication, assuming a phenomenological membrane surface area growth, and using a rigorous mathematical determination of vesicle shape by minimizing the membrane bending energy within the framework of the spontaneous-curvature model. Self-replication was found possible, within a certain range of membrane parameters. However, such approaches did not consider any inner cellular metabolism.

It is thus yet unclear whether chemical oscillations are required for the emergence of cellular self-replication by growth and division. If the answer proved to be yes, a consequence would be that such oscillations were most likely already present in early protocells. If the answer proved to be no, the question of benefits and costs brought by such oscillations should also be addressed.

Clarifying these points is the purpose of the present work. The relation between chemical oscillations and cellular division is investigated by comparing the two cases, with and without chemical oscillations, within a generic whole-cell modeling framework similar to that of [37], [38], or [7]. The whole-cell model consists of a membrane embedding an inner metabolism, which synthesizes not only all cytoplasmic constituents but also the membrane precursor. The balance of the osmotic pressure across the membrane results in an ordinary differential equation (ODE) system coupling the concentrations and the surface-area-to-volume ratio. The resulting asymptotic trajectories are used to determine the cell shape trajectory by minimizing the membrane bending energy within the spontaneous-curvature framework, simplifying the problem by restricting this minimization to an approximate predefined family of shapes.

Using CRNs with random topology and kinetics, previously observed asymptotic trajectories were typically stationary without any oscillatory behavior, which corresponds for instance to filamentous growth [7]. The present work shows that such a constant surface-area-to-volume ratio may correspond not only to filamentous growth, but also to different cell shape trajectories, resulting in cellular growth and division provided the membrane spontaneous curvature (which depends on the membrane structural properties) is above a minimum value (which depends on the embedded chemistry). This shows that biochemical oscillations are not required for cellular self-replication by growth and division.

Besides, examples of asymptotic oscillating trajectories are generated using a simple oscillating CRN example, which results in a variety of shape trajectories. This also results in cellular growth and division provided certain constraints are met: As in the stationary case, the spontaneous curvature should be above a minimum value for the cell to divide once the surface area and volume are doubled. Besides, the ratio of cellular doubling time to oscillation period should be an integer (which generalizes the synchronization condition first introduced in [38]); and, depending on the amplitude of the oscillations, the spontaneous curvature must be below a maximum value to avoid early division before the end of the cycle. This shows that, though not required, biochemical oscillations may also be compatible with cellular self-replication by growth and division, but that this comes at the expense of stringent constraints.

The former synchronization constraint requires a tight adjustment of chemical kinetic parameters, which should be highly sensitive to even slight environmental variations such as a mere temperature change. Additional feedback mechanisms (beyond the simplistic features of the considered oscillating CRN example) are required for such a system to be robust. It is thus likely that self-replication of early protocells did not rely upon chemical oscillations, because the oscillating case appears less robust than the stationary one. Biochemical oscillations may thus have emerged later in the evolution together with robust feedback mechanisms ensuring the system stability, because such oscillations may have presented other advantages (e.g., the metabolic advantage of producing more ATP per unit consumed glucose when going though the citric acid cycle rather than only through glycolysis).

This article is organized as follows. Section 2 briefly recalls the cell model and the resulting ordinary differential equation (ODE) system, and restates the synchronization constraint that must be met by any self-replicating systems: that the ratio of doubling time to oscillation period should be an integer. Section 3 presents the simplified family of shapes used for the minimization of the membrane bending energy, and explains how the cellular shape trajectory is deduced from the asymptotic trajectories, for concentrations and for the surface-area-to-volume ratio. Section 4 considers first the case of stationary trajectories, and second the case of oscillating trajectories using a simple example oscillating CRN. Section 5 is devoted to a discussion. Section 6 gives a conclusion.

2 Cell Model

2.1 Ordinary Differential Equation (ODE) System

The model presented in this section and illustrated in Figure 1 is similar to that of [37], [38], or [7]. It assumes that all concentrations are spatially homogeneous inside the cell and in the outside growth medium, and that the nutrient concentration cA remains constant. The nutrient A is imported into the cell by a surfacic process characterized by a reaction vector Fin. The membrane precursor is incorporated into a self-assembled membrane by a surfacic process characterized by a reaction vector Fout.

Figure 1. 

Schematic of the cell model. The nutrient A is in constant supply in the outside growth medium, is imported into the cell with a rate Fin per unit area, and is transformed into all other cell constituents including the membrane precursor. This precursor is incorporated into the growing membrane with a rate Fout per unit area.

Figure 1. 

Schematic of the cell model. The nutrient A is in constant supply in the outside growth medium, is imported into the cell with a rate Fin per unit area, and is transformed into all other cell constituents including the membrane precursor. This precursor is incorporated into the growing membrane with a rate Fout per unit area.

Concentration trajectories satisfy the following set of differential equations:
c˙=Sf+ρFinFoutλc
(1)
where c is the concentration vector inside the cell, S is the CRN stoichiometry matrix, f = f(c) is the reaction rate vector, ρ = A/V is the surface area-to-volume ratio, and λ = V˙/V is the instantaneous growth rate (with A and V being the cell surface area and volume, respectively). This system is underdetermined because it consists in N equations (where N is the number of chemical species inside the cell) for N + 2 variables (N concentrations, plus ρ and λ). Additional relations making the system fully determined are derived next:
  • 1. 
    Water is assumed to flow sufficiently fast across the membrane to balance quasi-instantaneously the osmotic pressure across the membrane. This assumption is consistent with its much higher permeability to water than to other solutes, as measured across various membranes [29]. This constant osmolarity may be expressed as uTc = cA, where u is the osmolarity vector such that ui is the number of particles in solution when dissolving one molecule of the chemical species indexed by i. The inside osmolarity is thus constant, uTċ = 0. Multiplying each side of Equation 1 by uT and rearranging terms, gives
    λ=uTSf+ρFinFoutuTc.
    (2)
  • 2. 
    An additional differential equation may be obtained for the time evolution of the surface-area-to-volume ratio ρ. By definition we have ρ = A/V, where A is the membrane surface area and V the cell volume. Its logarithmic derivative is given by ρ˙/ρ=A˙/AV˙/V=A˙/Aλ. Here A increases as the result of incorporation of the precursor into the growing membrane. The number of moles of precursor being incorporated per unit area and per unit time is given by Fout,1 which results in a relative rate of increase for A given by
    A˙A=FoutNme
    (3)
    where Nme is the number of molecules per unit area in the self-assembled membrane. This results in the following expression for the logarithmic derivative of ρ:
    ρ˙ρ=FoutNmeλ.
    (4)

Equations 1, 2, and 4 give an autonomous ODE system governing the time evolution of concentrations as well as that of the surface-area-to-volume ratio.

2.2 Doubling Time Must Be a Multiple of the Oscillation Period

If an oscillatory regime is reached, an important constraint arises for such a mechanism to be compatible with cellular division and self-replication. Self-replication requires that both volume and surface area be doubled in time Tdbl, that is, A(t + Tdbl) = 2A(t) and V(t + Tdbl) = 2V(t), so that at the end of one cycle, one ends up with two cells each identical to the original cell. The average growth rate is λavg = log(2)/Tdbl. Defining ζA(t) = exp(−λavgt)A(t) and ζV(t) = exp(−λavgt)V(t), we see that ζA(t) and ζV(t) are periodic with period Tdbl.2 Then ρt=A/V=ζA(t)/ζV(t) is also periodic, because it is the ratio of two periodic functions. Besides, self-replication requires that concentrations resume their original values after one cell cycle, that is, that concentrations also be periodic with period Tdbl.

It should be noted that the doubling time Tdbl is the duration over which the integral of A˙/A (as given by Equation 3) reaches log(2), and that the trajectories c(t), ρ(t) are determined by integration of the ODE system defined by Equations 1, 2, and 4. There is no a priori reason why the resulting oscillation period Tperiod should match Tdbl. But the above reasoning shows that such a match is a requirement for a self-replicating system.

To be more precise, the shortest period Tperiod common to ζA(t), ζV(t), ρ(t), and concentrations need not necessarily be Tdbl, but can also be one of its submultiples, that is, Tdbl = nTperiod, where n is an integer. The simplest case where n = 1 corresponds to the concentrations and surface-area-to-volume ratio resuming their original values in exactly one cell cycle, whereas n > 1 corresponds to cases where they display n oscillations (and resume their original values) in one cell cycle time. This a mere generalization of the synchronization condition first introduced by Surovstev et al. [38], in which only the case n = 1 was considered.

It should be noted that this requirement is not specific to the particular example described in this work, nor to any specific cellular model, but should hold for any self-replication model where both surface area and volume are doubled in one cycle time, while concentrations resume their original values after one cycle time.3

3 From Asymptotic Trajectories to Cellular Shape

3.1 Spontaneous-Curvature Model, with Predefined Family of Shapes

There has been an extensive body of work on the determination of a vesicle shape for given volume V and surface area A using the spontaneous-curvature framework (see [14, 34, 39, 33, 40], and references therein). Within this framework, the shape is the one that minimizes the membrane bending energy given by
Eb=κ2C1+C2C02dA+κGC1C2dA
(5)
where κ and κG are the bending rigidity and the Gaussian bending rigidity, respectively, C1 and C2 are the principal curvatures, and C0 is the spontaneous curvature. The second term is constant for topologically equivalent surfaces, and is thus omitted in the minimization problem. The exact resolution of this minimization problem entails (i) using a Lagrange multiplier approach to identify shapes with extremal energies as solutions of a third-degree differential equation (assuming axisymmetric shapes), (ii) determining which are local minima versus saddle points, and (iii) carrying out a stability analysis [14, 34, 39]. This approach has successfully explained the shape of erythrocytes and of various artificial vesicles [33, 40]. These results may be summarized in phase diagrams. However, such phase diagrams have not been fully explored yet.4

The present work relies upon a simpler approximate approach, comparable to that of Canham [12]. The shape is assumed to belong to a predefined family such as shown in Figure 2, and is parametrized by only three parameters: cell radius R, cylindrical length L, and constriction angle ω. Such a family may represent a single sphere (L = 0, ω = 0), an elongated rod (L > 0, ω = 0), an elongated rod with constriction L>0,0<ω<π2, two overlapping spheres L=0,0<ω<π2, or a double sphere L=0,ω=π2.5

Figure 2. 

Predefined family of shapes used as an approximation for the minimization of the bending energy. The shape is parametrized by a radius R, a cylindrical length L, and a constriction angle ω. In the most general case shown here 0<ω<π2, the cell is an elongated rod with constriction in the middle. The cell divides into two daughter cells when ω=π2 (two elongated rods connected by an infinitely narrow neck).

Figure 2. 

Predefined family of shapes used as an approximation for the minimization of the bending energy. The shape is parametrized by a radius R, a cylindrical length L, and a constriction angle ω. In the most general case shown here 0<ω<π2, the cell is an elongated rod with constriction in the middle. The cell divides into two daughter cells when ω=π2 (two elongated rods connected by an infinitely narrow neck).

Starting from an assumed spherical shape (with radius Rt0=3ρt0), the bending energy is computed for all compatible shapes within the predefined family, and the shape is determined as that minimizing the bending energy. This minimization is repeated over the various time points of the trajectories, to give the cell shape and its variation across the cell cycle.6

3.2 Minimization of the Bending Energy within the Predefined Family of Shapes

Instead of using the volume V and the surface area A, it proves convenient to use the reduced volume v and the normalized area a. Here v is defined as V divided by the volume of the sphere enclosed by A, that is,
v=6πVA32=6πρA.
(6)
It satisfies 0 ≤ v ≤ 1 with v = 1 for a sphere. a is defined by
a=AAt0.
(7)
Thus v varies from 1 at the beginning of the cycle (because of the assumed spherical shape) to 12 at the end of the cycle (because A is doubled and ρ is periodic and resumes its original value). a varies from 1 to 2 in one cell cycle.
Similarly, instead of R, L, and ω, it proves convenient to use the following set of transformed parameters: y=sinω,z=y+LR+1 (i.e., the cell half length divided by R), and R. Simple geometrical calculation shows that for a given reduced volume v, the parameters y and z are related by
y=3z12vz3213=fvz
(8)
and that for a given normalized area a, the values of R and z are related by
R=Rt0az.
(9)

For given v and a, Equations 89 thus define a continuum of possible shapes varying between two extremes:

  • 1. 

    Elongated rod: yel = 0 and f(v, zel) = 0.

  • 2. 

    Most constricted shape:

    • (a) 

      ycs = 1 and f(v, zcs) = 1 if v12 (two elongated rods attached by an infinitesimal neck).

    • (b) 

      ycs = zcs − 1 and f(v, zcs) = zcs − 1 if v>12 (two overlapping spheres).

The limit values at the beginning (t = t0) and at the end of the cell cycle (t = t0 + Tdbl) are given by:
  • 1. 

    a(t0) = 1, v(t0) = 1, y(t0) = 0, z(t0) = 1;

  • 2. 

    at0+Tdbl=2,vt0+Tdbl=12,yt0+Tdbl=1,zt0+Tdbl=2.

The continuum of possible shapes may be spanned by letting y vary between the two extremes, yel = 0 ≤ yycs, and by determining the corresponding z and R from Equations 8 and 9, respectively. For each of these possible shapes, simple geometrical calculation gives the bending energy as
Eb=2πκ1+y2C0R2+z1y1C0R2.
(10)
The actual cellular shape is the one that minimizes Eb.7

3.3 Determining the Reduced Volume and Normalized Area Trajectories

Defining η as the logarithmic derivative of v, logarithmic differentiation of Equation 6 gives
η=v˙v=ρ˙ρ12A˙A
(11)
Here ρ˙/ρ can be computed from the ρ trajectory, and A˙/A is given by Equation 3 and can be computed from the membrane precursor concentration trajectory. If an oscillatory regime is reached for asymptotic trajectories of c(t) and ρ(t), then η(t) is also periodic with the same period Tperiod.8
The trajectory v(t) is given by integration of η:
vt=expt0tηtdt
(12)
and the normalized area trajectory is given by integration of Equation 3:
at=expt0tFouttNmedt.
(13)

3.4 The Spontaneous Curvature C0 Should Be Sufficiently Large

C0 reflects the tendency of the cell to spontaneously take a spherical shape with radius 2/C0, for which the bending energy reaches its absolute minimum (Eb = 0). At the end of a cell cycle, the double sphere should have a lower bending energy than the elongated rod, else the cell will not divide even though both volume and surface area have been doubled. Simple geometrical calculation shows that, within the framework of our simplified family of shapes, the following inequality must be met:
C0Rt0>xthr=1233422+33+31.18
(14)
C0 must thus be higher than a minimum value, C0>C0,min=xthr/Rt0. Independently of the exact threshold value,9 the existence of such a minimum value for C0 is consistent with the rigorous phase diagram determinations of Seifert et al. [34]. It is also consistent with the previous work of Surovstev et al. [38], where it was found that a contractile ring was necessary to result in cellular division, because a null spontaneous curvature was implicitly assumed in their work.

C0 reflects a structural property of the membrane, and Rt0=3/ρt0 reflects the embedded chemistry and transmembrane kinetics (because the trajectory of ρ is the result of the differential equation system (1), (2), and (4)). The above constraint (the inequality (14)) could be satisfied as the result of a coevolution of the membrane and of the embedded metabolism.

4 Asymptotic Trajectories and Resulting Cellular Shapes

4.1 Stationary Case

In such a case, vt=expλt/2 and a(t) = exp(λt), where λ=log2/Tdbl is the constant growth rate. Applying the procedure described in Section 3.2 results in the cell shape trajectory shown in Figure 3, for various C0 values above the threshold value.

Figure 3. 

Stationary case where all concentrations and the surface-area-to-volume ratio are asymptotically constant. Cell shapes are shown at various time points across the cell cycle, for four C0R0 values.

Figure 3. 

Stationary case where all concentrations and the surface-area-to-volume ratio are asymptotically constant. Cell shapes are shown at various time points across the cell cycle, for four C0R0 values.

The shape pattern only depends on the product C0R0, where R0=3/ρ0 with ρ0 being the stationary surface-area-to-volume ratio. Only the size and time scales depend on the specific stationary asymptotic values: R0 for the size scale, and Tdbl=Nmelog2/Fout (obtained by integration of Equation 3 with constant Fout because the concentration vector is constant) for the time scale.

In the range xthr < C0R0 ≤ 1.05 × xthr, the shape trajectory is nearly identical to that of C0R0 = 1.05 × xthr (Figure 3a), with constriction only occurring towards the end of the cell cycle. In the range C0R0 ≥ 1.35 × xthr, it is nearly identical to that of C0R0 = 1.35 × xthr (Figure 3d), with constriction occurring early in the cell cycle. The average bending energy (averaging over one cell cycle) was found to reach a minimum for C0R0 ≈ 1.82 × xthr.

The cells do self-replicate by growth and division, although the surface-area-to-volume ratio remains constant. This is because spherical shapes are preferred for C0R0 > xthr. On the contrary, if C0R0 < xthr, then starting from a spherical shape with radius R0=3/ρ0, the cell keeps elongating without constriction towards a filament shape with radius R0=2/ρ0 (infinitely long cylinder).

4.2 Oscillating Case

Various simple oscillating CRNs have been proposed and investigated. The one proposed by Wilhelm and Heinrich [42] was chosen as an example in the present work. The rationale for this choice is detailed in Appendix A.1. This CRN is the following:
A+X2XX+YY+BYBXZZY

This system is conservative (all masses are equal), and the directionality of each reaction corresponds to the following ranking of Gibbs energies: gA > gX > gZ > gY > gB. Here A was chosen as the nutrient, and either Y or B (species having the lowest Gibbs energies) as the membrane precursor. Gibbs energies should be high for the nutrient and low for the membrane precursor for the growth rate to be high [8].

The first reaction A + X → 2X was taken as the nutrient import reaction. The associated reaction vector Fin has as single nonzero component Fin,X = KcAcX, where K is a surfacic kinetic constant and cA and cX are the concentrations of A and X respectively. A is only present in the outside growth medium (assumed sufficiently large so that cA remains approximately constant), and is transformed across the membrane into X under the activation of X itself. Such a scheme is reminiscent of the phosphoenoltransferase sugar (PTS) active transport system in E. coli, for which glucose is transformed across the membrane into glucophosphate under the activation of phosphoenolpyruvate (see Appendix A.1).

The incorporation of the membrane precursor is assumed unidirectional and kinetically controlled with coefficient Kme per unit area, so that Fout has as single nonzero component Fout,Y = KmecY if Y is the membrane precursor (or Fout,B = KmecB if B is the membrane precursor).10

All chemical species are assumed to contribute to one particle when dissolved in solution, ui = 1, i = A, X, Y, Z, or B.11

Applying Equations 1, 2, and 4 to this example CRN results in the following ODE system:

  • 1. 
    If Y is the membrane precursor,
    c˙X=ρKcAcXk1cXcYk3cXλcXc˙Y=k2cY+k4cZρKmecYλcYc˙Z=k3cXk4cZλcZc˙B=k1cXcY+k2cYλcBρ˙ρ=KmecYNmeλ
    (15)
    where the growth rate λ is given by
    λ=ρKcAcXKmecYcA.
    (16)
  • 2. 
    If B is the membrane precursor,
    c˙X=ρKcAcXk1cXcYk3cXλcXc˙Y=k2cY+k4cZλcYc˙Z=k3cXk4cZλcZc˙B=k1cXcY+k2cYρKmecBλcBρ˙ρ=KmecBNmeλ
    (17)
    where the growth rate λ is given by
    λ=ρKcAcXKmecBcA.
    (18)

In the above equations, cA (cX, cY, cZ, cB) is the concentration of A (of X, Y, Z, B), and k1 (k2, k3, k4) is the kinetic coefficient for the reaction X + YY + B (the reaction YB, XZ, ZY).

Appendix A.2 describes the numerical method used to identify four examples of parameter sets resulting in oscillatory behavior and meeting the synchronization constraint (that the ratio of doubling time to oscillation period must be an integer). These parameter sets are given in Table 2.

Figure 4 gives the asymptotic trajectories for concentrations and the surface-area-to-volume ratio for these examples. The doubling time ranges from 900 to 9000 s, and the typical cell radius (approximately given by 3/ρavg, assuming a spherical shape, where ρavg is the average surface-area-to-volume ratio, averaging over one cell cycle) from 1.5 to 15 μm. These values are biologically plausible.

Figure 4. 

Asymptotic trajectories for the concentrations (top) and the surface-area-to-volume ratio ρ (bottom), for examples Y1 (a), Y2 (b), B1 (c), and B2 (d). These asymptotic trajectories are shown over an arbitrarily chosen number four of doubling times to illustrate the oscillatory behavior. For examples Y1 and B1, the doubling time is equal to the period, whereas for examples Y2 and B2, it is equal to twice the period.

Figure 4. 

Asymptotic trajectories for the concentrations (top) and the surface-area-to-volume ratio ρ (bottom), for examples Y1 (a), Y2 (b), B1 (c), and B2 (d). These asymptotic trajectories are shown over an arbitrarily chosen number four of doubling times to illustrate the oscillatory behavior. For examples Y1 and B1, the doubling time is equal to the period, whereas for examples Y2 and B2, it is equal to twice the period.

The trajectories v(t) and a(t) are first integrated as described in Section 3.3 and are shown in Figure 5; and the cellular shape trajectories are then derived as described in Section 3.2, and are shown in Figure 6.

Figure 5. 

Reduced volume and normalized area trajectories over one doubling time, integrated from Equations 1213. t0 and the spontaneous curvature C0 were chosen so as to minimize the average bending energy (averaging over one doubling time) while meeting all constraints (see text). The corresponding t0 and C0 choices are given in the headings of Figure 6.

Figure 5. 

Reduced volume and normalized area trajectories over one doubling time, integrated from Equations 1213. t0 and the spontaneous curvature C0 were chosen so as to minimize the average bending energy (averaging over one doubling time) while meeting all constraints (see text). The corresponding t0 and C0 choices are given in the headings of Figure 6.

Figure 6. 

Cell shapes at various time points across the cell cycle. t0 and C0 were chosen so as to minimize the average bending energy (averaging over one cell cycle) while meeting all constraints (see text). Note that the scale is preserved within each subfigure corresponding to the cell cycle of each example, but not across different subfigures: Figure 4 shows that ρ(t0) ranges from ≈0.2 μm−1 (example Y1) to ≈2.5 μm−1 (example B1), corresponding to a sphere radius (at the beginning of the cell cycle) ranging from ≈1.2 μm (example B1) to ≈15 μm (example Y1).

Figure 6. 

Cell shapes at various time points across the cell cycle. t0 and C0 were chosen so as to minimize the average bending energy (averaging over one cell cycle) while meeting all constraints (see text). Note that the scale is preserved within each subfigure corresponding to the cell cycle of each example, but not across different subfigures: Figure 4 shows that ρ(t0) ranges from ≈0.2 μm−1 (example Y1) to ≈2.5 μm−1 (example B1), corresponding to a sphere radius (at the beginning of the cell cycle) ranging from ≈1.2 μm (example B1) to ≈15 μm (example Y1).

Contrary to the stationary case, the trajectories v(t) and a(t) depend on the choice t0 of the beginning of a cell cycle (relative to an oscillation period). v(t) may decrease monotonically (stationary behavior or mild oscillations), or not, with η(t) spanning both negative and positive values (case of stronger oscillations as seen in numerical examples). Such a nonmonotonic variation of v(t) imposes restrictions on the precise starting point t0 of the cell cycle to ensure v(t) ≤ 1 at any time during the cell cycle, and may additionally impose an upper limit on the spontaneous curvature C0:

  • 1. 

    Appendix A.3 explains why t0 must satisfy t0,min < t0 < t0,max in order to ensure that v(t) ≤ 1 at all times across the cell cycle, and gives the integral relations from which t0,min and t0,max can be determined.

  • 2. 

    Besides, if at some critical point before cycle end, say tc < t0 + Tdbl, the reduced volume reaches vtc=1/2, then the double sphere should have a higher bending energy than the elongated rod, so that the cell may not divide too early into two smaller cells (smaller than the original cell). The same reasoning as above leads to the following inequality: C0Rc < xthr, where Rc is the sphere radius for the double sphere having a(tc) < 2 as normalized area. Equation 9 gives Rc=Rt0atc/2, because z = 2 for a double sphere. The above inequality may thus be rewritten as C0Rt0<xthr2/atc. Whether v(t) may actually reach 1/2 before cycle end or not depends on the amplitude of η(t) oscillations as detailed in Appendix A.4. The existence of such a critical point tc constrains C0 and R(t0) to be closely related, as C0R(t0) should satisfy xthr<C0Rt0<xthr2/atc, and 2/atc may be only slightly above unity. This is the case of examples Y1, Y2, and B2, but not of example B1, which displays milder oscillations.

Figure 6 shows the resulting cellular shapes, for t0 and C0 chosen so as to minimize the average bending energy (averaging over one cell cycle) while meeting the above constraints. For all four examples, this minimum was found for t0 close to its upper bound t0,max. This t0 is also the origin from which the trajectories shown in Figure 4 were plotted. For examples Y1, Y2, and B2, this minimum was found for C0 close to its upper bound, while for example B1 (for which C0 is only constrained by a lower bound), it was found for C0 close to twice its lower bound. It can be seen that the cells do indeed self-replicate by growth and division, and that the various oscillation patterns observed in Figure 4 result in a variety of cellular shape trajectories across the cell cycle.

5 Discussion

5.1 Consistency of Oscillatory Behavior with a Previous Theoretical Proof on Possible Cell Trajectories

For the example CRN of Section 4.2, the focus was on selected parameter sets resulting in oscillatory asymptotic trajectories (for concentrations and for the surface-area-to-volume ratio). Other parameter sets were found to result in constant strictly positive values (stationary behavior), in the divergence of ρ towards infinity (cell collapsing onto itself), or in the extinction of the membrane precursor. These outcomes are consistent with the proof in [7]: ρ is either unbounded or bounded. In the latter case, all variables that are solutions of the ODE system defined by Equations 1, 2, and 4 are bounded (concentrations are bounded by the constant osmolarity), which assures the existence of a fixed point. The membrane precursor concentration at this fixed point is either null (no growth) or strictly positive. In the latter case, the fixed point is either attracting (stationary behavior), oscillating (this work), or bounded chaotic (yet to be observed numerically).

With randomly generated chemical reaction networks, no oscillating behavior was observed [7]. This is not surprising, because both network topology and reaction kinetics were random, whereas oscillatory behavior requires both specific topologies and specific kinetics.

5.2 Limitations Arising from the Simplified Predefined Family of Shapes

Quantitative findings of this work are most likely influenced by the arbitrary choice of a predefined family of shapes. For example, the minimum value of xthr ≈ 1.18 for C0R(t0) depends on this choice, as illustrated by the alternative value for a different shape comparison in footnote 9. But the existence of such a minimum is expected to be quite general, and the chosen predefined family of shapes helps understand why. C0 reflects the tendency of a cell to adopt a spherical shape with radius 2/C0. If the starting cell is too small (i.e., its curvature is too high), it may not divide after having doubled in size, but continue growing instead.

The same holds for the maximum value of C0R(t0) in the case of oscillations with large amplitude such that the reduced volume may reach 1/2 before cycle end. The precise value of this threshold most likely depends on the choice of a predefined family of shapes, but the existence of such a maximum is expected to be quite general. If the starting cell is too large, it may divide early during the cell cycle, provided such a division be possible.

Whether the amplitude of biochemical oscillations in modern evolved cells is low or high (and whether an upper bound for C0R(t0) would be applicable, as is the case for examples Y1, Y2, and B2, but not for B1) is still an open question.12 Fairly large oscillations have been reported in [25], with concentrations in human cells varying across the cell cycle by a factor > 2 for H+, NAD, and NADH, and > 4 for NADP and NADPH. But this data was only available for a very limited sample of the total metabolite pool. Besides, there is always the possibility that the variation of concentrations across the cell cycle is due not only to nonlinear dynamics in the metabolic reaction network (assuming constant metabolic kinetics, which corresponds to a fixed gene expression pattern), but also to more complex physico-chemical interactions such as those involved in gene regulation mechanisms, the net effect of which would amount to varying certain metabolic kinetics across the cell cycle (by varying the level of enzyme synthesis).

5.3 Applicability of Findings from This Work to Other Oscillating Chemical Systems or Whole-Cell Models

All four numerical examples of Section 4.2 rely upon a single example of an oscillating CRN. However, this example is only a key to illustrate principles and constraints that have wider applicability.

The constraint of having the doubling time equal to a multiple of the oscillation period is not only applicable to any oscillating CRN embedded within the present model, but it also holds for any self-replicating system and any whole-cell model where both surface area and volume are doubled in one cycle time while concentrations resume their original values after one cycle time (see in particular footnote 3).

Regardless of the specific numbers (see the discussion in Section 5.2 above), the constraint of having the spontaneous curvature and cell size adapted (at least through the requirement of having C0R(t0) above a minimum value, and possibly also below a maximum value) is expected to hold not only for any other oscillating chemical reaction system, but also and more generally for other whole-cell models. It is always conceivable that a cell wall or a cytoskeleton may contribute to alleviating such constraints. However, such structures do not play the same role in any type of cell (e.g., eukaryotes versus prokaryotes) or inner compartment (e.g., mitochondria).

Such a co-adaptation of C0 and R(t0) suggests a coevolution of the membrane and of the embedded metabolism. This is because C0 depends on membrane structural properties, while R(t0) depends on the biochemical and transmembrane kinetics.

5.4 Oscillations Appear More Detrimental than Beneficial

Monitoring trajectories over several generations (with cell division whenever y = 1; and upon division: normalized area divided by 2, reduced volume multiplied by 2, concentrations and surface-area-to-volume trajectories unchanged) reveals the following patterns:

  • 1. 

    When the constraints on the spontaneous curvature are not met, growth and division may still occur over the course of a few generations, and then either division stops (because C0R(t0 + pTdbl)—where p is the number of generations—is too low), in which case, the cell still continues to elongate without dividing, or bursts (i.e., suffers reduced volume v > 1, because C0R(t0 + pTdbl) is too high).

  • 2. 

    When the constraint on the doubling time (being a multiple of the period) is not met, growth and division may still occur for a few generations, but the resulting cells are more irregular in size, shape, and composition, with the same possible end fates as above.

  • 3. 

    When the constraints on the doubling time and on the spontaneous curvature are met, harmonious growth and division occurs over a large number of generations (50 for example Y2, over 500 for examples Y1, B1, and B2), but in all tested numerical examples, even slight deviations from perfect periodicity make the system ultimately deviate from the ideal self-replicator trajectory, and cells end up either stopping growing or bursting, as above.

Cells based on oscillating chemical reaction networks are thus inherently fragile compared to their stationary counterparts, because (i) the above synchronization cannot be perfect, and (ii) even minor environmental fluctuations (e.g., change in temperature, nutrient concentration, or outside osmolarity) could result in system parameters drifting away from the synchronization condition. This lack of stability is obviously a disadvantage.

Although the requirement of having the doubling time a multiple of the oscillation period appears to be a strong constraint, it should be noted that this outcome (of monitoring trajectories over several generations) could not be theoretically predicted. In particular, the possibility of self-stabilization could not be excluded a priori (potentially with irregular cell size and shape; see point 2 above).

Additional feedback mechanisms (whereby a change in kinetic constant would itself induce another change compensating the first one) would be required for such a system to be robust. Because of such required additional complexity, and because it has previously been observed that cells based upon randomly generated CRNs typically reach a non-oscillating stationary regime [7], self-replication of early protocells may not have relied upon chemical oscillations.

Biochemical oscillations may thus have emerged later in the evolution, for other reasons. One such reason is a potential metabolic advantage: Autocatalytic mechanisms such as those of the citric acid cycle are advantageous because they significantly increase the number of ATP molecules produced per consumed molecule of glucose, compared to depending only on glycolysis. A side consequence of such autocatalytic mechanisms is the possibility of biochemical oscillations. Another possible reason has to do with the peptidoglycan cell wall. Such a wall presumably gave an evolutionary advantage because of the increased mechanical robustness. But a side consequence of such a rigid wall might have been the risk of having the cell remain in a metastable elongated rod shape rather than dividing (even though the bending energy of the double cell is lower than that of the elongated rod) because the bending-energy barrier height to be overcome increases with the bending rigidity (coefficient κ in Equation 5). It is thus possible that the appearance of such a cell wall might have coincided with the emergence of more complex division machinery, requiring biochemical oscillations for precise clocking and cell size control.

5.5 Limitations Arising from Protocell Model Assumptions

The fact that our conclusions are influenced by an oversimplistic cell model cannot be ignored:

  • 1. 

    Our model assumes that chemical reactions taking place inside the cell can be successfully modeled by a set of deterministic ODEs. There is extensive experimental evidence of stochastic effects in living cells, especially for chemical species in very low concentrations [15], and more complex stochastic models have also been proposed to account for them [43].

  • 2. 

    Our model also assumes that all chemical concentrations are spatially homogeneous inside and outside the cell, and that the nutrient concentration remains constant at all time. This is not the case in real biological cells, where spatial organization is known to play an important role [18].

Whether stochastic effects or spatial organization play a determinant role in explaining fundamental mechanisms of cell growth and division and influence the conclusions of this work, remains an open question.

6 Conclusion

This simple cell model shows that biochemical oscillations are not required for cellular self-replication by growth and division, although a greater variety of cellular shapes can be obtained with biochemical oscillations. In both cases, the spontaneous curvature should be higher than a minimum value for the double sphere to be more stable than an elongated rod having the same surface area and the same volume. Besides, additional constraints must be met in the case of chemical oscillations: (i) the doubling time should be equal to, or a multiple of, the oscillation period; and (ii), depending on the amplitude of oscillations, the spontaneous curvature should be lower than a maximum value. The former synchronization constraint suggests additional feedback mechanisms must be at play to make cells robust self-replicators. The latter spontaneous-curvature constraint suggests a coevolution of the membrane and of the metabolism, to adapt the membrane spontaneous curvature to the surface-area-to-volume ratio, which is determined by chemical kinetics. Because of these additional constraints, the most stringent of which is the synchronization condition, early self-replicating protocells may not have relied upon chemical oscillations. Biochemical oscillations may have appeared later because of different evolutionary advantages. Beyond protocells, this work also stresses the importance of taking into account cellular shape and membrane growth in any whole-cell model.

Acknowledgments

The authors would like to thank Alain Hamel, Laurent Schwartz, and Jean-Yves Trosset for stimulating discussions, as well as anonymous reviewers for constructive comments.

Notes

1 

Here, the scalar Fout is the component of the vector Fout along the membrane precursor.

2 

ζAt+Tdbl=expλavgtexpλavgTdblAt+Tdbl=expλavgt×12×At+Tdbl=expλavgtAt=ζAt, and similarly ζVt+Tdbl=ζVt.

3 

In such a most general case, the shortest period need not be the same for ζA(t), ζA(t), ρ(t), and concentrations. This is the reason why Tperiod was defined here as the shortest period common to all. In the present cell model, the concentrations, ζA(t), ζV(t), and ρ(t) are all coupled by the ODE system given by Equations 13, and all these quantities oscillate with the same shortest period Tperiod.

4 

See in particular Figure 18 in [33]. On a (v, c0) diagram, where v is the reduced volume and c0 is the scaled spontaneous curvature defined as the product of the spontaneous curvature C0 and the radius of the sphere enclosed by A, trajectories for examples of the present work span areas of the phase diagram where the shape of lowest energy has yet to be determined.

5 

One of the two principal curvatures is non-continuous at the junction between the two half cells, which may result in an infinite contribution to the total membrane energy. The present approach thus only provides a crude approximation of the exact shape trajectory, for which the constriction pattern should be rounded and without such discontinuity.

6 

The bending rigidity of the membrane is implicitly assumed to be sufficiently low so that thermal fluctuations are sufficient for the cell to relax to its shape of lowest membrane energy regardless of the height of any energy barrier to be overcome. Else, the cell might stay in a metastable state.

7 

In all investigated numerical examples, this minimum always lay at one of the two extremes (elongated rod or most constricted shape), but not in between. The bending energy either varied monotonically between these two extremes or exhibited a maximum in between. When, along the (v, a) trajectory, the shape having minimum bending energy varies from one extreme to the other, no energy barrier is to be overcome if the bending energy varies monotonically between these two extremes. Conversely, if the bending energy exhibits its maximum in between, the resulting energy barrier may be crossed as the result of thermal fluctuations; see footnote 6.

8 

The logarithmic derivative of A = ζA × exp(λavgt) is A˙/A=ζ˙A/ζA+λavg. Besides, ρ˙/ρ=ζ˙A/ζAζ˙V/ζV. This gives η=ζ˙V/ζV3/2ζ˙A/ζA1/2λavg, which is periodic because the derivative of a periodic function is periodic.

9 

The precise value of this threshold depends on the choice of predefined shapes. To illustrate this point, it can be easily shown that an infinite string of spherical beads (with radius R(t0)) connected by infinitely narrow necks is more stable than the infinite cylinder having the same area and the same volume if C0Rt0>7/4.

10 

Such a single unidirectional kinetic step is a rough simplification, as much more complex mechanisms have been evidenced in artificial protocell membranes [13]. Likewise, the flip-flop of lipids between the two layers of the membrane is assumed sufficiently fast [35].

11 

For this example, the osmolarity vector u is parallel to the mass vector m, which results in uTS = 0 because the CRN is conservative. This simplifies the expression of Equation 2 for the growth rate.

12 

It is recalled that this refers to the amplitude of oscillations in η(t) obtained by integration of Equation 11. In the case of one oscillation period per doubling time, η(t) typically exhibits two lobes of opposite signs, the area under the negative lobe being larger (in absolute value) than that under the positive lobe (because the integral over one doubling time is negative, equal to 1/2). By a high amplitude we mean that the area under the negative lobe is larger (in absolute value) than 2×1/2=2 (see Appendix A.4).

13 

The existence of a limit cycle for the simple chemostat model as in references given in Table 1 does not necessarily translate into the existence of a limit cycle for the cell model described in Section 2, because the ODE system is different. A detailed investigation of the dependence of cell trajectories on the initial conditions has been left for further study.

14

Replacing the reaction B + XY + D by the reaction XY reduces the number of moieties to unity, but this modified system loses its nondegenerate fixed point.

15 

Setting all derivatives in the equation system (19) to zero, the last equation gives l = ϕdY. Feeding this into the equation before last gives dB=dX+ψ1/ϕ>ψ1/ϕ. Feeding this into the constant osmolarity gives dA=dX+dY+dZ+dB>dB>ψ1/ϕ, or equivalently, the inequality (24).

16 

This ensures that the inequality (24) is verified for any tested parameter set.

17 

The Fourier transformation was carried out over a minimum of 100 periods after reaching the asymptotic regime.

References

1
Adamala
,
K.
, &
Szostak
,
J. W.
(
2013
).
Competition between model protocells driven by an encapsulated catalyst
.
Nature Chemistry
,
5
(
6
),
495
501
.
2
Angeli
,
D.
,
Leenheer
,
P. D.
, &
Sontag
,
E.
(
2007
).
A Petri net approach to the study of persistence in chemical reaction networks
.
Mathematical Biosciences
,
210
(
2
),
598
618
.
3
Bennett
,
B.
,
Kimball
,
E.
,
Gao
,
M.
,
Osterhout
,
R.
,
Dien
,
S. V.
, &
Rabinowitz
,
J.
(
2009
).
Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli
.
Nature Chemical Biology
,
5
(
8
),
593
599
.
4
Bennett
,
B.
,
Yuan
,
J.
,
Kimball
,
E.
, &
Rabinowitz
,
J.
(
2008
).
Absolute quantitation of intracellular metabolite concentrations by an isotope ratio-based approach
.
Nature Protocols
,
3
(
8
),
1299
1311
.
5
Bigan
,
E.
,
Paulevé
,
L.
,
Steyaert
,
J.-M.
, &
Douady
,
S.
(
2015
).
Necessary and sufficient conditions for protocell growth
.
Journal of Mathematical Biology
,
73
(
6–7
),
1
38
.
6
Bigan
,
E.
,
Steyaert
,
J.-M.
, &
Douady
,
S.
(
2015
).
Chemical schemes for maintaining different compositions across a semi-permeable membrane with application to proto-cells
.
Origins of Life and Evolution of Biospheres
,
45
(
4
),
439
454
.
7
Bigan
,
E.
,
Steyaert
,
J.-M.
, &
Douady
,
S.
(
2015
).
Filamentation as primitive growth mode?
Physical Biology
,
12
(
6
),
066024
.
8
Bigan
,
E.
,
Steyaert
,
J.-M.
, &
Douady
,
S.
(
2015
).
Minimal conditions for proto-cell stationary growth
.
Artificial Life
,
21
(
2
),
166
192
.
9
Božič
,
B.
, &
Svetina
,
S.
(
2004
).
A relationship between membrane properties forms the basis of a selectivity mechanism for vesicle self-reproduction
.
European Biophysics Journal
,
33
(
7
),
565
571
.
10
Božič
,
B.
, &
Svetina
,
S.
(
2007
).
Vesicle self-reproduction: The involvement of membrane hydraulic and solute permeabilities
.
The European Physical Journal E
,
24
(
1
),
79
90
.
11
Budin
,
I.
,
Debnath
,
A.
, &
Szostak
,
J. W.
(
2012
).
Concentration-driven growth of model protocell membranes
.
Journal of the American Chemical Society
,
134
(
51
),
20812
20819
.
12
Canham
,
P. B.
(
1970
).
The minimum energy of bending as a possible explanation of the biconcave shape of the human red blood cell
.
Journal of Theoretical Biology
,
26
(
1
),
61
81
.
13
Chen
,
I. A.
, &
Szostak
,
J. W.
(
2004
).
A kinetic study of the growth of fatty acid vesicles
.
Biophysical Journal
,
87
(
2
),
988
998
.
14
Deuling
,
H.
, &
Helfrich
,
W.
(
1976
).
The curvature elasticity of fluid membranes: A catalogue of vesicle shapes
.
Journal de Physique
,
37
(
11
),
1335
1345
.
15
Elowitz
,
M. B.
,
Levine
,
A. J.
,
Siggia
,
E. D.
, &
Swain
,
P. S.
(
2002
).
Stochastic gene expression in a single cell
.
Science
,
297
(
5584
),
1183
1186
.
16
Fanelli
,
D.
, &
McKane
,
A. J.
(
2008
).
Thermodynamics of vesicle growth and instability
.
Physical Review E
,
78
(
5
),
051406
.
17
Field
,
R. J.
, &
Noyes
,
R. M.
(
1974
).
Oscillations in chemical systems. IV. Limit cycle behavior in a model of a real chemical reaction
.
The Journal of Chemical Physics
,
60
(
5
),
1877
1884
.
18
Goodsell
,
D. S.
(
2009
).
Escherichia coli
.
Biochemistry and Molecular Biology Education
,
37
(
6
),
325
332
.
19
Gray
,
P.
, &
Scott
,
S.
(
1984
).
Autocatalytic reactions in the isothermal, continuous stirred tank reactor: Oscillations and instabilities in the system A + 2B → 3B; BC
.
Chemical Engineering Science
,
39
(
6
),
1087
1097
.
20
Huijbregts
,
R.
,
de Kroon
,
A.
, &
de Kruijff
,
B.
(
2000
).
Topology and transport of membrane lipids in bacteria
.
Biochimica et Biophysica Acta (BBA)—Reviews on Biomembranes
,
1469
(
1
),
43
61
.
21
Lotka
,
A. J.
(
1920
).
Undamped oscillations derived from the law of mass action
.
Journal of the American Chemical Society
,
42
(
8
),
1595
1599
.
22
Macía
,
J.
, &
Solé
,
R. V.
(
2007
).
Protocell self-reproduction in a spatially extended metabolism–vesicle system
.
Journal of Theoretical Biology
,
245
(
3
),
400
410
.
23
Martín
,
O.
,
Peñate
,
L.
,
Alvaré
,
A.
,
Cárdenas
,
R.
, &
Horvath
,
J.
(
2009
).
Some possible dynamical constraints for life's origin
.
Origins of Life and Evolution of Biospheres
,
39
(
6
),
533
544
.
24
Mavelli
,
F.
, &
Ruiz-Mirazo
,
K.
(
2013
).
Theoretical conditions for the stationary reproduction of model protocells
.
Integrative Biology
,
5
(
2
),
324
341
.
25
Moreira
,
J. d. V.
,
Hamraz
,
M.
,
Abolhassani
,
M.
,
Bigan
,
E.
,
Pérès
,
S.
,
Paulevé
,
L.
,
Nogueira
,
M. L.
,
Steyaert
,
J.-M.
, &
Schwartz
,
L.
(
2016
).
The redox status of cancer cells supports mechanisms behind the Warburg effect
.
Metabolites
,
6
(
4
),
33
.
26
Neidhardt
,
F.
(
1996
).
Escherichia coli and Salmonella: Cellular and molecular biology
.
Washington
:
ASM Press
.
27
Novák
,
B.
, &
Tyson
,
J. J.
(
2008
).
Design principles of biochemical oscillators
.
Nature Reviews Molecular Cell Biology
,
9
(
12
),
981
991
.
28
Prigogine
,
I.
, &
Lefever
,
R.
(
1968
).
Symmetry breaking instabilities in dissipative systems. II
.
The Journal of Chemical Physics
,
48
(
4
),
1695
1700
.
29
Sacerdote
,
M.
, &
Szostak
,
J.
(
2005
).
Semipermeable lipid bilayers exhibit diastereoselectivity favoring ribose
.
Proceedings of the National Academy of Sciences of the U.S.A.
,
102
(
17
),
6004
6008
.
30
Schlögl
,
F.
(
1972
).
Chemical reaction models for non-equilibrium phase transitions
.
Zeitschrift für Physik
,
253
(
2
),
147
161
.
31
Schnakenberg
,
J.
(
1979
).
Simple chemical reaction systems with limit cycle behaviour
.
Journal of Theoretical Biology
,
81
(
3
),
389
400
.
32
Schuster
,
S.
, &
Höfer
,
T.
(
1991
).
Determining all extreme semi-positive conservation relations in chemical reaction systems: A test criterion for conservativity
.
Journal of the Chemical Society, Faraday Transactions
,
87
(
16
),
2561
2566
.
33
Seifert
,
U.
(
1997
).
Configurations of fluid membranes and vesicles
.
Advances in Physics
,
46
(
1
),
13
137
.
34
Seifert
,
U.
,
Berndl
,
K.
, &
Lipowsky
,
R.
(
1991
).
Shape transformations of vesicles: Phase diagram for spontaneous-curvature and bilayer-coupling models
.
Physical Review A
,
44
(
2
),
1182
.
35
Shirt-Ediss
,
B.
,
Ruiz-Mirazo
,
K.
,
Mavelli
,
F.
, &
Solé
,
R. V.
(
2014
).
Modelling lipid competition dynamics in heterogeneous protocell populations
.
Scientific Reports
,
4
(
5675
).
36
Shirt-Ediss
,
B.
,
Solé
,
R. V.
, &
Ruiz-Mirazo
,
K.
(
2015
).
Emergent chemical behavior in variable-volume protocells
.
Life
,
5
(
1
),
181
211
.
37
Surovstev
,
I. V.
,
Morgan
,
J. J.
, &
Lindahl
,
P. A.
(
2007
).
Whole-cell modeling framework in which biochemical dynamics impact aspects of cellular geometry
.
Journal of Theoretical Biology
,
244
(
1
),
154
166
.
38
Surovtsev
,
I. V.
,
Zhang
,
Z.
,
Lindahl
,
P. A.
, &
Morgan
,
J. J.
(
2009
).
Mathematical modeling of a minimal protocell with coordinated growth and division
.
Journal of Theoretical Biology
,
260
(
3
),
422
429
.
39
Svetina
,
S.
, &
Žekš
,
B.
(
1989
).
Membrane bending energy and shape determination of phospholipid vesicles and red blood cells
.
European Biophysics Journal
,
17
(
2
),
101
111
.
40
Svetina
,
S.
, &
Žekš
,
B.
(
2002
).
Shape behavior of lipid vesicles as the basis of some cellular processes
.
The Anatomical Record
,
268
(
3
),
215
225
.
41
Walde
,
P.
,
Wick
,
R.
,
Fresta
,
M.
,
Mangone
,
A.
, &
Luisi
,
P. L.
(
1994
).
Autopoietic self-reproduction of fatty acid vesicles
.
Journal of the American Chemical Society
,
116
(
26
),
11649
11654
.
42
Wilhelm
,
T.
, &
Heinrich
,
R.
(
1995
).
Smallest chemical reaction system with Hopf bifurcation
.
Journal of Mathematical Chemistry
,
17
(
1
),
1
14
.
43
Wilkinson
,
D. J.
(
2009
).
Stochastic modelling for quantitative description of heterogeneous biological systems
.
Nature Reviews Genetics
,
10
(
2
),
122
133
.
44
Zhu
,
T. F.
, &
Szostak
,
J. W.
(
2009
).
Coupled growth and division of model protocell membranes
.
Journal of the American Chemical Society
,
131
(
15
),
5705
5713
.

Appendix

A.1 Choice of Oscillating CRN

A number of simple oscillating chemical reaction networks have been proposed. A non-exhaustive selection is presented in Table 1, along with selection criteria, as explained below. For all of them, high-energy compounds are burnt into low-energy compounds and intermediate compounds are generated in this process. These systems are oscillating in the sense that if the concentrations of the high-energy compounds are kept constant (chemostat model), the concentrations of intermediate compounds oscillate, while the concentrations of the low-energy compounds increase indefinitely. Such oscillations are only observed within a certain range of kinetic parameters and nutrient concentration, which have been characterized in references given in Table 1.

Table 1. 

Candidate oscillating chemical reaction networks. The one proposed by Wilhelm and Heinrich [42] was chosen in the present work.

System name Chemical reactions Limit cycle Single moiety Short siphon Reference 
Lotka-Volterra A + X → 2X No Yes Yes [21
X + Y → 2Y 
YB 
 
Schlögl A + 2X → 3X No (only bistability) Yes Yes [30
3X → 2X + B 
XB 
 
Brusselator AX Yes No (2 moieties) No (a siphon containing X also contains A[28
2X + Y → 3X 
B + XY + D 
XE 
 
Wilhelm A + X → 2X Yes Yes Yes [42
X + YY + B 
YB 
XZ 
ZY 
System name Chemical reactions Limit cycle Single moiety Short siphon Reference 
Lotka-Volterra A + X → 2X No Yes Yes [21
X + Y → 2Y 
YB 
 
Schlögl A + 2X → 3X No (only bistability) Yes Yes [30
3X → 2X + B 
XB 
 
Brusselator AX Yes No (2 moieties) No (a siphon containing X also contains A[28
2X + Y → 3X 
B + XY + D 
XE 
 
Wilhelm A + X → 2X Yes Yes Yes [42
X + YY + B 
YB 
XZ 
ZY 

The selection criteria are the following:

  • 1. 

    Limit cycle: All systems exhibit a limit cycle except the Lotka-Volterra and the Schlögl systems. A system exhibiting a limit cycle has been preferred in this work because the resulting oscillations are independent of the initial conditions, which is not the case for the Lotka-Volterra system. The Schlögl system exhibits bistability but is not compatible with oscillations.13

  • 2. 

    Single moiety: Moieties correspond to conserved quantities in conservative chemical reaction networks [32]. In a growing cell, every moiety must be fed, else chemical species belonging to the unfed moiety become extinct by dilution [5]. A system exhibiting a single moiety has been preferred in this work because the single nutrient feed reduces the number of parameters in the resulting ODE system. The Brusselator requires two nutrient species, A and B, while all other systems only require A as single nutrient species (the Brusselator also has two low-energy compounds, D and E, while all other systems only have B as a single low-energy compound).14

  • 3. 

    Short siphon: Siphons correspond to extinction sets in chemical reaction networks, that is, to sets of species that can go extinct independently of species outside the set [2]. For some reaction networks, the only siphon is the full set of species. Other networks exhibit siphons that are shorter than the full set of species. A growing cell contains a number of constituents that are only found inside, and not in the outside growth medium. For the nutrient to be available alone in the outside nutrient (i.e., isolated from all other species), it is necessary that the set of all chemical constituents, minus the nutrient, be a siphon. Else, the nutrient in the outside growth medium would spontaneously react into other constituents, independently of the presence of a cell [6]. Furthermore, to avoid this above-mentioned caveat, even in the presence of minute quantities of cell constituents (e.g., those having leaked from a lysed cell), the ultimate chemical insulation can be achieved if the pass reaction (connecting the nutrient with the siphon consisting of all cell constituents) is catalyzed by the membrane or only occurs upon crossing the membrane, as is the case for the phosphoenoltransferase sugar (PTS) active transport system in E. coli [6, 7]. In the PTS system, glucose (glu) is transformed into glucophosphate (g6p) upon crossing the membrane, and glu is not present inside the cytoplasm. The corresponding pass reaction is glu + pep → g6p + pyr, where pep (pyr) designates phosphoenolpyruvate (pyruvate). A system compatible with such a chemical insulation scheme has been preferred in this work.

All chemical reaction networks in Table 1 are realistic in the sense that (i) they are conservative (a set of strictly positive masses can be assigned to each species; in particular, there is no reaction of the kind X → 2X in which mass is created out of nothing), and (ii) a set of Gibbs energies can be assigned to each species in such a way that reactions flow from high to low Gibbs energies.

The system proposed by Wilhelm and Heinrich [42] was chosen in the present work because it meets all the above requirements.

We stress the fact that the list of chemical reaction networks in Table 1 is by no means exhaustive. Other commonly studied CRNs include the Gray-Scott model [19] (equivalent to the Schlögl system without its second reaction, and which may lead to limit cycle behavior within a continuous-flow stirred tank reactor), the Oregonator (an alternative to the Brusselator [17], which also exhibits two moieties), or the CRNs proposed by Schnakenberg [31].

A.2 Identification of Example Sets of Parameters Resulting in Oscillatory Behavior

A.2.1 Making the ODE System Dimensionless to Facilitate Numerical Exploration of the Parameter Space

Making the ODE systems given by Equations 1516 or 1718 dimensionless by changing variables facilitates the numerical exploration of the parameter space. Defining a dimensionless time t′ = αt with α = k3; dimensionless concentrations dj = βcj, j = A, X, Y, Z, or B, with β = k1/k3; and a dimensionless surface-area-to-volume ratio r = γρ with γ = K/k1, we obtain the following dimensionless ODE systems:

  • 1. 
    If Y is the membrane precursor, then
    d˙X=rdAdXdXdYdXldXd˙Y=ψ1dY+ψ2dZrθdYldYd˙Z=dXψ2dZldZd˙B=dXdY+ψ1dYldBr˙r=ϕdYl
    (19)
    where the dimensionless instantaneous growth rate l is given by
    l=rdAdXθdYdA.
    (20)
  • 2. 
    If B is the membrane precursor, then
    d˙X=rdAdXdXdYdXldXd˙Y=ψ1dY+ψ2dZldYd˙Z=dXψ2dZldZd˙B=dXdY+ψ1dYrθdBldBr˙r=ϕdBl
    (21)
    where the dimensionless instantaneous growth rate l is given by
    l=rdAdXθdBdA.
    (22)

The dimensionless parameters are given by
ψ1=k2k3ψ2=k4k3θ=KmeKk1k3ϕ=Kmek1Nme.
(23)

Thus making the ODE system dimensionless reduces the number of parameters from 8 ({ki}i = 1,2,3,4,cA, K, Kme, and Nme) down to 5 (ψ1, ψ2, dA, θ, and ϕ).

A.2.2 Numerical Exploration of the Parameter Space
The differential equation systems 19–20 and 21–22 are too complex to be handled analytically. The determination of the fixed point involves a fourth-degree polynomial equation, and the determination of the Jacobian eigenvalues at this fixed point involves a fifth-degree polynomial equation, involving not only the system parameters (ψ1, ψ2, dA, θ, and ϕ), but also the fixed point itself. There is, however, one simple constraint that can be derived analytically: If Y is chosen as membrane precursor, the following inequality must be verified to ensure the existence of a nondegenerate fixed point:15
ψ1<dAϕ
(24)

If B is chosen as membrane precursor, no such simple constraint could be identified. The inequality (24) helps restrict the parameter space to be explored in the search for sets resulting in oscillatory behavior. It should also be noted that the goal of this study is not to fully characterize these differential equation systems, but simply to identify examples of oscillatory behavior, which will be used to identify qualitative characteristics and constraints that should be applicable to any oscillating cell.

The following procedure was used to identify such examples:

  • 1. 

    B as membrane precursor: All five (ψ1, ψ2, dA, θ, and ϕ) dimensionless parameters were set to unity and let vary by ±1 order of magnitude.

  • 2. 

    Y as membrane precursor: All four (ψ2, dA, θ, and ϕ) dimensionless parameters were set to unity and let vary by ±1 order of magnitude, and for each such choice, the remaining dimensionless parameter ψ1 was let vary by ±1 order of magnitude around an average value two orders of magnitude below the upper bound defined by the inequality (24).16

This resulted in 35 = 243 combinations per choice of membrane precursor. The corresponding trajectories were numerically integrated and analyzed to characterize their asymptotic behavior.

The following initial conditions were adopted: all dimensionless concentrations equal to one-fourth of the outside nutrient concentration; dj (t = 0) = dA/4, j = X, Y, Z, or B (this choice satisfies the balance of osmolarity given by uTc = cA using the notation of Section 2); and the dimensionless surface-area-to-volume ratio equal to unity, r(t = 0) = 1.

Numerical integration was carried out using the scipy.integrate.odeint native solver of the Scipy scientific library of Python (Scipy version 0.17.0, Python version 2.7.11). All trajectories fell in four classes: convergence towards constant strictly positive values (stationary fixed point), divergence of ρ towards infinity (cell collapsing onto itself), convergence of the membrane precursor concentration towards zero (cell ceasing to grow), and oscillatory behavior of all concentrations and of the surface-area-to-volume ratio.

Oscillations were detected and characterized by Fourier transformation of the asymptotic trajectories,17 and by detecting peaks with amplitude higher than 2% of the zeroth harmonic. Tperiod was determined from the first harmonic peak position. The doubling time Tdbl was computed by integration of Equation 3 over one doubling time. An oscillatory behavior was observed for eight parameter combinations if the membrane precursor was Y, and only one parameter combination if the membrane precursor was B. Focusing on Y as membrane precursor, the range of values leading to oscillatory behavior was sharply distributed for all parameters except ψ1. These eight oscillating combinations are such that
dA=10forall8combinationsϕ=0.1forall8combinationsθ=0.1for5combinations,1for3combinationsψ2=0.1for6combinations,1for2combinations.

The same above-described ±1 order-of-magnitude variation was reiterated starting from the above set of most frequent parameters: dA = 10, ϕ = 0.1, θ = 0.1, ψ2 = 0.1, with ψ1 = 1 if B is the membrane precursor, or with ψ1 varying by ±1 order of magnitude around an average value two orders of magnitude below the upper bound defined by the inequality (24) if Y is the membrane precursor.

An oscillatory behavior was observed for 41 (respectively, 11) combinations if Y (B) is the membrane precursor, with Tdbl/Tperiod ranging from 0.20 to 2.34 (from 0.25 to 7.54). The parameter combinations with Tdbl/Tperiod closest to 1 and 2, respectively, were selected as examples. Tdbl/Tperiod was further adjusted to exactly 1 and 2, respectively, by adjusting the single parameter ϕ (which was numerically found to have the largest influence on this ratio). The resulting examples are listed in Table 2.

Table 2. 

Examples of parameter sets resulting in oscillatory behavior, and meeting the synchronization constraint of having the doubling time Tdbl equal to, or a multiple of, the period Tperiod. Two examples were generated, for each of the two choices of the membrane precursor (Y or B).

Membrane precursor Y B 
k1 (M−1 ⋅ s−18.207 1.060 × 101 5.460 1.197 × 101 
k2 (s−18.207 × 10−5 1.060 × 10−6 5.460 × 10−4 1.197 × 10−2 
k3 (s−18.207 × 10−3 1.060 × 10−2 5.460 × 10−3 1.197 × 10−1 
k4 (s−18.207 × 10−3 1.060 × 10−3 5.460 × 10−3 1.197 × 10−2 
cA (M) 1 × 10−2 
K (dm ⋅ M−1 ⋅ s−11 × 10−4 1 × 10−5 1 × 10−5 1 × 10−4 
Kme (dm ⋅ s−11 × 10−8 
Nme (mol ⋅ dm−21 × 10−7 
Tdbl (s) 3.921 × 103 9.051 × 103 1.348 × 103 9.078 × 102 
Tperiod (s) 3.921 × 103 4.526 × 103 1.348 × 103 4.539 × 102 
TdblTperiod 
Label Y1 Y2 B1 B2 
Membrane precursor Y B 
k1 (M−1 ⋅ s−18.207 1.060 × 101 5.460 1.197 × 101 
k2 (s−18.207 × 10−5 1.060 × 10−6 5.460 × 10−4 1.197 × 10−2 
k3 (s−18.207 × 10−3 1.060 × 10−2 5.460 × 10−3 1.197 × 10−1 
k4 (s−18.207 × 10−3 1.060 × 10−3 5.460 × 10−3 1.197 × 10−2 
cA (M) 1 × 10−2 
K (dm ⋅ M−1 ⋅ s−11 × 10−4 1 × 10−5 1 × 10−5 1 × 10−4 
Kme (dm ⋅ s−11 × 10−8 
Nme (mol ⋅ dm−21 × 10−7 
Tdbl (s) 3.921 × 103 9.051 × 103 1.348 × 103 9.078 × 102 
Tperiod (s) 3.921 × 103 4.526 × 103 1.348 × 103 4.539 × 102 
TdblTperiod 
Label Y1 Y2 B1 B2 
A.2.3 Conversion of Dimensionless Parameters to Realistic Physico-chemical Parameters

To convert the five dimensionless parameters (ψ1, ψ2, dA, θ, and ϕ) back to eight realistic physico-chemical parameters ({ki}i=1,2,3,4, cA, K, Kme, and Nme), realistic values are first assigned to three physico-chemical parameters, cA, Nme, and Kme:

  • 1. 

    cA = 1 × 10−2 M. This is equal not only to the nutrient concentration, but also to the total osmolarity. The chosen value is in the range of protocell experimental work [41, 44]. Although the total metabolite concentration in modern evolved cells is significantly higher (around 3 × 10−1 M for E. coli [26]), the present cell only contains four constituents instead of several thousands for modern evolved cells, so that this chosen concentration results in realistic values for metabolite concentrations.

  • 2. 

    Nme = 1 × 10−7 mol ⋅ dm−2. This value is inferred from reported experimental data for E. coli chosen as an arbitrary reference organism. The total cell dry weight (DW) reported in [3] is 300 fg, out of which an approximate 10% corresponds to the membrane [26]. The reported cell volume is 2.3 × 10−3 L ⋅ g−1 (per gram of DW) [4], corresponding to a unit cell volume of 6.9 × 10−16 L. The assumed cell geometry approximately matching this volume is that of a cylinder with a length of 3 μm and a diameter of 0.5 μm, which corresponds to an approximate 5-μm2 membrane surface area. The membrane is primarily constituted of phosphatidylethanolamine with molecular weight 744 g ⋅ mol−1 [20]. Dividing the membrane DW by this molecular weight and by the membrane surface area gives Nme=30fg/744gmol1×5μm21×107moldm2.

  • 3. 

    Kme = 1 × 10−8 dm ⋅ s−1. This value is inferred from a typical doubling time of Tdbl = 3600 s. Integrating Equation 3 over one cell cycle gives log(2)/Tdbl = Kmecavg/Nme, where cavg is the average membrane precursor concentration (averaging over one cell cycle). Making the crude approximation cavgcA/4 to reflect the fact that the cell has four constituents, we have Kme ≈ (log(2)/Tdbl)(Nme/cavg) ≈ 1 × 10−8 dm ⋅ s−1.

With these three physico-chemical parameters cA, Nme, and Kme on the one hand, and the five dimensionless parameters ψ1, ψ2, dA, θ, and ϕ on the other hand, the remaining five physico-chemical parameters can be given by the following relations, which are derived from the equation system (23) and from the definition of β = k1/k3:
k1=KmeϕNmek2=ψ1KmeϕNmecAdAk3=KmeϕNmecAdAk4=ψ2KmeϕNmecAdAK=KmeθdAcA.
(25)

A.3 Constraints on the Choice of t0

A.3.1 η(t) Trajectories

The ηt=v˙/v (see Equation 11) trajectories are shown in Figure 7. Although in principle v(t) might decrease monotonically in the case of mild oscillations, this is not what is seen in the present numerical examples, where η(t) spans both negative and positive values. Such a nonmonotonic variation of v(t) imposes restrictions on the precise starting point t0 of the cell cycle:

  • 1. 

    This starting point should satisfy η(t = 0+) ≤ 0, so that v(t = 0+) ≤ 1. This imposes t0t0,min, where t0,min is a zero of η(t) with negative first derivative. Such t0,min was chosen as time origin in the η(t) trajectories shown in Figure 7.

  • 2. 

    It should also be chosen such that at any time during the trajectory, v(t) ≤ 1. This imposes t0t0,max, where t0,max is defined by t0,mint0,maxηtdt=log2/2n, where n=Tdbl/Tperiod (see proof in Section A.3.2 below). This t0,max is annotated in the η(t) trajectories shown in Figure 7.

There thus exists a range of possible starting points, each associated with its own shape trajectory.

Figure 7. 

Asymptotic trajectories for the logarithmic derivative of the reduced volume, η=v˙v, shown for one doubling time. For all four examples, the time was counted in units of Tdbl, which makes the plots comparable in terms of variation amplitude. The starting point t0 of the cell cycle (when the cell is assumed spherical) must satisfy t0,min < t0 < t0,max (to ensure that the reduced volume remains lower than or equal to unity). t0,min was chosen as the time origin on the shown trajectories. Besides, the spontaneous curvature C0 has not only a lower bound (to ensure cell division when both surface area and volume are doubled), but also has an upper bound when the magnitude of oscillations is so high that t0,int > t0,max (this upper bound ensures that the cell does not divide early during the cell cycle). This is the case for examples Y1, Y2, and B2. When oscillations are milder so that t0,int < t0,max, and when t0 > t0,int, this constraint is alleviated. This is the case for example B1.

Figure 7. 

Asymptotic trajectories for the logarithmic derivative of the reduced volume, η=v˙v, shown for one doubling time. For all four examples, the time was counted in units of Tdbl, which makes the plots comparable in terms of variation amplitude. The starting point t0 of the cell cycle (when the cell is assumed spherical) must satisfy t0,min < t0 < t0,max (to ensure that the reduced volume remains lower than or equal to unity). t0,min was chosen as the time origin on the shown trajectories. Besides, the spontaneous curvature C0 has not only a lower bound (to ensure cell division when both surface area and volume are doubled), but also has an upper bound when the magnitude of oscillations is so high that t0,int > t0,max (this upper bound ensures that the cell does not divide early during the cell cycle). This is the case for examples Y1, Y2, and B2. When oscillations are milder so that t0,int < t0,max, and when t0 > t0,int, this constraint is alleviated. This is the case for example B1.

A.3.2 Proof That v(t) ≤ 1 ⇒ t0t0,max, with t0,max Such That t0,mint0,maxηtdt=log22n

The proof proceeds in two steps:

  • 1. 

    For any time t,vt+Tperiod=21/2n×vt.

    Proof : We have tt+Tdblηtdt=logvt+Tdbl/vt=log1/2=log2/2 because A is doubled every Tdbl and ρ is periodic. Besides, tt+Tdblηtdt=n×tt+Tperiodηtdt because η(t) is periodic with period Tperiod and Tdbl = nTperiod, so that tt+Tperiodηtdt=log2/2n, or equivalently vt+Tperiod=21/2n×vt. In particular, vt0,max=21/2n×vt0,min=vt0,min+Tperiod.

  • 2. 

    If t0 > t0,max, then v(t) exceeds unity at time t = t0,min + Tperiodduring the cell cycle.

    Proof: If t0 > t0,max, then v(t0) < v(t0,max) because η(t = t0+) ≤ 0, and v(t0,min + Tperiod) = v(t0,max) > v(t0) = 1.

A.4 Possible Upper Bound for C0 Depending on the Amplitude of Oscillations

If, at some critical point before cycle end, tc < t0 + Tdbl, the reduced volume reaches vtc=12, the double sphere should have a higher bending energy than the elongated rod, so that the cell may not divide into two smaller cells (smaller than the original cell). Whether v(t) may actually reach 12 before cycle end or not depends on the amplitude of η(t) oscillations. Over the course of one period, v(t) reaches a local minimum at the zero of η(t) having a positive first derivative, denoted t1, and over the course of one cycle, v(t) reaches its last and lowest local minimum at t = t1 + (n − 1)Tperiod. Defining t0,int by t0,intt1ηtdt=log22n, we have the following properties, which are proved in turn below:

  • 1. 

    If t0 > t0,int, then vt=t1+n1Tperiod>12 and v(t) does not reach 12 before cycle end. Conversely, if t0 < t0,int, v(t) reaches 12 before cycle end.

  • 2. 

    If log2n<t0,mint1ηtdt<0, then t0,int < t0,max, and v(t) does not reach 12 if t0 is chosen such that t0,int < t0 < t0,max. Conversely, if t0,mint1ηtdt<log2n, then v(t) reaches 12 for any choice of t0, t0,min < t0 < t0,max.

A.4.1 Proof That If t0 > t0,int, Then v (t = t1 + (n − 1)Tperiod)>12

We recall that t0,int is defined by t0,intt1ηtdt=log22n. If t0 > t0,int, then vt1=exp(t0t1ηtdt)>exp(t0,intt1ηtdt)=212n. If n > 1, v(t) further decreases over the course of the remaining n − 1 periods, to reach its lowest local minimum vt1+n1Tperiod=vt1×2n12n>212n×2n12n=12. The reduced volume cannot reach 12 before cycle end.

A.4.2 Proof That t0,int < t0,maxlog2n<t0,mint1ηtdt<0

We first note that t0,int<t0,maxt0,intt0,maxηtdt<0, because η(t) < 0 for t0,min < t < t1 and both t0,int and t0,max fall in that interval. Second, we decompose t0,intt0,maxηtdt=t0,intt1ηtdt+t1t0,minηtdt+t0,mint0,maxηtdt=log22n+t1t0,minηtdtlog22n=log2n+t1t0,minηtdt=log2nt0,mint1ηtdt. This gives t0,intt0,maxηtdt<0log2n<t0,mint1ηtdt.

Author notes

Laboratoire de Biochimie, École Polytechnique, CNRS, Université Paris Saclay Palaiseau 91128, France. E-mail: plateau@bioc.polytechnique.fr