On the Relation between Chemical Oscillations and Self-Replication

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.

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 c_A remains constant. The nutrient A is imported into the cell by a surfacic process characterized by a reaction vector F_in. The membrane precursor is incorporated into a self-assembled membrane by a surfacic process characterized by a reaction vector F_out.

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.

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 T_dbl, that is, $A$(t + T_dbl) = 2$A$(t) and $V$(t + T_dbl) = 2$V$(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)/T_dbl. 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 T_dbl.² 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 T_dbl.

It should be noted that the doubling time T_dbl 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 T_period should match T_dbl. But the above reasoning shows that such a match is a requirement for a self-replicating system.

To be more precise, the shortest period T_period common to ζ_$A$(t), ζ_$V$(t), ρ(t), and concentrations need not necessarily be T_dbl, but can also be one of its submultiples, that is, T_dbl = nT_period, 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.³

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$⁠.⁵

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.⁶

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

• 1. 

  Elongated rod: y_el = 0 and f(v, z_el) = 0.

• 2. 

  Most constricted shape:

  • (a) 

    y_cs = 1 and f(v, z_cs) = 1 if $v≤12$ (two elongated rods attached by an infinitesimal neck).

  • (b) 

    y_cs = z_cs − 1 and f(v, z_cs) = z_cs − 1 if $v>12$ (two overlapping spheres).

• 1. 

  a(t₀) = 1, v(t₀) = 1, y(t₀) = 0, z(t₀) = 1;

• 2. 

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

C₀ 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.

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 C₀ values above the threshold value.

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

In the range x_thr < C₀R₀ ≤ 1.05 × x_thr, the shape trajectory is nearly identical to that of C₀R₀ = 1.05 × x_thr (Figure 3a), with constriction only occurring towards the end of the cell cycle. In the range C₀R₀ ≥ 1.35 × x_thr, it is nearly identical to that of C₀R₀ = 1.35 × x_thr (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 C₀R₀ ≈ 1.82 × x_thr.

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 C₀R₀ > x_thr. On the contrary, if C₀R₀ < x_thr, 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).

This system is conservative (all masses are equal), and the directionality of each reaction corresponds to the following ranking of Gibbs energies: g_A > g_X > g_Z > g_Y > g_B. 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 F_in has as single nonzero component F_in,X = Kc_Ac_X, where K is a surfacic kinetic constant and c_A and c_X are the concentrations of A and X respectively. A is only present in the outside growth medium (assumed sufficiently large so that c_A 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 K_me per unit area, so that F_out has as single nonzero component F_out,Y = K_mec_Y if Y is the membrane precursor (or F_out,B = K_mec_B if B is the membrane precursor).¹⁰

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

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=ρKcAcX−k1cXcY−k3cX−λcXc˙Y=−k2cY+k4cZ−ρKmecY−λcYc˙Z=k3cX−k4cZ−λcZc˙B=k1cXcY+k2cY−λcBρ˙ρ=KmecYNme−λ$

  (15)

  where the growth rate λ is given by

  $λ=ρKcAcX−KmecYcA.$

  (16)
• 2. 
  If B is the membrane precursor,

  $c˙X=ρKcAcX−k1cXcY−k3cX−λcXc˙Y=−k2cY+k4cZ−λcYc˙Z=k3cX−k4cZ−λcZc˙B=k1cXcY+k2cY−ρKmecB−λcBρ˙ρ=KmecBNme−λ$

  (17)

  where the growth rate λ is given by

  $λ=ρKcAcX−KmecBcA.$

  (18)

In the above equations, c_A (c_X, c_Y, c_Z, c_B) is the concentration of A (of X, Y, Z, B), and k₁ (k₂, k₃, k₄) is the kinetic coefficient for the reaction X + Y → Y + B (the reaction Y → B, X → Z, Z → Y).

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.

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.

Contrary to the stationary case, the trajectories v(t) and a(t) depend on the choice t₀ 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 t₀ 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 C₀:

• 1. 

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

• 2. 

  Besides, if at some critical point before cycle end, say t_c < t₀ + T_dbl, 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: C₀R_c < x_thr, where R_c is the sphere radius for the double sphere having a(t_c) < 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 t_c constrains C₀ and R(t₀) to be closely related, as C₀R(t₀) 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 t₀ and C₀ 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 t₀ close to its upper bound t_0,max. This t₀ 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 C₀ close to its upper bound, while for example B1 (for which C₀ is only constrained by a lower bound), it was found for C₀ 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.

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.

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 x_thr ≈ 1.18 for C₀R(t₀) depends on this choice, as illustrated by the alternative value for a different shape comparison in footnote ⁹. But the existence of such a minimum is expected to be quite general, and the chosen predefined family of shapes helps understand why. C₀ 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 C₀R(t₀) 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 C₀R(t₀) would be applicable, as is the case for examples Y1, Y2, and B2, but not for B1) is still an open question.¹² 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).

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 C₀R(t₀) 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 C₀ and R(t₀) suggests a coevolution of the membrane and of the embedded metabolism. This is because C₀ depends on membrane structural properties, while R(t₀) depends on the biochemical and transmembrane kinetics.

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 C₀R(t₀ + pT_dbl)—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 C₀R(t₀ + pT_dbl) 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.

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.

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.

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.

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.

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.¹³

• 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).¹⁴

• 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].

Making the ODE systems given by Equations 15–16 or 17–18 dimensionless by changing variables facilitates the numerical exploration of the parameter space. Defining a dimensionless time t′ = αt with α = k₃; dimensionless concentrations d_j = βc_j, j = A, X, Y, Z, or B, with β = k₁/k₃; and a dimensionless surface-area-to-volume ratio r = γρ with γ = K/k₁, we obtain the following dimensionless ODE systems:

• 1. 
  If Y is the membrane precursor, then

  $d˙X=rdAdX−dXdY−dX−ldXd˙Y=−ψ1dY+ψ2dZ−rθdY−ldYd˙Z=dX−ψ2dZ−ldZd˙B=dXdY+ψ1dY−ldBr˙r=ϕdY−l$

  (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=rdAdX−dXdY−dX−ldXd˙Y=−ψ1dY+ψ2dZ−ldYd˙Z=dX−ψ2dZ−ldZd˙B=dXdY+ψ1dY−rθdB−ldBr˙r=ϕdB−l$

  (21)

  where the dimensionless instantaneous growth rate l is given by

  $l=rdAdX−θdBdA.$

  (22)

Thus making the ODE system dimensionless reduces the number of parameters from 8 ({k_i}_{i = 1,2,3,4,}c_A, K, K_me, and N_me) down to 5 (ψ₁, ψ₂, d_A, θ, and ϕ).

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 (ψ₁, ψ₂, d_A, θ, and ϕ) dimensionless parameters were set to unity and let vary by ±1 order of magnitude.

• 2. 

  Y as membrane precursor: All four (ψ₂, d_A, θ, and ϕ) dimensionless parameters were set to unity and let vary by ±1 order of magnitude, and for each such choice, the remaining dimensionless parameter ψ₁ 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).¹⁶

This resulted in 3⁵ = 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; d_j (t = 0) = d_A/4, j = X, Y, Z, or B (this choice satisfies the balance of osmolarity given by u^Tc = c_A 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.

The same above-described ±1 order-of-magnitude variation was reiterated starting from the above set of most frequent parameters: d_A = 10, ϕ = 0.1, θ = 0.1, ψ₂ = 0.1, with ψ₁ = 1 if B is the membrane precursor, or with ψ₁ 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.

To convert the five dimensionless parameters (ψ₁, ψ₂, d_A, θ, and ϕ) back to eight realistic physico-chemical parameters ({k_i}_i=1,2,3,4, c_A, K, K_me, and N_me), realistic values are first assigned to three physico-chemical parameters, c_A, N_me, and K_me:

• 1. 

  c_A = 1 × 10⁻² 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⁻¹ 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. 

  N_me = 1 × 10⁻⁷ mol ⋅ dm⁻². 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⁻³ L ⋅ g⁻¹ (per gram of DW) [4], corresponding to a unit cell volume of 6.9 × 10⁻¹⁶ 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-μm² membrane surface area. The membrane is primarily constituted of phosphatidylethanolamine with molecular weight 744 g ⋅ mol⁻¹ [20]. Dividing the membrane DW by this molecular weight and by the membrane surface area gives $Nme=30fg/744g⋅mol−1×5μm2≈1×10−7mol⋅dm−2$⁠.

• 3. 

  K_me = 1 × 10⁻⁸ dm ⋅ s⁻¹. This value is inferred from a typical doubling time of T_dbl = 3600 s. Integrating Equation 3 over one cell cycle gives log(2)/T_dbl = K_mec_avg/N_me, where c_avg is the average membrane precursor concentration (averaging over one cell cycle). Making the crude approximation c_avg ≈ c_A/4 to reflect the fact that the cell has four constituents, we have K_me ≈ (log(2)/T_dbl)(N_me/c_avg) ≈ 1 × 10⁻⁸ dm ⋅ s⁻¹.

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 t₀ of the cell cycle:

• 1. 

  This starting point should satisfy η(t = 0⁺) ≤ 0, so that v(t = 0⁺) ≤ 1. This imposes t₀ ≥ t_0,min, where t_0,min is a zero of η(t) with negative first derivative. Such t_0,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 t₀ ≤ t_0,max, where t_0,max is defined by $∫t0,mint0,maxηtdt=−log2/2n$⁠, where $n=Tdbl/Tperiod$ (see proof in Section A.3.2 below). This t_0,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.

The proof proceeds in two steps:

• 1. 

  For any time t,$vt+Tperiod=2−1/2n×vt$⁠.

  Proof : We have $∫tt+Tdblηt′dt′=logvt+Tdbl/vt=log1/2=−log2/2$ because $A$ is doubled every T_dbl and ρ is periodic. Besides, $∫tt+Tdblηt′dt′=n×∫tt+Tperiodηt′dt′$ because η(t) is periodic with period T_period and T_dbl = nT_period, so that $∫tt+Tperiodηt′dt′=−log2/2n$⁠, or equivalently $vt+Tperiod=2−1/2n×vt$⁠. In particular, $vt0,max=2−1/2n×vt0,min=vt0,min+Tperiod$⁠.

• 2. 

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

  Proof: If t₀ > t_0,max, then v(t₀) < v(t_0,max) because η(t = t₀⁺) ≤ 0, and v(t_0,min + T_period) = v(t_0,max) > v(t₀) = 1.

If, at some critical point before cycle end, t_c < t₀ + T_dbl, 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 t₁, and over the course of one cycle, v(t) reaches its last and lowest local minimum at t = t₁ + (n − 1)T_period. Defining t_0,int by $∫t0,intt1ηtdt=−log22n$⁠, we have the following properties, which are proved in turn below:

• 1. 

  If t₀ > t_0,int, then $vt=t1+n−1Tperiod>12$ and v(t) does not reach $12$ before cycle end. Conversely, if t₀ < t_0,int, v(t) reaches $12$ before cycle end.

• 2. 

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

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

We first note that $t0,int<t0,max⇔∫t0,intt0,maxηtdt<0$⁠, because η(t) < 0 for t_0,min < t < t₁ and both t_0,int and t_0,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ηtdt−log22n=−log2n+∫t1t0,minηtdt=−log2n−∫t0,mint1ηtdt.$ This gives $∫t0,intt0,maxηtdt<0⇔−log2n<∫t0,mint1ηtdt$⁠.