## 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 *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}.

**is the concentration vector inside the cell,**

*c***is the CRN stoichiometry matrix,**

*S***=**

*f***(**

*f***) is the reaction rate vector, ρ = $A$/$V$ is the surface area-to-volume ratio, and λ = $V\u02d9$/$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**

*c**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
*u*^{T}=*c**c*_{A}, whereis the osmolarity vector such that*u**u*_{i}is the number of particles in solution when dissolving one molecule of the chemical species indexed by*i*. The inside osmolarity is thus constant,*u*^{T}= 0. Multiplying each side of Equation 1 by*ċ**u*^{T}and rearranging terms, gives$\lambda =uTSf+\rho Fin\u2212FoutuTc.$(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 $\rho \u02d9/\rho =A\u02d9/A\u2212V\u02d9/V=A\u02d9/A\u2212\lambda $. 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
*F*_{out},^{1}which results in a relative rate of increase for $A$ given bywhere$A\u02d9A=FoutNme$(3)*N*_{me}is the number of molecules per unit area in the self-assembled membrane. This results in the following expression for the logarithmic derivative of ρ:$\rho \u02d9\rho =FoutNme\u2212\lambda .$(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 *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(−λ_{avg}*t*)$A$(*t*) and ζ_{$V$}(*t*) = exp(−λ_{avg}*t*)$V$(*t*), we see that ζ_{$A$}(*t*) and ζ_{$V$}(*t*) are periodic with period *T*_{dbl}.^{2} Then $\rho 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\u02d9/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.^{3}

## 3 From Asymptotic Trajectories to Cellular Shape

### 3.1 Spontaneous-Curvature Model, with Predefined Family of Shapes

_{G}are the bending rigidity and the Gaussian bending rigidity, respectively,

*C*

_{1}and

*C*

_{2}are the principal curvatures, and

*C*

_{0}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<\omega <\pi 2$, two overlapping spheres $L=0,0<\omega <\pi 2$, or a double sphere $L=0,\omega =\pi 2$.^{5}

Starting from an assumed spherical shape (with radius $Rt0=3\rho 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

*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*≤ 1 with

*v*= 1 for a sphere.

*a*is defined by

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

*R*,

*L*, and ω, it proves convenient to use the following set of transformed parameters: $y=sin\omega ,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

*a*, the values of

*R*and

*z*are related by

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\u226412$ (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).

- (a)

*t*=

*t*

_{0}) and at the end of the cell cycle (

*t*=

*t*

_{0}+

*T*

_{dbl}) are given by:

- 1.
*a*(*t*_{0}) = 1,*v*(*t*_{0}) = 1,*y*(*t*_{0}) = 0,*z*(*t*_{0}) = 1; - 2.
$at0+Tdbl=2,vt0+Tdbl=12,yt0+Tdbl=1,zt0+Tdbl=2$.

*y*vary between the two extremes,

*y*

_{el}= 0 ≤

*y*≤

*y*

_{cs}, 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

*E*

_{b}.

^{7}

### 3.3 Determining the Reduced Volume and Normalized Area Trajectories

*v*, logarithmic differentiation of Equation 6 gives

**(**

*c**t*) and ρ(

*t*), then η(

*t*) is also periodic with the same period

*T*

_{period}.

^{8}

*v*(

*t*) is given by integration of η:

### 3.4 The Spontaneous Curvature *C*_{0} Should Be Sufficiently Large

*C*

_{0}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 (

*E*

_{b}= 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:

*C*

_{0}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

*C*

_{0}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.

*C*_{0} reflects a structural property of the membrane, and $Rt0=3/\rho 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\u2212\lambda t/2$ and *a*(*t*) = exp(λ*t*), where $\lambda =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*_{0} values above the threshold value.

The shape pattern only depends on the product *C*_{0}*R*_{0}, where $R0=3/\rho 0$ with ρ_{0} being the stationary surface-area-to-volume ratio. Only the size and time scales depend on the specific stationary asymptotic values: *R*_{0} 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*_{0}*R*_{0} ≤ 1.05 × *x*_{thr}, the shape trajectory is nearly identical to that of *C*_{0}*R*_{0} = 1.05 × *x*_{thr} (Figure 3a), with constriction only occurring towards the end of the cell cycle. In the range *C*_{0}*R*_{0} ≥ 1.35 × *x*_{thr}, it is nearly identical to that of *C*_{0}*R*_{0} = 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*_{0}*R*_{0} ≈ 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*_{0}*R*_{0} > *x*_{thr}. On the contrary, if *C*_{0}*R*_{0} < *x*_{thr}, then starting from a spherical shape with radius $R0=3/\rho 0$, the cell keeps elongating without constriction towards a filament shape with radius $R0=2/\rho 0$ (infinitely long cylinder).

### 4.2 Oscillating Case

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* → 2*X* was taken as the nutrient import reaction. The associated reaction vector *F*_{in} has as single nonzero component *F*_{in,X} = *Kc*_{A}*c*_{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*_{me}*c*_{Y} if *Y* is the membrane precursor (or *F*_{out,B} = *K*_{me}*c*_{B} if *B* is the membrane precursor).^{10}

All chemical species are assumed to contribute to one particle when dissolved in solution, *u*_{i} = 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,where the growth rate λ is given by$c\u02d9X=\rho KcAcX\u2212k1cXcY\u2212k3cX\u2212\lambda cXc\u02d9Y=\u2212k2cY+k4cZ\u2212\rho KmecY\u2212\lambda cYc\u02d9Z=k3cX\u2212k4cZ\u2212\lambda cZc\u02d9B=k1cXcY+k2cY\u2212\lambda cB\rho \u02d9\rho =KmecYNme\u2212\lambda $(15)$\lambda =\rho KcAcX\u2212KmecYcA.$(16) - 2. If
*B*is the membrane precursor,where the growth rate λ is given by$c\u02d9X=\rho KcAcX\u2212k1cXcY\u2212k3cX\u2212\lambda cXc\u02d9Y=\u2212k2cY+k4cZ\u2212\lambda cYc\u02d9Z=k3cX\u2212k4cZ\u2212\lambda cZc\u02d9B=k1cXcY+k2cY\u2212\rho KmecB\u2212\lambda cB\rho \u02d9\rho =KmecBNme\u2212\lambda $(17)$\lambda =\rho KcAcX\u2212KmecBcA.$(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*_{1} (*k*_{2}, *k*_{3}, *k*_{4}) 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/\rho 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*_{0} 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*_{0} 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*_{0}:

- 1.
Appendix A.3 explains why

*t*_{0}must satisfy*t*_{0,min}<*t*_{0}<*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*_{0}+*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*_{0}*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*_{0}and*R*(*t*_{0}) to be closely related, as*C*_{0}*R*(*t*_{0}) 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*_{0} and *C*_{0} 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*_{0} close to its upper bound *t*_{0,max}. This *t*_{0} 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*_{0} close to its upper bound, while for example B1 (for which *C*_{0} is only constrained by a lower bound), it was found for *C*_{0} 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 *x*_{thr} ≈ 1.18 for *C*_{0}*R*(*t*_{0}) 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. *C*_{0} 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*_{0}*R*(*t*_{0}) 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*_{0}*R*(*t*_{0}) 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 *C*_{0}*R*(*t*_{0}) 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*_{0} and *R*(*t*_{0}) suggests a coevolution of the membrane and of the embedded metabolism. This is because *C*_{0} depends on membrane structural properties, while *R*(*t*_{0}) 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

*C*_{0}*R*(*t*_{0}+*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*_{0}*R*(*t*_{0}+*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.

### 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

Here, the scalar *F*_{out} is the component of the vector *F*_{out} along the membrane precursor.

$\zeta At+Tdbl=exp\u2212\lambda avgtexp\u2212\lambda avgTdblAt+Tdbl=exp\u2212\lambda avgt\xd712\xd7At+Tdbl=exp\u2212\lambda avgtAt=\zeta At$, and similarly $\zeta Vt+Tdbl=\zeta Vt$.

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 *T*_{period} 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 1–3, and all these quantities oscillate with the same shortest period *T*_{period}.

See in particular Figure 18 in [33]. On a (*v*, *c*_{0}) diagram, where *v* is the reduced volume and *c*_{0} is the scaled spontaneous curvature defined as the product of the spontaneous curvature *C*_{0} 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.

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.

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.

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

The logarithmic derivative of $A$ = ζ_{$A$} × exp(λ_{avg}*t*) is $A\u02d9/A=\zeta \u02d9A/\zeta A+\lambda avg$. Besides, $\rho \u02d9/\rho =\zeta \u02d9A/\zeta A\u2212\zeta \u02d9V/\zeta V$. This gives $\eta =\zeta \u02d9V/\zeta V\u22123/2\zeta \u02d9A/\zeta A\u22121/2\lambda avg$, which is periodic because the derivative of a periodic function is periodic.

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*(*t*_{0})) connected by infinitely narrow necks is more stable than the infinite cylinder having the same area and the same volume if $C0Rt0>7/4$.

For this example, the osmolarity vector ** u** is parallel to the mass vector

**, which results in**

*m*

*u*^{T}

**= 0 because the CRN is conservative. This simplifies the expression of Equation 2 for the growth rate.**

*S*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 $\u22121/2$). By a high amplitude we mean that the area under the negative lobe is larger (in absolute value) than $2\xd71/2=2$ (see Appendix A.4).

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.

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

Setting all derivatives in the equation system (19) to zero, the last equation gives *l* = ϕ*d*_{Y}. Feeding this into the equation before last gives $dB=dX+\psi 1/\varphi >\psi 1/\varphi $. Feeding this into the constant osmolarity gives $dA=dX+dY+dZ+dB>dB>\psi 1/\varphi $, or equivalently, the inequality (24).

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

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

## References

*Escherichia coli*

*A*+ 2

*B*→ 3

*B*;

*B*→

*C*

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

System name | Chemical reactions | Limit cycle | Single moiety | Short siphon | Reference |

Lotka-Volterra | A + X → 2X | No | Yes | Yes | [21] |

X + Y → 2Y | |||||

Y → B | |||||

Schlögl | A + 2X → 3X | No (only bistability) | Yes | Yes | [30] |

3X → 2X + B | |||||

X → B | |||||

Brusselator | A → X | Yes | No (2 moieties) | No (a siphon containing X also contains A) | [28] |

2X + Y → 3X | |||||

B + X → Y + D | |||||

X → E | |||||

Wilhelm | A + X → 2X | Yes | Yes | Yes | [42] |

X + Y → Y + B | |||||

Y → B | |||||

X → Z | |||||

Z → Y |

System name | Chemical reactions | Limit cycle | Single moiety | Short siphon | Reference |

Lotka-Volterra | A + X → 2X | No | Yes | Yes | [21] |

X + Y → 2Y | |||||

Y → B | |||||

Schlögl | A + 2X → 3X | No (only bistability) | Yes | Yes | [30] |

3X → 2X + B | |||||

X → B | |||||

Brusselator | A → X | Yes | No (2 moieties) | No (a siphon containing X also contains A) | [28] |

2X + Y → 3X | |||||

B + X → Y + D | |||||

X → E | |||||

Wilhelm | A + X → 2X | Yes | Yes | Yes | [42] |

X + Y → Y + B | |||||

Y → B | |||||

X → Z | |||||

Z → Y |

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* → 2*X* 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 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*_{3}; dimensionless concentrations *d*_{j} = β*c*_{j}, *j* = *A*, *X*, *Y*, *Z*, or *B*, with β = *k*_{1}/*k*_{3}; and a dimensionless surface-area-to-volume ratio *r* = γρ with γ = *K*/*k*_{1}, we obtain the following dimensionless ODE systems:

- 1. If
*Y*is the membrane precursor, thenwhere the dimensionless instantaneous growth rate$d\u02d9X=rdAdX\u2212dXdY\u2212dX\u2212ldXd\u02d9Y=\u2212\psi 1dY+\psi 2dZ\u2212r\theta dY\u2212ldYd\u02d9Z=dX\u2212\psi 2dZ\u2212ldZd\u02d9B=dXdY+\psi 1dY\u2212ldBr\u02d9r=\varphi dY\u2212l$(19)*l*is given by$l=rdAdX\u2212\theta dYdA.$(20) - 2. If
*B*is the membrane precursor, thenwhere the dimensionless instantaneous growth rate$d\u02d9X=rdAdX\u2212dXdY\u2212dX\u2212ldXd\u02d9Y=\u2212\psi 1dY+\psi 2dZ\u2212ldYd\u02d9Z=dX\u2212\psi 2dZ\u2212ldZd\u02d9B=dXdY+\psi 1dY\u2212r\theta dB\u2212ldBr\u02d9r=\varphi dB\u2212l$(21)*l*is given by$l=rdAdX\u2212\theta 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 (ψ_{1}, ψ_{2}, *d*_{A}, θ, and ϕ).

##### A.2.2 Numerical Exploration of the Parameter Space

_{1}, ψ

_{2},

*d*

_{A}, θ, 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}

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},*d*_{A}, θ, and ϕ) dimensionless parameters were set to unity and let vary by ±1 order of magnitude. - 2.
*Y*as membrane precursor: All four (ψ_{2},*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 ψ_{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 3^{5} = 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*^{T}** c** =

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

^{17}and by detecting peaks with amplitude higher than 2% of the zeroth harmonic.

*T*

_{period}was determined from the first harmonic peak position. The doubling time

*T*

_{dbl}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

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, ψ_{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.

Membrane precursor | Y | B | ||

k_{1} (M^{−1} ⋅ s^{−1}) | 8.207 | 1.060 × 10^{1} | 5.460 | 1.197 × 10^{1} |

k_{2} (s^{−1}) | 8.207 × 10^{−5} | 1.060 × 10^{−6} | 5.460 × 10^{−4} | 1.197 × 10^{−2} |

k_{3} (s^{−1}) | 8.207 × 10^{−3} | 1.060 × 10^{−2} | 5.460 × 10^{−3} | 1.197 × 10^{−1} |

k_{4} (s^{−1}) | 8.207 × 10^{−3} | 1.060 × 10^{−3} | 5.460 × 10^{−3} | 1.197 × 10^{−2} |

c_{A} (M) | 1 × 10^{−2} | |||

K (dm ⋅ M^{−1} ⋅ s^{−1}) | 1 × 10^{−4} | 1 × 10^{−5} | 1 × 10^{−5} | 1 × 10^{−4} |

K_{me} (dm ⋅ s^{−1}) | 1 × 10^{−8} | |||

N_{me} (mol ⋅ dm^{−2}) | 1 × 10^{−7} | |||

T_{dbl} (s) | 3.921 × 10^{3} | 9.051 × 10^{3} | 1.348 × 10^{3} | 9.078 × 10^{2} |

T_{period} (s) | 3.921 × 10^{3} | 4.526 × 10^{3} | 1.348 × 10^{3} | 4.539 × 10^{2} |

$TdblTperiod$ | 1 | 2 | 1 | 2 |

Label | Y1 | Y2 | B1 | B2 |

Membrane precursor | Y | B | ||

k_{1} (M^{−1} ⋅ s^{−1}) | 8.207 | 1.060 × 10^{1} | 5.460 | 1.197 × 10^{1} |

k_{2} (s^{−1}) | 8.207 × 10^{−5} | 1.060 × 10^{−6} | 5.460 × 10^{−4} | 1.197 × 10^{−2} |

k_{3} (s^{−1}) | 8.207 × 10^{−3} | 1.060 × 10^{−2} | 5.460 × 10^{−3} | 1.197 × 10^{−1} |

k_{4} (s^{−1}) | 8.207 × 10^{−3} | 1.060 × 10^{−3} | 5.460 × 10^{−3} | 1.197 × 10^{−2} |

c_{A} (M) | 1 × 10^{−2} | |||

K (dm ⋅ M^{−1} ⋅ s^{−1}) | 1 × 10^{−4} | 1 × 10^{−5} | 1 × 10^{−5} | 1 × 10^{−4} |

K_{me} (dm ⋅ s^{−1}) | 1 × 10^{−8} | |||

N_{me} (mol ⋅ dm^{−2}) | 1 × 10^{−7} | |||

T_{dbl} (s) | 3.921 × 10^{3} | 9.051 × 10^{3} | 1.348 × 10^{3} | 9.078 × 10^{2} |

T_{period} (s) | 3.921 × 10^{3} | 4.526 × 10^{3} | 1.348 × 10^{3} | 4.539 × 10^{2} |

$TdblTperiod$ | 1 | 2 | 1 | 2 |

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}, *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^{−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.
*N*_{me}= 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-μm^{2}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/744g\u22c5mol\u22121\xd75\mu m2\u22481\xd710\u22127mol\u22c5dm\u22122$. - 3.
*K*_{me}= 1 × 10^{−8}dm ⋅ s^{−1}. 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*_{me}*c*_{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^{−8}dm ⋅ s^{−1}.

*c*

_{A},

*N*

_{me}, and

*K*

_{me}on the one hand, and the five dimensionless parameters ψ

_{1}, ψ

_{2},

*d*

_{A}, θ, 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 β =

*k*

_{1}/

*k*

_{3}:

#### A.3 Constraints on the Choice of *t*_{0}

##### A.3.1 η(*t*) Trajectories

The $\eta t=v\u02d9/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*_{0} of the cell cycle:

- 1.
This starting point should satisfy η(

*t*= 0^{+}) ≤ 0, so that*v*(*t*= 0^{+}) ≤ 1. This imposes*t*_{0}≥*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*_{0}≤*t*_{0,max}, where*t*_{0,max}is defined by $\u222bt0,mint0,max\eta tdt=\u2212log2/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.

##### A.3.2 Proof That *v*(*t*) ≤ 1 ⇒ *t*_{0} ≤ *t*_{0,max}, with *t*_{0,max} Such That $\u222bt0,mint0,max\eta tdt=\u2212log22n$

The proof proceeds in two steps:

- 1.
*For any time t,*$vt+Tperiod=2\u22121/2n\xd7vt$.Proof : We have $\u222btt+Tdbl\eta t\u2032dt\u2032=logvt+Tdbl/vt=log1/2=\u2212log2/2$ because $A$ is doubled every

*T*_{dbl}and ρ is periodic. Besides, $\u222btt+Tdbl\eta t\u2032dt\u2032=n\xd7\u222btt+Tperiod\eta t\u2032dt\u2032$ because η(*t*) is periodic with period*T*_{period}and*T*_{dbl}=*nT*_{period}, so that $\u222btt+Tperiod\eta t\u2032dt\u2032=\u2212log2/2n$, or equivalently $vt+Tperiod=2\u22121/2n\xd7vt$. In particular, $vt0,max=2\u22121/2n\xd7vt0,min=vt0,min+Tperiod$. - 2.
*If t*_{0}>*t*_{0,max}*, then v*(*t*)*exceeds unity at time t*=*t*_{0,min}+*T*_{period}*during the cell cycle.*Proof: If

*t*_{0}>*t*_{0,max}, then*v*(*t*_{0}) <*v*(*t*_{0,max}) because η(*t*=*t*_{0}^{+}) ≤ 0, and*v*(*t*_{0,min}+*T*_{period}) =*v*(*t*_{0,max}) >*v*(*t*_{0}) = 1.

#### A.4 Possible Upper Bound for *C*_{0} Depending on the Amplitude of Oscillations

If, at some critical point before cycle end, *t*_{c} < *t*_{0} + *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*_{1}, and over the course of one cycle, *v*(*t*) reaches its last and lowest local minimum at *t* = *t*_{1} + (*n* − 1)*T*_{period}. Defining *t*_{0,int} by $\u222bt0,intt1\eta tdt=\u2212log22n$, we have the following properties, which are proved in turn below:

- 1.
If

*t*_{0}>*t*_{0,int}, then $vt=t1+n\u22121Tperiod>12$ and*v*(*t*) does not reach $12$ before cycle end. Conversely, if*t*_{0}<*t*_{0,int},*v*(*t*) reaches $12$ before cycle end. - 2.
If $\u2212log2n<\u222bt0,mint1\eta tdt<0$, then

*t*_{0,int}<*t*_{0,max}, and*v*(*t*) does not reach $12$ if*t*_{0}is chosen such that*t*_{0,int}<*t*_{0}<*t*_{0,max}. Conversely, if $\u222bt0,mint1\eta tdt<\u2212log2n$, then*v*(*t*) reaches $12$ for any choice of*t*_{0},*t*_{0,min}<*t*_{0}<*t*_{0,max}.

##### A.4.1 Proof That If *t*_{0} > *t*_{0,int}, Then *v* (*t* = *t*_{1} + (n − 1)T_{period})>$12$

We recall that *t*_{0,int} is defined by $\u222bt0,intt1\eta tdt=\u2212log22n$. If *t*_{0} > *t*_{0,int}, then $vt1=exp(\u222bt0t1\eta tdt)>exp(\u222bt0,intt1\eta tdt)=2\u221212n$. If *n* > 1, *v*(*t*) further decreases over the course of the remaining *n* − 1 periods, to reach its lowest local minimum $vt1+n\u22121Tperiod=vt1\xd72\u2212n\u221212n>2\u221212n\xd72\u2212n\u221212n=12$. The reduced volume cannot reach $12$ before cycle end.

##### A.4.2 Proof That *t*_{0,int} < *t*_{0,max} ⇔ $\u2212log2n<\u222bt0,mint1\eta tdt<0$

We first note that $t0,int<t0,max\u21d4\u222bt0,intt0,max\eta tdt<0$, because η(*t*) < 0 for *t*_{0,min} < *t* < *t*_{1} and both *t*_{0,int} and *t*_{0,max} fall in that interval. Second, we decompose $\u222bt0,intt0,max\eta tdt=\u222bt0,intt1\eta tdt+\u222bt1t0,min\eta tdt+\u222bt0,mint0,max\eta tdt=\u2212log22n+\u222bt1t0,min\eta tdt\u2212log22n=\u2212log2n+\u222bt1t0,min\eta tdt=\u2212log2n\u2212\u222bt0,mint1\eta tdt.$ This gives $\u222bt0,intt0,max\eta tdt<0\u21d4\u2212log2n<\u222bt0,mint1\eta tdt$.

## Author notes

Laboratoire de Biochimie, École Polytechnique, CNRS, Université Paris Saclay Palaiseau 91128, France. E-mail: [email protected]