Evolution of cell size control is canalized towards adders or sizers by cell cycle structure and selective pressures

  1. Felix Proulx-Giraldeau
  2. Jan M Skotheim
  3. Paul François  Is a corresponding author
  1. Department of Physics, McGill University, Canada
  2. Department of Biology, Stanford University, United States
  3. Chan Zuckerberg Biohub, San Francisco, United States

Abstract

Cell size is controlled to be within a specific range to support physiological function. To control their size, cells use diverse mechanisms ranging from ‘sizers’, in which differences in cell size are compensated for in a single cell division cycle, to ‘adders’, in which a constant amount of cell growth occurs in each cell cycle. This diversity raises the question why a particular cell would implement one rather than another mechanism? To address this question, we performed a series of simulations evolving cell size control networks. The size control mechanism that evolved was influenced by both cell cycle structure and specific selection pressures. Moreover, evolved networks recapitulated known size control properties of naturally occurring networks. If the mechanism is based on a G1 size control and an S/G2/M timer, as found for budding yeast and some human cells, adders likely evolve. But, if the G1 phase is significantly longer than the S/G2/M phase, as is often the case in mammalian cells in vivo, sizers become more likely. Sizers also evolve when the cell cycle structure is inverted so that G1 is a timer, while S/G2/M performs size control, as is the case for the fission yeast S. pombe. For some size control networks, cell size consistently decreases in each cycle until a burst of cell cycle inhibitor drives an extended G1 phase much like the cell division cycle of the green algae Chlamydomonas. That these size control networks evolved such self-organized criticality shows how the evolution of complex systems can drive the emergence of critical processes.

Editor's evaluation

This paper develops evolutionary simulations to identify the type of molecular networks that can give rise to size control. The authors propose an evolutionary framework to find which factors select for particular mechanisms in cell size control. They show that the evolution of a specific cell size control mechanism is dependent on the cell cycle structure.

https://doi.org/10.7554/eLife.79919.sa0

Introduction

Cell size is fundamental to cell physiology and function because it sets the scale of subcellular compartments, cellular biosynthetic capacity, metabolism, mechanical properties, surface-to-volume ratios, and molecular transport (Chan and Marshall, 2010; Ginzberg et al., 2015; Neurohr et al., 2019; Zatulovskiy and Skotheim, 2020). While different types of cells vary enormously in size to perform their functions, cells within a particular type are generally uniform in size indicating that cell growth may be accurately coupled to division and differentiation processes. On a phenomenological level, there are many commonalities in how cells regulate their size even though the molecules controlling cell division vary across the tree of life with the most striking differences separating eukaryotes and bacteria. This extreme molecular diversity in the regulatory proteins raises the question as to what are the common features of the control systems that evolved to implement size control.

Most generally, cell size control can be viewed as a return map where the division size is a function of the cell size at birth. The examination of proliferating cells in laboratory conditions has revealed a variety of size control phenomena that can be characterized quantitatively by plotting the size of a cell at birth against the amount of mass added before it divides (Amir, 2014; Facchetti et al., 2017; Jun and Taheri-Araghi, 2015). A ‘sizer’ has a slope of –1 so that all variation in cell mass at birth is compensated for in one cell cycle, whereas an ‘adder’ has a slope of 0 so that each cell adds the same amount of mass during the cell cycle regardless of initial size. In the case of an adder, control is weaker so that multiple cell cycles are required for a particularly large or small cell to return to the average cell size. Importantly, the slope relating size at birth with the amount of growth in the cell cycle is a metric that quantifies the amount of size control occurring in a particular condition.

Studies of cell size control have revealed a diverse set of phenomena. Fission yeast and mouse epidermal stem cells exhibit ‘sizers’, and most bacteria, archaea, and cultured human cell lines exhibit behavior closer to adders (Cadart et al., 2018; Eun et al., 2018; Jun et al., 2018; Sveiczer et al., 1996; Westfall and Levin, 2017; Willis and Huang, 2017; Xie and Skotheim, 2020). Thus, while diverse size control behaviors have been observed, adders have been observed more often than sizers. This raises the question of why adders are more frequently observed if sizers, by definition, are more effective at controlling cell size (Barber et al., 2017; Willis et al., 2020).

To address the question of why adders are the most often observed form of cell size control, we used evolutionary algorithms (Holland, 1992) to identify commonalities between networks evolved to control cell size. Evolutionary algorithms are a class of machine learning techniques aiming at mimicking evolutionary processes (Crombach, 2021; François, 2014; François and Hakim, 2004; Xiong et al., 2019). Because of the nature of evolution, results of evolutionary computations are often more efficient and more creative than expected (Lehman et al., 2020). Furthermore, solutions found by evolutionary algorithms are constrained by their evolutionary paths followed and present similar characteristics to biologically evolved systems (Schaerli et al., 2018). Cell size is regulated through the cell cycle control network that governs transitions from one phase of the cell cycle to the next. The division cycle can be broken up into distinct phases that are characterized by different molecular activities (Morgan, 2007). While it is typically considered that there are 4 phases of the cell cycle (G1, S, G2, and M), we here consider a two phase model based on a G1 phase and a composite S/G2/M phase. This is because size control in general has been associated with either the G1/S transition or mitosis at the end of the cell cycle.

We start with a simple model of the cell cycle with a timer for an S/G2/M phase (Figure 1) that we evolve to optimize homeostatic cell size control (Figure 1B). We discovered that different control mechanisms could perform cell size control based on protein quantity or concentration. Simulations in which size control takes place in G1 phase converge toward an adder mechanism for the entire cell cycle and identified an active quantity sensing mechanism similar to dilution-based mechanisms previously identified experimentally (Chen et al., 2020; Qu et al., 2019; Schmoller et al., 2015). The relative durations of G1 and S/G2/M were important in determining size control properties. A relatively shorter S/G2/M phase favors sizer mechanisms, while longer S/G2/M phases favor adders. Moreover, inverting the model so that cell size controls S/G2/M and G1 is a timer, like in fission yeast, results in more sizer-like control. Thus, we anticipate adders arise when cell size is controlled at a point intermediate in the cell cycle, like the G1/S transition, while sizers will appear when cell size regulates a point later in the cell cycle, as is the case when G1 is proportionally longer, or control takes place at the transition to mitosis. We finally identify a self-organized mechanism based on fluctuation sensing where size control occurs on average over multiple cycles. While there is no one-to-one correspondence between a specific size control mechanism and a given evolutionary pressure, our work identifies clear evolutionary principles that shed light on the diverse cell size control phenomena previously observed experimentally.

Implementation of the cell cycle seed network and evolution algorithm.

(A) Schematic representation of the coupling between cell size and cell cycle progression. The transition between G1 (red) and S phases of the cell cycle at the G1/S transition (orange) can evolve to depend on cell size, while the duration of the S/G2/M (blue) phase is independent of size. Cell division takes place instantaneously following mitosis (purple). (B) Schematic representation of the φ-evo algorithm implementation. The network generating tool takes an initial network topology as its starting point for evolution as well as a user-defined fitness function. φ-evo then goes through successive epochs of mutation and selection to extract a final optimized network. At each selection step, the fittest half of the networks are retained and duplicated for evolution in the subsequent epoch. Interactions permitted to be mutated by φ-evo include transcriptional activation (green arrow), transcriptional repression (red arrow) and protein-protein interactions responsible for complex formation (black arrow). (C) Schematic of our seed network topology implementing a simplified relaxation oscillator. Cell cycle state (G1 or S/G2/M) is encoded via a binary switch called ‘S/G2/M Switch’ that is 0 in G1 and 1 in S/G2/M. Transition between G1 and S is controlled by the quantity of an inhibitor of the G1/S transition that we call I. The lower the quantity I, the higher the chance of progression through the G1/S transition. This interaction is represented as the grey arrow in the network topology and cannot be mutated by φ-evo. After progressing through G1/S, cells enter S/G2/M which we model as a pure timer of fixed duration with some uniform noise. Cell volume grows exponentially and is divided symmetrically following mitosis. We then follow one of the daughter cells and disregard the other one. (D) Phase-space representation of the initial relaxation oscillator. The X-coordinate shows the S/G2/M Switch variable, and the Y-coordinate shows the concentration of [I]. The oscillator runs counterclockwise with the left branch (x=0) corresponding to G1 and the right branch (x=1) corresponding to S/G2/M. The G1/S transition and division events are instantaneous in our simulations but are smoothly represented here for visualization purposes. From these transitions, we can extract the approximate shape of the [I] nullcline that we plot under the oscillator with a dashed line. We note that the position of the G1/S transition in phase-space will vary as a function of the volume of the cell as it depends on the quantity of I rather than its concentration [I].

Results

Initial cell-cycle model

In general, there are two classes of mechanisms that cells use to control their size that can be separated in terms of whether cell division or cell growth per se is regulated by cell size. Note that in this work, we will use mass, size, and volume interchangeably. In the first class, it is crucial that the growth rate per unit mass of a cell depends on cell size so that cells that are significantly larger than the optimum cell size grow slower (Cadart et al., 2019; Ginzberg et al., 2018; Miettinen and Björklund, 2016; Nordholt et al., 2020; Tzur et al., 2009). Such slower growing cells are then outcompeted by cells closer to the optimum size even when divisions occur purely by chance (Conlon and Raff, 2003). While size-dependent growth mechanisms exist and do support size homeostasis, such mechanisms rely on inefficient growth in all the cells away from the optimum size (Ginzberg et al., 2018; Miettinen and Björklund, 2016; Nordholt et al., 2020). To avoid such inefficient growth, many types of cells use active size control mechanisms to accelerate progression through the cell cycle in larger cells (Zatulovskiy and Skotheim, 2020). In our simulations, we keep cell growth rates constant over a physiological range of cell sizes. This allows us to focus on the common features of the molecular networks in which increasing cell size drives changes in molecular activities to trigger cell division. We assume that cell volume V grows at a rate λ(V)×V, so that growth is exponential when λ is a constant. Volume is divided by 2 at each division after which we follow one of the two daughter cells. The growth rate sets the time scale for the system dynamics as it defines the doubling time τ=ln(2)/λ. Any interdivision time shorter than τ will see the cell volume shrink at the next generation while any interdivision time larger than τ will see the cell volume grow. We also use the chemistry square bracket convention such that any protein X’s concentration is denoted by [X] . Correspondingly, its quantity is denoted by X only and is defined as X=[X]×V.

We initialize our network evolution simulations with a very simplified model of the cell-cycle (Figure 1A). We model two independent phases of a symmetrically dividing cell, G1 and S/G2/M, separated by a commitment point at the end of G1 and division at the end of S/G2/M (Figure 1A, Chandler-Brown et al., 2017). We encode this cell cycle state information via a binary switch variable we call ‘S/G2/M Switch’ that is 0 in G1 and 1 in S/G2/M. In all simulations, we follow an inhibitor model (Heldt et al., 2018; Schmoller et al., 2015; Zatulovskiy et al., 2020) and assume that the probability of passing the G1/S transition is controlled by the quantity of a transcription regulator I. One way the quantity rather than the concentration of a molecule could be sensed is through its titration against a fixed cellular quantity such as the genome, which is part of a general class of titration-based cell size sensing mechanisms (Amodeo et al., 2015; Heldt et al., 2018; Si et al., 2019; Wang et al., 2009). A lower quantity of this inhibitor I means a higher probability of a G1/S transition at the current time step of the simulation (Appendix 1—figure 2). Like all other proteins, the quantity I is produced with a rate proportional to volume, degraded at a constant rate, diluted by cell growth, and equally partitioned between mother and daughter cells at division (see Materials and Methods). We found that due to the volume scaling assumption, [I]’s concentration alone was largely independent of volume and could not trigger a size-dependent G1/S transition, which is why we opted for the quantity of I instead (Appendix 1—figure 3). Upon passing the G1/S transition, we assume cells are committed to division and there is a fixed time delay before they divide thus modeling S/G2/M as a timer (Chandler-Brown et al., 2017; Doncic et al., 2015). We initially fixed the timer duration to be roughly equal to 50% of the doubling time τ with some uniform noise such that G1 and S/G2/M durations would be the same at equilibrium. Regulation of the quantity I during the cell cycle thus controls the precise timing of the G1/S transition, but it is not always perfect since the transition is probabilistic. This, along with the noise in S/G2/M timer duration, creates natural cell to cell variability in volume that needs to be compensated for by the evolved mechanism. We note that we initialized most of our simulations with one added interaction in which production of [I] is activated by the S/G2/M Switch variable to reset its concentration to a higher level before the next generation. We initially ran simulations without this specific interaction but found that it systematically appears in the initial stages of evolution simulations. We therefore included it in the initial network to speed up our simulations. We refer the reader to the Appendix 1 for more details. The models used in this study are publicly available (Proulx-Giraldeau and François, 2022).

Typical dynamics of this simple cell-cycle model are represented in Figure 1. These dynamics are similar to models of cell cycles based on relaxation oscillators (Cross, 2003; Tsai et al., 2008). The left and right slow branches correspond to G1 and S/G2/M, respectively, and the fast horizontal branches represent G1/S and division. An intermediate fictitious nullcline is shown as a line that connects the average concentration [I] at G1/S and at division. Starting with cell birth, the system goes down the left-G1 branch because of degradation, then jumps to the right-S/G2/M branch below the threshold for the G1/S transition, stays there while moving up due to production by the S/G2/M Switch, until it jumps back to the left branch at the end of the timer phase. We note that there is no explicit volume control in this initial model since the only control comes from the quantity of I which does not initially depend in any way on the volume. This initial quantity sensing oscillator does not perform size control and instead results in unstable growth where size deviations are amplified at each generation instead of being corrected (Appendix 1—figure 3) as had been previously described for a size scaling inhibitor dilution model (Barber et al., 2017; Willis et al., 2020). Thus, the network needs to evolve some other interactions and/or parameters to go beyond a simple G1 inhibitor driven by production in S/G2/M to create a viable cell lineage.

Evolution of quantity-based size control mechanisms

To examine how networks controlling cell size could evolve, we ran evolutionary simulations that optimize both the number of divisions NDiv and the coefficient of variation of the size distribution at birth CVBirth (see Materials and Methods for algorithm details). Figure 2A illustrates the behavior and results of a typical evolutionary run, with axes defined by both fitnesses used. Simulations successfully evolving size control mechanisms typically follow the same pattern. Networks initially cluster in two regions: region [ii] where cells have low CVBirth but grow too small and die after a few divisions, and region [i] where cells grow too big and reach our cut-off for fitness 2 (y-axis). Notice that our Pareto evolutionary algorithm maximizes network diversity, so that those two clusters are at first maintained during evolution (rank 1 Pareto networks Warmflash et al., 2012). As the number of epochs increases, networks in cluster [ii] have more and more divisions, but still grow too big, so that those cells are therefore penalized (see details in Appendix 1). At some point in evolution (around epoch 700 for this particular simulation), some weak control mechanism suddenly evolves, preventing cells from becoming too big without imposing a tight control on the average volume (see also Figure 5 for an explicit example of how this is done). Thus, fitness 2 collapses and the number of divisions is optimized simultaneously. Cells later optimize the control to give a lower CVBirth . The optimal networks, at the right most end of this line, both maximize NDiv and minimize CVBirth .

Evolution of feedback-based size control.

(A) Typical 2D fitness trajectory for an evolutionary run. Individual networks are dots color coded by their epoch within the evolutionary trajectory. Fitness function of the number of divisions of a cell lineage during a time interval of fixed length (F1(NDiv), X-coordinate) and fitness function of the coefficient of variation of the volume distribution at birth (F2(CVBirth), Y-coordinate). Optimal model behavior is located in the bottom right corner of the figure where networks produce cell lineages with many offspring and strong size control. First, there are several epochs without any size control; networks cluster in two regions of the Pareto front corresponding to volume going to the maximum allowed value (cluster [i]) or to the minimum value (cluster [ii]). Both cases are highly penalized in their fitness score. Evolution goes back and forth between the [i] and [ii] clusters with a slow increase in the number of divisions (X-coordinate). Eventually, some volume control evolves and networks transition in the [iii] cluster where their CVBirth is slowly optimized further until the end of the run. (B) Core network topology of the evolved Model A1 network that employs a feedback-based mechanism described in detail in panels C-F. (C) Concentration of [I] at the G1/S transition (Y-coordinate, left axis, orange) and average concentration of the repressor protein [R] in S/G2/M (Y-coordinate, right axis, light blue) as a function of the volume of the cell at G1/S (X-coordinate), that is, the beginning of S phase. We see here that [R] acts as a direct size sensor of the volume at G1/S. (D) Quantity of inhibitor I at birth as a function of volume of the cell at birth, which is independent of size due to titration by [R] during S/G2/M. (E) Trajectories of quantity of inhibitor I in G1 as a function of time. Trajectories are color coded as a function of the volume at G1/S during the previous generation’s cell cycle. For visualization purposes, trajectories are offset vertically to all begin at the average quantity of I at birth (t=0) shown to be on average independent of volume in panel D. Larger cell volumes lead to greater titration of [I] in G1 by [R]. In turn, this ensures that G1 duration of the daughter cell cycle is shorter, which underpins the size control mechanism. (F) Characteristic dynamics of Model A1. Circles indicate volume at G1/S. Extrinsic perturbations are applied to the model by temporarily changing the division ratios which kicks the system out of equilibrium at the subsequent cycle such that VNextG1/S=(1±0.5)VG1/S. The volume relaxation back to its homeostatic value takes ~2–3 generations, almost insensitive to the fact that the perturbation is applied towards higher or lower volumes. (G) Amount of volume added ΔV in G1 (red), S/G2/M (dark blue), and over the whole cycle (purple) as a function of their initial volume at the beginning of these phases, that is, birth for G1 and cycle, and G1/S for S/G2/M, with the slope of linear fits indicated in legend. Slope of 1 corresponds to a Timer, slope of –1 to a Sizer and slope of 0 to an Adder. (H) Network topology of Model A2, a second evolved network that is similar to Model A1 albeit with different kinetic parameters and 2 additional interactions (see text). (I) Characteristic dynamics of Model A2. Extrinsic perturbations are applied like in (F). The volume relaxation back to its homeostatic value takes ~3–4 generations when applied towards the higher volumes but only 1 generation when applied towards the lower volumes. (J) Amount of volume added for different periods of the cell cycle for Model A2.

Evolution simulations are in part reproducible and most often lead to similar network topologies. The evolution trajectory leading to Model A1 is a typical example (Figure 2B, see also variations of this network in Model A2 in Figure 2H and models A3-6 in Appendix 1—figures 811). The minimal network common to all those models is very simple. One gene, R, is added to the seed network and is both repressed and titrated by I forming the network motif known as a Mixed Feedback Loop (François and Hakim, 2005). Size control can then be understood intuitively as follows.[I] represses [R], which is thus only produced in the narrow window of the cycle when [I] is low, i.e., when the cell is close to the G1/S transition and in early S/G2/M. But, since the quantity I=[I]×V is fixed at G1/S by design, the concentration [I] is inversely proportional to the volume of the cell at the G1/S transition (VG1/S) as shown in our simulations (Figure 2C). Because of this, the [I] dependent synthesis rate of [R] and therefore its subsequent concentration are (linear) functions of the volume of VG1/S (Figure 2C), allowing for the cell to keep a memory of its volume at G1/S via the [R] variable (this holds even once [R] is constantly degraded for the remainder of the cycle). This has two effects. First, during S/G2/M, I synthesis rate is proportional to volume by hypothesis (and thus to VG1/S), and I is titrated by [R], also proportional to VG1/S. Both effects even out so that cells are born with a fixed quantity of inhibitor I that is independent of volume (Figure 2D). Second, after division, production of I is 0 by hypothesis, but I still is titrated in G1 by the remaining [R] (still proportional to VG1/S). Because I quantity at the beginning of G1 is size independent, this ensures that daughter cells reach the I quantity threshold of G1/S earlier if they were born larger, thus ensuring size control. To confirm this, we examine the change in quantity of I as a function of time spent in G1 and of VG1/S, and we see the slope of these two quantities is volume-dependent (Figure 2E). We notice that this evolved size mechanism likely is the simplest possible allowed by our formalism: on the one hand, it entirely captures the volume dependency in a single variable [R], and on the other hand, it ensures proper scaling of I both at birth and at G1/S for size control with the help of a single titration. Notice that such simple control also explains the sudden evolutionary “jump” of Pareto front around epoch 700 on Figure 2A, which corresponds to when the Mixed Feedback Loop motif first appears.

These evolved size control networks, while relying on quantity sensing, are conceptually similar to the budding yeast network relying on concentration sensing of the cell cycle inhibitor Whi5 since there is a constant quantity of I present right after division (just like Whi5). In budding yeast, the Whi5 protein is passively diluted in G1 to increase the stochastic rate of progression through the G1/S transition. The time spent in G1 depends on the initial concentration of Whi5 at birth, which scales as 1V , to promote a sizer mechanism. Here, the concentration of [I] at birth scales as 1V , but so does the threshold concentration of [I] regulating the G1/S transition. This is precisely why an active titration mechanism is required to obtain G1 size control in our setup. Such homeostatic control ensures cell size returns to its steady state distribution following an artificial perturbation as soon as a volume deviation is detected at G1/S (Figure 2F). When we plot the amount of volume added ΔV for different phases of the cell cycle as a function of the initial volume at the beginning of these phases, we find an approximate adder over the whole cycle that results from weak sizer in G1 followed by a timer in S/G2/M as has been found in budding yeast (Chandler-Brown et al., 2017; Di Talia et al., 2007; Soifer et al., 2016, Figure 2G).

While we chose one simple model to illustrate the control mechanism common to our set of evolved networks (Model A1 shown in Figure 2B), other evolved networks were more elaborate but illustrated a similar principle. For example, Model A2 contains extra interactions for the volume sensing gene [R], where [R] is repressed by the S/G2/M Switch (meaning its production is completely shut down in S/G2/M leading to sawtooth-like dynamics). Furthermore, [R] represses the synthesis of [I], adding another layer of repression to promote size control beyond the previously described titration by [R] (Figure 2H). If we perturb cell size to examine the dynamics of the return to steady state and look at the added volume during the cell cycle, we again see overall a weak adder behavior similar to that found in Model A1 (Figure 2I–J). We give additional examples of similarly evolved networks in Appendix 1—figures 811 where we can see the sensing and the feedback mechanism being implemented in slightly different ways. Yet, despite these mechanistic differences in feedback regulation the resulting function of the evolved networks were similar as indicated by their CVBirth (Appendix 1—table 1).

Quantifying size control

To study the mechanisms implicated in cell size control, we modify the control at G1/S and introduce the control volume VC . This control variable is independent from the biochemical network and is maintained fixed allowing us to disconnect the actual cell volume V from the biochemical network and by forcing the G1/S transition to be triggered once IC=[I]×VC is low enough. We then numerically integrate the differential equations of the model and measure the period T(VC) of the simulated cell-cycle for this control volume. Use of the control volume allows us to break the size feedback system and distinguish its input, VC , from its output, the induced cycle period T (Angeli et al., 2004). We compute T(VC) for Models A1 and A2 and compare their responses with the analytical curves of the archetypical timer, adder, and sizer (Figure 3A and C). Those curves intersect at the point where the induced period is exactly equal to the population doubling time (τ=ln(2)/λ), which defines the equilibrium volume achieved by our cell size control network corresponding to VG1/S. Examination of the size control in different cell cycle phases indicates the contributions of G1 and S/G2/M to the overall system behavior (Figure 3B and D). We note that we later examine statistics of ensembles of evolved models but that Models A1 and A2 are both typical examples of evolved feedback-based models.

Characterizing and comparing evolved size control mechanisms.

(A) Average period T(VC) of the oscillator of Model A1 as a function of the control volume VC at the G1/S transition. Period is normalized by τ the doubling time of the cell, and volume is rescaled by VG1/S, which corresponds to T(VG1/S)=τ. Normalized periods larger than 1 indicate cell lineages that grow over time whereas normalized periods smaller than 1 indicate lineages that shrink over time. Periods for the sizer (red), adder (orange) and timer (dark blue) are shown for comparison. The S/G2/M timer period is incompressible and prevents a perfect sizer from existing in the large volume range as indicated by the red dotted line. Model A1 follows approximately the adder archetype over a large range of control volumes. (B) Added volumes ΔV for different phases of the cell cycles for simulations of Model A1. Individual dots correspond to different cell cycles for a simulation at steady-state. The full line corresponds to the extrapolation from the T(VC) curve shown in A for a restricted range of VC relevant to the scatter. The black cross, star and square indicate the average added volumes corresponding to when the system senses a volume corresponding to VG1/S at the G1/S transition. We see that the model is predicted to follow an adder over a large range of volumes. (C) Average period T(VC) of the oscillator of Model A2, with similar conventions as for panel A. We note that the T(VC) curve of this model is closer to the sizer at lower volumes and closer to a weak adder/timer at higher volumes relative to VG1/S. (D) Added volumes ΔV for different phases of the cell cycles for simulations of Model A2, with similar conventions as for panel B. Here we see the predicted sizer behavior at lower volumes and the weak adder/timer behavior at higher volumes.

Quantifying precisely how the cell cycle period depends on the control volume at G1/S allows us to see that the T(VC) for Model A1 overlaps with the theoretical adder curve over a broad range of volumes. In contrast, Model A2 behaves as a sizer for volumes smaller than the equilibrium volume. However, at higher volumes it behaves more like an adder/timer similar to Model A1. Model A2’s equilibrium is thus ‘tuned’ by evolution to be in a regime corresponding to the minimum of the ΔVCycle curve extrapolated from the T(VC) as shown in Figure 3D, precisely when the system transits from a sizer at small volumes to an adder-like behavior at higher volumes. Similar behavior was observed experimentally in budding yeast (Chandler-Brown et al., 2017; Delarue et al., 2017). Thus, both Models A1 and A2 approximate an adder near the equilibrium size, but their behavior differs further from equilibrium where for smaller volumes, Model A1 is still an adder while Model A2 is a sizer.

Modulating cell cycle structural constraints selects for adders or sizers

So far, our evolutionary algorithm selects networks that implement adders rather than sizers near the equilibrium size. This is surprising because sizers are in principle better than adders at controlling cell size and reducing the CV of the size distribution at birth, which is one of our fitness functions. That adders are more frequently observed in nature than sizers (Cadart et al., 2018; Eun et al., 2018; Jun et al., 2018; Westfall and Levin, 2017; Willis and Huang, 2017; Zatulovskiy and Skotheim, 2020), is consistent with our evolution simulations, but adds to the mystery as to why this takes place.

To gain insight into the underlying reason for the prevalence of adders, we considered what might be exceptional in the cases where sizers occur. The best studied, and highly accurate sizer, is found in the fission yeast S. pombe (Fantes, 1977; Sveiczer et al., 1996). In contrast to budding yeast and human cells, where the size control takes place largely in G1 phase, fission yeast exerts size control later in the cell cycle at the G2/M transition. That size control took place later in the cell cycle in fission yeast suggested that this structural feature of its cell cycle might be responsible for its stronger size control. To test this, we inverted our seed network so that the G1 phase was a timer, while cell size control could evolve in S/G2/M (see Materials and Methods). We compare these results to the evolutionary simulations starting with the ‘control’ seed network with G1 size control and an S/G2/M timer (Figure 2). We note that S. pombe growth rates have been reported to deviate from exponential (Wood and Nurse, 2015). While slower than exponential growth would aid cell size control, our analysis here is restricted to exponential growth.

To determine how the seed network structure influences the subsequent evolution of cell size control, we performed 120 independent evolutionary simulations for the two network structures initialized with the Model A1 topology and parameters (Figure 4). Sixty simulations were performed using Pareto optimization and another 60 simulations were performed using individual fitness optimization based on the number of cell divisions NDiv . For each simulation’s most fit model after 500 epochs, we calculate the CV of the volume distribution at birth, CVBirth , and the slope of the linear fit of the volume added in the entire cell cycle as a function of the cell size at birth, Slope ΔVCycle (we remind the reader that a slope of –1 corresponds to a sizer, 0 to an adder, and 1 to a timer). We chose to use 500 epochs in our simulations because in our previous experience this was sufficient for networks to evolve to be near the optimum, but not so much that they were forced to extensively explore the effects of neutral mutations near the optimum. Model A1’s initial and Slope before evolution are indicated by the dashed black line. In the control experiment where G1 is a sizer and S/G2/M is a timer, most evolutionary simulations with two fitness functions (Pareto) yield models close to the adder regime with a low CVBirth . However, when only the number of divisions (NDiv) is used as a fitness, the evolutionary simulations are closer to the sizer regime, albeit with a slightly higher CVBirth (Figure 4A, Welch’s t-test p<10–4). When the cell cycle structure is inverted so that G1 is a timer and S/G2/M is a sizer like it is in the fission yeast S. pombe, we found that more sizer-like networks evolve than in the control experiment (Figure 4B, Welch’s t-test on agglomerated data p=0.03). This shows that having a network structure like the fission yeast S. pombe promotes sizers, while having the size control portion of the cell cycle earlier, as in the budding yeast S. cerevisiae, promotes adders. Thus, performing simulations using cell cycle network structures of these two yeasts results in the evolution of size control mechanisms that reflect those that are naturally occurring.

Distinct network constraints and selection pressures bias size control evolution towards adders or sizers.

Summary statistics for evolutionary simulations each having 500 epochs. Model A1 shown in Figure 2A-G was used as the initial seed network. 60 simulations were performed using Pareto optimization of the number of divisions (NDiv) and the CV of cell size at birth (CVBirth), are labeled Pareto and are shown in full colors. 60 more simulations were performed using only the number of divisions as the fitness function, are labeled NDiv and are shown in colored outlines only. Scatter plots show the coefficient of variation of the size distribution at birth (CVBirth, Y-coordinate) as a function of the fitted added volume slope over the whole cycle as a function of volume at birth (Slope ΔVCycle , X-coordinate) for the most fit models evolved during each of the 120 independent simulations. Horizontal box plots above the scatter plots display the distributions of the added volume slopes for the Pareto and NDiv simulations. Timer (dark blue), adder (orange) and sizer (red) slopes are shown respectively at 1, 0, and –1 for comparison. Vertical box plots on the right of the scatter plots show the distributions of CVBirth for the Pareto and NDiv simulations. Asterisks represent p-values for the Welch’s t-Test between the distributions. For reference, NS indicates p>0.05, * indicates p<0.05, ** indicates p<102 , *** indicates p<103 and **** indicates p<104 . The values of CVBirth and Slope ΔVCycle for the initial seed Model A1 are shown as a black square in the scatter plot or as a dashed black line in the box plots. Each panel explores different cell cycle structures which are summarized by the pie charts. Cycles begin on the left of the pie charts and rotate clockwise, indicating the order of the sizer (red) and timer (dark blue) phases. The labels indicate each phase’s duration at equilibrium as a percentage of the doubling time τ. (A) Identical evolutionary parameters as for Model A1 evolution shown in Figure 2A–G. G1 performs size control and has a duration 0.46τ at equilibrium and S/G2/M is a timer of duration 0.54τ. (B) Evolution results for a cell cycle structure where the sizer and the timer phases of the cell cycle are inverted akin to S. pombe. G1 is a timer of duration 0.54τ and S/G2/M performs size control and has duration 0.46τ at equilibrium. (C) Evolution results for a G1 size control of average duration 0.64τ at equilibrium where S/G2/M is a timer of duration 0.36τ. (D) Evolution results for a G1 size control of average duration 0.28τ at equilibrium where S/G2/M is a timer of duration 0.72τ.

The general notion that having cell size control in G1 results in adder-like mechanisms, while control later in the cell cycle results in more sizer-like mechanisms fits most observations of human cell lines grown in culture, budding yeast, and fission yeast. However, a recent study examining mouse epidermal stem cells growing and dividing in the skin found both a strong sizer and that this sizer was largely due to the size-dependent regulation of G1 (Mesa et al., 2018; Xie and Skotheim, 2020). This raised the question as to how a network performing size control in G1 could result in a sizer for the entire cell cycle. One important difference between mammalian cells grown in culture and the mouse epidermal stem cells growing in an animal (in vivo) is the change in the relative durations of the G1 and S/G2/M phases of the cell cycle. While the S/G2/M phase of the cell cycle is similar in duration in cultured cells and the epidermal stem cells in vivo at ~12 hr, the G1 phase extends ~fivefold from ~10 hr in culture to ~50 hr in vivo (Cadart et al., 2018; Xie and Skotheim, 2020). This suggests the hypothesis that the overall size control behavior can be dominated by the relatively longer cell cycle phase, as is likely the case for the G1 phase of epidermal stem cells. To test this hypothesis, we performed evolutionary simulations with size control in G1 and a timer in S/G2/M but where we changed the duration of the S/G2/M timer phases of the cell cycle to be significantly shorter or longer than the G1 phase at equilibrium. When a timer in S/G2/M is relatively shorter compared to G1, we generally see more sizer-like behavior can evolve (Figure 4C, Welch’s t-test on agglomerated data p=4 x 10–3), while when it is relatively longer, we see more adder-like or even timer-like behavior (Figure 4D, Welch’s t-test on agglomerated data p<10–4).

We next considered the effect of changing the amount of noise in the timer phase of the cell cycle. To do this, we examined the evolution of networks performing size control in G1 and where the S/G2/M phase with an increasing amount of noise. Increasing the noise in the timer progressively reduced the amount of size control done by the network (Appendix 1—figure 5). This is likely because the fixed duration of S/G2/M allows the system to accurately reset protein concentrations for the subsequent cell cycle to promote accurate G1 control (Willis et al., 2020). We also examined the effects of adding noise to the cellular growth rate and to volume partitioning at division and found similar results (Appendix 1—figures 67).

Taken together, our simulations show how the structural features of the cell cycle are important for determining what type of size control ultimately evolves. G1 control is more conducive to the evolution of adders, while S/G2/M control is more conducive to sizers. Moreover, size control can be dominated by the cell cycle phase of longest duration and is modulated by the specific selection criteria.

A two-step evolutionary pathway for cell size control

From our series of evolution simulations, we found a somewhat paradoxical inverse correlation between the added volume slope quantifying the degree of size control (sizer vs. adder) and the CVBirth (Figure 4B–C). This was surprising because it means that there is a broader distribution of volumes in the sizer regime where control should be more effective in theory. To better understand this inverted correlation between Slope ΔVCycle and CVBirth , we revisited our evolutionary simulations to examine the evolutionary pathways through which the networks progressed through the simulated epochs.

A typical evolutionary pathway for a Pareto simulation of the S. pombe-like network structure presented in Figure 4B is shown in detail in Figure 5. Here, Model A1 topology is conserved throughout evolution although individual parameter values change. In the early stages of the evolution (epoch 650), we typically see dynamics where small cell size triggers an overshoot to a large cell size, which is then reduced through a series of rapid divisions (Figure 5A). Because of this small-size-triggered overshoot, the system behaves more like a sizer when the average behavior is analyzed (Figure 5B). However, the high degree of variability in the cell cycles also leads to a broad distribution of volumes (Figure 5C). The variability in volume is attenuated in later epochs where there are fewer and smaller volume overshoots triggered by small cell size (Figure 5D–I). Since small cell size no longer triggers a dramatic amount of cell growth, the slope of ΔVCycle is increased and the system converges towards a weak adder in which the distribution of volume at birth is more Gaussian and the CVBirth is lower (Figure 5H–I). Taken together, these analyses suggest a two-step evolutionary pathway, consistent with the evolutionary dynamics first seen in Figure 2A. First, a strong but imprecise sizer mechanism evolves where, because of noise in the system, small variations in volume lead to a dramatic overcorrection and overshoot of the target volume. The variability in volume produced by this overshoot is then reduced by attenuating the strength of the size control response. Indeed, the overall weaker size control allows the system to respond more mildly to size deviations, thus yielding a lower CVBirth overall which we select for. Thus, selecting for a smaller CVBirth , i.e. better size control, can end up selecting for adders rather than sizers. This paradoxical result is consistent with the fact that when we select only for the number of cell divisions (NDiv), one sees that more sizer-like behavior can evolve (Figure 4A–C). The typical behavior before optimization of CVBirth is a strong sizer as illustrated in Figure 5A.

System and evolutionary dynamics of cell size control networks.

Snapshots of an evolutionary simulation of 2500 epochs initialized with the Model A1 network topology along with an S. pombe-like cell cycle structure with a timer in G1 followed by a sizer in S/G2/M (see Figure 4B). Pareto fitness optimization was performed using NDiv and CVBirth as fitness functions. Rows indicate simulation results for the fittest networks from evolutionary epochs 650 (Panels A-C), 1000 (Panels D-F), and 2500 (Panels G-I). Network topology remains the same throughout the evolutionary simulation and is shown on the right. Evolutionary dynamics continually reduce the selected for CVBirth and proceed through a noisy sizer to a less noisy adder. (A) Typical dynamics of the most fit model from epoch 650. (B) Added volumes ΔV for different phases of the cell cycle for the most fit model from epoch 650. Fitted slopes are indicated in the legend. Fits for the S/G2/M and Cycle added volumes were split in two separate the size control for small and large cells. (C) Size distributions at birth (red), G1/S (orange), and division (purple) for the most fit model from epoch 650. The coefficient of variation of the volume distribution at birth CVBirth=0.279. (D) Typical dynamics of the most fit model from epoch 1000. (E) Added volumes ΔV for different phases of the cell cycle for the most fit model from epoch 1000. Fitted slopes are indicated in the legend. (F) Size distributions at birth (red), G1/S (orange), and division (purple) for the most fit model from epoch 1000. The coefficient of variation of the volume distribution at birth CVBirth=0.090. (G) Typical dynamics of the most fit model from epoch 2500. (H) Added volumes ΔV for different phases of the cell cycle for the most fit model from epoch 2500. Fitted slopes are indicated in the legend. (I) Size distributions at birth (red), G1/S (orange), and division (purple) for the most fit model from epoch 2500. We see here that the sizer behavior from epoch 1000 was abandoned for a weaker adder overall yielding lower CVBirth=0.063.

Fluctuation sensing and the evolution of self-organized criticality

One of the main features of smaller cells is that they have fewer proteins and mRNA. If some aspects of protein synthesis and degradation are subject to Poisson fluctuations, we expect such fluctuations to produce larger concentration fluctuations in smaller cells. For example, let us assume that the balance of synthesis and degradation of a generic protein results into a Poisson distribution with parameter ρVδ , where ρ is the synthesis rate in number of proteins per unit of time for a reference volume of 1, V is the volume, and δ is the degradation rate. The average concentration of this protein in an exponentially growing cell will be ρδ , which is independent of the cell volume V as expected from the production rate scaling. However, following the Bienaymé formula, the variance in the concentration is ρVδ×1V2=ρδV (Figure 6A), which decreases with volume. This result makes intuitive sense because bigger cells have to produce more proteins to keep concentrations constant, so that the fluctuations in the relative number of proteins (and thus concentration) are smaller (see Jia et al., 2021 for a complete analytical study of how in general variance scales differently from mean when volume varies). Thus, if the cell could sense the size of concentration fluctuations in some way, it would be able to harness the cell size-dependence of such stochastic fluctuations to regulate cell division and control cell size.

An evolved concentration fluctuation sensing size control model exhibits self-organized criticality.

(A) Size-dependent molecular noise arises due to Poissonian fluctuations in molecule number. Consider a protein production-degradation scheme for protein quantity X with production rate ρ and degradation rate δ contained in a volume V. At equilibrium, the concentration of [X] will be given by a distribution with mean ρδ and variance 1V . The effect of V on the molecular noise is shown on the time trajectories for simulations in a constant volume of V=1 (light blue), V=5 (medium blue), and V=50 (dark blue). Corresponding concentration distributions are shown on the right-hand side of the panel. (B) Network topology of Model B which evolved to sense fluctuations. Here, the G1/S transition is controlled by the concentration of [I] and not its quantity. (C) Characteristic cell cycle dynamics of Model B. Trajectories of [I], [A], and S/G2/M Switch are rescaled with arbitrary units (AU) for visualization purposes. Below, we zoom-in on three cycles to show how a low-volume induced burst in [A] leads to a massive production of [I] inducing a temporarily prolonged G1 phase. Subsequently, cells become bigger and display lower molecular noise inducing a comparatively shorter G1 phase. (D) Volume distributions at birth (red), G1/S (orange), and division (purple) for Model B. The coefficient of variation of the volume distribution at birth CVBirth=0.235. (E) Amount of volume added ΔV in G1 (red), S/G2/M (dark blue), and over the whole cycle (purple) as a function of their initial volume at the beginning of these phases, i.e., birth for G1 and cycle, and G1/S for S/G2/M, with the slope of linear fits indicated in legend. (F) Complementary cumulative distribution functions of the cycle duration (CCDF; probability that the cell cycle duration is larger than the value on the X-axis) for three models discussed in the main text: Model A1 (light blue), Model A2 (dashed dark blue), and Model B (red). The light grey line indicates the doubling time τ. We see that Model B exhibits a long tail past the doubling time, which is consistent with a power-law scaling of the cycle duration probability. We find a criticality indicative scaling exponent of –3.05 for the CCDF after fitting the tail of the distributions of 5 independent realizations of the dynamics of Model B. (G) Box plots of the cycle length distributions as a function of control volume VC at birth. Cycle lengths are normalized by the doubling time τ. Here, VC sets the molecular noise level to be equivalent to that of an exponentially growing cell born at VC but whose volume is reset to VC at each division. Note the very long tail of the distributions at small VC . (H) Position of the G1 attractor for inhibitor [I] as a function of activator [A]. The black dashed line corresponds to the level of [A] which triggers a transition between two modes of growth and division as shown by the position of the G1 attractor for [I] becoming equal to the concentration required to induce the G1/S transition. The two modes of growth are labeled [i] and [ii] and are also indicated in panels I-K. (I) Activator protein [A] is synthesized in bursts whose amplitude and duration are a function of volume. We define the burst duration as the total time during which [A] >0 for a cycle. The burst amplitude corresponds to the average level of [A] during each G1 phase. Each burst is then color-coded as a function of the birth volume of the cell that induced it. We use a divergent colormap whose center value (light yellow) corresponds to the average volume of the cells at birth and is indicated by a notch on the colorbar. Here, [i] corresponds to the deterministic regime when volume is high and [ii] corresponds to the noisy regime when volume is low. Note that the average volume of the cells at birth is positioned close to the black dashed line. (J) Phase-space representation of the relaxation oscillator in the deterministic regime [i]. The X-coordinate shows the S/G2/M Switch variable, and the Y-coordinate shows the concentration of [I]. Here, when volume is high, the position of the G1 attractor is below the [I] concentration at which the G1/S transition happens. Thus, G1/S takes place and cells are in the cell cycle with a period of ~0.85. (K) Phase-space representation of the noisy regime [ii]. Here, when volume is low, the position of the G1 attractor becomes greater than the [I] concentration at which the G1/S transition happens, and cells remain temporarily stuck in a prolonged G1 state and are unable to trigger the G1/S transition.

To test if we can evolve networks that control cell size through sensing Poissonian fluctuations, we first initialized our simulation with a network similar to that shown in Figure 1C, but with an added self-activating gene A that can activate the production of the I inhibitor. We then ran evolutionary simulations using the cell cycle structure of a sizer controlling G1 and timer in S/G2/M, but where the G1/S transition is regulated by the concentration [I] instead of its quantity. We also used Pareto fitness optimization of NDiv and of CVBirth . Importantly, we use the stochastic version of our equations with the molecular noise modeled using a Langevin noise term with a variance inversely proportional to the volume as explained above. We then extracted the most fit network and optimized it further using Pareto optimization of NDiv and of the ΔVCycle slope to push the models towards the sizer regime. The most fit network of those combined evolutionary simulations is presented in Figure 6B and demonstrates that we can indeed evolve fluctuation-based cell size control (Model B).

The mechanism for size control that evolved based on size-dependent fluctuation sensing is remarkably similar to what we observed for models without size-dependent fluctuations (Figure 5A). For large volumes, the cycle has a constant period which corresponds to approximately 85% of the doubling time τ. This ensures that in the high-volume regime, the system shrinks over time. When the volume is small however, fluctuations allow the concentration [A] to cross the threshold of the highly non-linear transcriptional activation of [I] by [A]. This results in a massive increase of [I] that needs to be degraded to progress further into the cell cycle. Thus, the low volume regime occasionally leads to a considerable increase in G1 length and a correspondingly very large cell at division. These very large cells then reliably and deterministically re-enter multiple, rapid cell cycles with short G1 until the cell is small again and the concentration fluctuations again become large enough to trigger the activation of [I] by [A]. This mechanism thus appears very similar to the early sizer mechanism observed in other quantity sensing simulations shown in Figure 5. However, here the mechanism is based only on size-dependent fluctuations in protein concentration and the overall behavior is closer to an adder (Figure 6E).

The system dynamics that evolved to perform fluctuation-based cell size control produce volume distributions that are long-tailed due to the stochastic occurrence of occasional exceptionally long G1 phases. Interestingly, the probability distribution of cell cycle durations follows a power law (Figure 6F), which is due to the very broad distributions of G1 duration at lower cell volumes. A more controlled analysis specifying the initial conditions showed that cell cycles get increasingly long, and their distributions widen with decreasing control volume VC (Figure 6G). Such non-Gaussianity is the hallmark of critical behavior, suggesting that the evolution of fluctuation-based cell size control is based on self-organized criticality (SOC). SOC is defined as a system where an order parameter feeds back on a control parameter (Sornette et al., 1995; Vidiella et al., 2021). The canonical example of SOC is the sandpile to which grains of sand are added on top. As the sand accumulates, the slope steepens, and the angle of the pile (control parameter) increases. Eventually, this triggers avalanches (order parameter) that feedback to dramatically reduce the angle of the pile. This ensures that the system dynamically tunes itself at the critical value of the angle of the pile where avalanches can occur.

We conclude that our evolved size control network exhibits SOC based on several observations. Starting from a high volume, multiple divisions at a rate faster than it takes to double the biomass reduce cell volume V just like the addition of grains of sand gradually increases the slope of the pile. Then, for small enough volumes, bursts of [A] drive an extended G1 that greatly increases cell size, which, like the sandpile avalanches, resets the system’s control parameter (volume of the cell or angle of the sandpile). Interestingly, evolution tuned the system to be near a bifurcation (Figure 6I). If we consider the deterministic regime, in which the fluctuations in [A] are small, the cycle is unperturbed and oscillates with a period roughly equal to 85% of the doubling time (Figure 6J). In contrast, if we consider the noisy regime, in which the fluctuations in [A] are large, the cycle disappears, and the system stays locked in a prolonged G1 state with a high value of [I] which is akin to a bifurcation destroying the cycle (Figure 6K). This bifurcation takes place because the position of the G1 attractor for [I] becomes larger with increasing [A] and eventually overcomes the concentration required to induce the stochastic G1/S transition (Figure 6H). Then, the system remains stuck in a state where G1/S cannot be triggered, and cells effectively exit the cycle. As growth occurs, noise dies down and so does the position of the G1 attractor, eventually becoming smaller than the [I] concentration required to induce the G1/S transition which allows cells to re-enter the cell cycle. Thus, the system is critical from a dynamical systems standpoint and also fits the general observation that SOC systems tune themselves to be right at the point where the order parameter is non zero, but infinitesimal (Sornette et al., 1995). In our case, the bifurcation corresponds exactly to the point where [A] can sufficiently activate the production of [I] to prevent the G1/S transition.

Discussion

The last decade saw an explosion of time lapse microscopy studies measuring how cells control their size. These studies revealed diverse phenomena that are characterized by the correlation between cell size at birth and cell size at division. Size control ranged from sizers, where the size at division is uncorrelated from the size at birth, to adders, which add a constant volume in each cell division cycle, to timers, whose cell cycle duration is size-independent (Cadart et al., 2018; Eun et al., 2018; Jun et al., 2018; Willis and Huang, 2017; Wood and Nurse, 2015; Zatulovskiy and Skotheim, 2020). The presence of these diverse phenomena raises the question as to why the underlying control networks evolve one rather than another type of cell size control?

To explore the evolution of cell size control networks subject to distinct selection pressures, we used computational evolution simulations. We initially examined the evolution of a seed cell cycle model consisting of G1 and S/G2/M phases of similar duration, where the G1 phase was free to evolve size-dependence, but the S/G2/M phase was constrained as a timer. Our simulations reliably evolved a control mechanism based on a Mixed Feedback Loop (François and Hakim, 2005, Figure 2). This network is centered on a cell cycle regulator (I) that inhibits the G1/S transition in proportion to its quantity. [I] is titrated away into an inactive complex by an increasing amount of another protein [R] that is synthesized in proportion to cell size. This results in a size-dependent decrease in the effective cell cycle inhibitor (free [I]). Thus, our evolved network implements an effective dilution of a cell cycle inhibitor that is conceptually similar to the well-described inhibitor dilution models of budding yeast, human cells, and Arabidopsis plants (D’Ario et al., 2021; Schmoller et al., 2015; Xie and Skotheim, 2020; Zatulovskiy et al., 2020). We note that we did not allow the synthesis of our proteins, such as [I], to be size-independent as has been found for budding yeast (Chen et al., 2020; Schmoller et al., 2015; Swaffer et al., 2021) as this could result in a one-step implementation of cell size control through the pure dilution of a cell cycle inhibitor. It is therefore interesting that given the constraint that all proteins be made in proportion to cell size, the network still evolved an effective ‘dilution’ of the active form of the cell cycle inhibitor molecule [I] .

Our evolution simulations gave insight into factors that bias evolution towards sizer or adder type control mechanisms (Figure 4). First, it is worth noting that our evolution simulations were not deterministic. There was no one-to-one correspondence between a given evolutionary pressure and any one specific cell size control mechanism. Rather, our claims represent an average behavior observed over the course of many simulations. Size control, as measured by the CV at a particular point in the cell cycle, has contribution both from the slope of the correlation between cell size and the amount of cell growth, and from the amount of noise characterizing the differences between cells that are initially the same size (Di Talia et al., 2007). It is therefore possible that a low noise adder can produce a lower CV than a higher noise sizer. This is reflected in the evolutionary paths of some of our simulations, which traverse from a noisy sizer to a less noisy adder (Figure 5). However, we anticipate even noisy sizers will be better than adders at controlling cell size in response to large deviations away from the steady state distribution. This is because sizers will always return the cell size to be within the steady state distribution within a cell cycle. We note that these generic results of how sizers and adders can govern cell size homeostasis can be derived from more traditional analytical methods (Barber et al., 2017; Willis et al., 2020). However, our evolution simulations are particularly useful because the molecular networks that evolved give non-trivial insights into how the observed size homeostasis dynamics can be regulated (e.g. via a Mixed Feedback Loop or using a system close to criticality). They are also suggestive of evolutionary pathways: despite different evolutionary modalities and control types, a natural step in many of our simulated evolutions is a system with strong sizers at very small volume only (Figures 5 and 6). This is practically reflected in a strongly negative slope on the very left side of the ΔVCycle plots, and a positive slope at higher volume corresponding to timers (Figures 5B and 6E). Similar non-monotonicity of ΔVCycle has been identified in models of various realism and complexity (Chandler-Brown et al., 2017; Delarue et al., 2017) and we provide here an evolutionary explanation for such an effect. We thus predict that this will be observed in systems where CV at birth does not need to be tightly controlled.

In the selection of a size controlling G1 network followed by a timer in S/G2/M, we observed a prevalence of adders that is consistent with the prevalence of adders reported in the literature. While fewer in number, sizers have also been observed. That the most accurate sizers have been observed in the fission yeast S. pombe (Fantes, 1977; Sveiczer et al., 1996; Wood and Nurse, 2015), and that this organism performs cell size control at G2/M rather than at G1/S led us to explore the effect of cell cycle structure on the evolution of cell size control. We found that controlling cell size later in the cycle in S/G2/M biases evolution away from adders and towards sizers. In retrospect, this result can be rationalized since any size deviations incurred earlier during the timer period can be compensated for by the end of the cycle with the sizer. However, when the order is inverted, any size deviations escaping a G1 control mechanism would only be amplified by exponential volume growth during the S/G2/M timer period. A second recent case exhibiting sizer control was found in mouse epidermal stem cells, which exhibit a greatly elongated G1 phase and a relatively short S/G2/M phase (Mesa et al., 2018; Xie and Skotheim, 2020). We found that if we increased the relative duration of G1 in our simulations by shortening the S/G2/M timer, we also see a bias towards sizer control. In essence, by extending G1 to a larger and larger fraction of the cell cycle the control system is gradually approaching a size control taking place at the end of the cell cycle, that is, an S/G2/M size control. Taken together, these simulations suggest the principle that having size-dependent transitions later in the cell cycle selects for sizers, while having such transitions earlier selects for adders.

In addition to identifying cell cycle structural features that canalize evolution towards sizers and adders as described above, we also observed an intriguing mechanism relying on molecular fluctuations. In this case, small cells would trigger an abnormally long G1 that would result in very large cells that decrease in size through a series of rapid cell divisions. This type of size control is reminiscent of that found in the green algae Chlamydomonas where a series of rapid, size-reducing cell divisions cease when cells go below a target size (Heldt et al., 2020). In our case, small size results in larger concentration fluctuations due to Poisson noise in the number of molecules. These concentration fluctuations, when large enough, are then able to trigger a burst of G1/S inhibitor that leads to an extended G1 phase and massive cell size growth before another series of rapid cell divisions is initiated (Figure 6). Interestingly, the system thus performs statistical size control over many generations. Intriguingly, this size control mechanism exhibits hallmarks of self-organized criticality (SOC). Just like adding grains of sand to a pile eventually triggers avalanches, the consistent decrease of cell size in the rapid division cycles eventually triggers a greatly extended G1 phase. To our knowledge, this is the first example where self-organized criticality is obtained in artificially evolved models of gene networks performing a well-defined function and is consistent with the idea that evolution of complex systems can favor the emergence of critical processes.

Materials and methods

Mathematical formalism

Request a detailed protocol

To model gene networks, we follow a standard ODE based formalism, where we simulate dynamics of the concentrations of proteins. We use Hill functions for transcriptional interactions, and standard mass action kinetics for protein-protein interactions. We also assume that all proteins are degraded at a constant rate. For most of the simulations presented in the paper, we use deterministic ODEs for simplicity. Importantly, cell-to-cell variability arises from the precise timing of cell cycle progression events. This allows for a natural way to generate noise on cell volume that should then be compensated for by the evolved network. In the last part of the paper, we explicitly include Langevin noise for biochemical reactions that are modeled using a classical tau-leaping formalism (Gillespie, 2007). Thus, each biochemical reaction takes place with a rate that corresponds to the deterministic rate, to which we add one white Gaussian noise with a variance equal to that rate. For example, given a deterministic biochemical rate k and a time interval of size Δt, we consider a tau-leaping change of kΔt+N(0,kΔt) where N(0,kΔt) is a random gaussian variable of mean 0 and variance kΔt.

Volume influences protein dynamics in three ways. First, protein production rates are generally proportional to cell volume so that proteins reach and maintain a constant concentration that is independent of the cell volume (Chen et al., 2020; Elliott and McLaughlin, 1978; Newman et al., 2006; Swaffer et al., 2021). We note that we are not allowing the cell to employ proteins such as Whi5 in budding yeast whose production is independent of cell size so that its concentration is a direct readout of cell size (Schmoller et al., 2015; Swaffer et al., 2021). We chose to do this because we want to explore how cell size control can be done by a network with multiple feedbacks rather than just the concentration of a single protein with a special dedicated synthesis mechanism. Thus, the only deterministic influence of volume on concentration dynamics is on the dilution rate, which is proportional to the cell growth rate λ(V) (see details in Appendix 1). At cell division, we also assume that proteins are equally partitioned between the daughter cells, that is, the concentration is the same before and after division. Note that we scale all our variables so that a concentration of one arbitrary unit corresponds roughly to 1000 proteins in a 100fL cell (Milo et al., 2010). Additionally, we scale the time variable so that 1 arbitrary time unit corresponds roughly to 30 min (Di Talia et al., 2007).

In this study, we chose a hierarchical way of introducing noise in the system, starting with the biggest contributing factor and incrementally adding additional sources of noise in subsequent analyses. All simulations presented include noise (stochastic control of G1/S transition and timing of S/G2/M, see below) in the cell cycle phases, whose CV has been found to be as high as 50% (Di Talia et al., 2007). Then, we introduced protein production noise via Langevin noise because the CV of regulatory protein concentrations is typically 20–30% (Newman et al., 2006). Importantly, the cell volume also contributes to stochastic effects, which are larger in smaller cells with fewer molecules. Thus, for stochastic simulations, we include a multiplicative 1V contribution to the added Gaussian noise term (see more complete description in the Appendix 1).

We also checked that our results are largely invariant when adding other sources of noise (see Appendix 1—figures 57). In these simulations, we also included noise in cell growth rate (CV ~15%; e.g. Di Talia et al., 2007), and in mass partitioning at cytokinesis (CV ~10%; e.g. Zatulovskiy et al., 2020).

Evolutionary procedure

Request a detailed protocol

To evolve networks regulating cell size, we use the φ-evo software (Henry et al., 2018) with a modified numerical integrator accounting for volume dynamics and volume dependencies as described above (Figure 1B). φ-evo simulates the Darwinian evolution of a population of gene networks. A network is encoded with the help of a bipartite graph connecting biochemical species (typically proteins) and interactions between them (we use a custom-made Python library). Networks are converted into an ensemble of stochastic differential equations using a Python to C interpreter. This code is then compiled and integrated on the fly to compute the behavior of the networks.

Each selection step in the algorithm is referred to as an epoch rather than the more commonly used term generation because we use the term generation to refer to cell divisions in the simulations. At each epoch of the algorithm, each gene network is simulated, and its fitness computed. Based on the fitness function(s) (see below), half of the networks are selected and duplicated, while the other half is discarded to maintain a constant population size. The duplicated networks are then randomly mutated. From the most to least probable, mutations consist in random changes of parameters of the network, random removal of interactions, and random additions of interactions or new proteins. Absolute mutation rates are adjusted as a function of the number of evolutionary epochs so that all networks in a population are mutated on average once per epoch. This implements a numerical equivalent of the biological Drake’s rule that mutation rates adjust with genome size (Lynch, 2007). Practically, this prevents the known phenomenon of code-bloating in evolutionary simulations (Foster, 2001) and also means that the total number of epochs is a good proxy of the number of mutations (in random directions) needed to evolve the best networks. All of this is easily made with our customized Python library encoding networks. For more details on technical aspects and implementations of the φ-evo software, we refer the reader to Henry et al., 2018.

Realistic evolutionary processes select for multiple phenotypes in parallel. While trade-offs between those phenotypes are non-trivial, it has been observed that phenotypes typically define an evolutionarily Pareto front (Shoval et al., 2012; Warmflash et al., 2012). We thus perform network selection using a Pareto mode (Warmflash et al., 2012), in which two distinct fitness functions are computed. During the selection step, networks are first Pareto ranked. For example, consider two networks A and B and two fitness functions f1 and f2. f1A refers to the fitness of network A calculated with function f1. Assuming fitness functions are to be maximized, we say network A Pareto-dominates network B if both f1A > f1B and f2A ≥ f2B. Rank 1 networks are networks which are not dominated by any other networks, Rank 2 networks are networks dominated only by Rank 1 networks, and Rank 3 networks are only dominated by Rank 1 and Rank 2 networks and so on. The algorithm then selects half of the population of Rank 1 networks using a fitness sharing algorithm to maximize population diversity (see details in Warmflash et al., 2012). One advantage of Pareto selection is the increased flexibility of the evolutionary process. Multiple fitness functions can provide different optimization paths in parameter space, which prevents the selection process from getting stuck in a local optimum of a single fitness function. We also perform a few simulations with only one fitness function, in which case networks are simply ranked based on their fitness.

We impose two evolutionary selection pressures in the form of two fitness functions. The first fitness function is simply the number of cell divisions during a long period, which we call NDiv . This is consistent with the classical definition of fitness as optimizing the number of offspring and is to be maximized by the algorithm. The second fitness function is the coefficient of variation of the volume distribution at birth for those NDiv generations, which we call CVBirth and is to be minimized by the algorithm. This penalizes broad distributions of volume at birth, which are detrimental to cell size homeostasis, which is what we aim to examine here. We further imposed fitness penalties to prevent way too small or too big cells, see Appendix 1. There, we also study alternative fitness functions, such as least-square residual function to minimize volume variation about a target size, and the fitted slope of the amount of volume added at each cycle to be minimized to drive models toward being a sizer.

Varying cell cycle structure

Request a detailed protocol

We also ran evolutionary simulations with different cell cycle structures. For evolutionary simulations where the G1/S transition was controlled by the concentration of [I] , we simply change the probability to pass the G1/S transition to depend on concentration [I] instead of its quantity I. For evolutionary simulations with a cell cycle structure similar to that found in the fission yeast S. pombe, we invert the cell cycle network structure. In this case, I quantity controls division and the Switch is turned on for a fixed amount of time in G1. In terms of the relaxation oscillator, this means that the left branch is now S/G2/M and the right branch is G1.

Appendix 1

In section Mathematical implementation, we first describe how to modify deterministic and stochastic ODEs to account for a exponentially growing cell volume. In section Size control, we describe the three size control archetypes, namely the timer, the adder, and the sizer, as well as our implementation of the initial seed cell cycle model. Then, in section Evolutionary algorithm, we give details on the φ-evo evolutionary algorithm. In section Analysis of sources of noise, we present evolutionary simulations investigating the effects of noise on the evolved networks. Finally, in sections Model descriptions and Additional models, we provide parameter values and equations for the models presented in the main text and for additional models produced by our evolutionary simulations.

Mathematical implementation

Deterministic

All biochemical concentrations can be described as a quantity of molecules of a biochemical species divided by the volume that contains it. The fundamental equation describing all concentrations is thus:

(1) [X](t)=X(t)V(t)

Here, X(t) is the quantity of molecules of an arbitrary biochemical species X as a function of time and V(t) is the volume of the cell containing said species as a function of time. We will use the bracket notation [X](t) for concentrations. In this project, we will model all biochemical species directly at the concentration level and assume proteins are uniformly distributed in an exponentially growing cell volume. The absolute growth rate of the cell λ is chosen to be constant over a large viable range of cell volume and is otherwise 0. Consequently, volume grows over time with a fixed absolute growth rate λ=0.25 over a viable volume range following the equation:

(2) dVdt=λV(t)

Deterministic rate equations describing the dynamics of the biochemical species at the concentration level have to be adjusted to take into consideration a time-varying volume. From Equation 1 and the derivative chain rule, we get:

d[X]dt=ddtX(t)V(t)=1V(t)dXdt-X(t)V(t)2dVdt

which we can combine with Equation 2 to give:

(3) d[X]dt=1V(t)(dXdt|V=Cst)-λ[X](t)=f([X](t))-λ[X](t)

The first term in the Equation 3, f([X](t)), corresponds to the usual biochemical reaction rates that occur when the volume of the cell is fixed. The second term, -λ[X](t), is a dilution term that we can interpret as an effective degradation of the concentration [X](t) due to the exponential growth of the cell volume over time.

To further simplify the expression for f([X](t)), we make an assumption about the production rates of all biochemical species in our models. Let us consider the rate equation describing the dynamics of the quantity of an arbitrary protein X with generic production rate ρ and degradation rate δ contained in a fixed cell volume V.

dX(t)dt=ρ-δX(t)

Here, f([X](t)) is given by:

f([X](t))=1VdX(t)dt=ρV-δ[X](t)

We assume that all the proteins in our models are constitutively expressed by the cell. In other words, the production rates ρ are linear functions of the cell volume ρ=ρ(V(t))=ρ0V(t). This ensures that the protein production rates scale with the volume such that concentrations stay constant over time which is a general feature of most proteins in S. cerevisiae (Chen et al., 2020; Newman et al., 2006; Swaffer et al., 2021). This yields: f([X](t))=ρ0-δ[X](t)

Taken altogether with Equation 3, the deterministic dynamics of a constitutively expressed arbitrary protein concentration [X](t) contained in an exponentially growing cell volume are given by:

(4) d[X]dt=ρ0-(δ+λ)[X](t)

Stochastic

To simulate molecular noise, we follow a classical tau-leaping formalism (Gillespie, 2007). Specifically, we choose the Euler-Maruyama implementation to generate approximate solutions to stochastic differential equations (Kloeden and Platen, 1992).

As we have done before in the deterministic case, let us first consider the quantity of an arbitrary protein X(t) in order to extract the equation for the concentration [X](t). Let’s assume X(t) is changing via a single biochemical reaction rate dXdt=g(X(t)) over the time interval [0,T] given X(t=0)=X0. We will show later how this approach can be generalized to include multiple reaction rates. We begin by partitioning the time interval in N equal segments of length Δt such that 0<t1<t2<...<tN=T with tn=nΔt, n{1,N}.

The Euler-Murayama approximate solution to the stochastic differential equation at the discrete time points tn is then recursively given by the following equation for n{1,N-1} where the single biochemical reaction is assumed to happen with a Poisson rate (corresponding to the deterministic rate), which adds one white Gaussian noise to the differential equations with a variance equal to that rate:

(5) Xn+1=Xn+g(Xn)Δt+|g(Xn)|ΔtN(0,1)

Here, Xn=X(tn), g(Xn) is the drift term, |g(Xn)| is the diffusion term and N(0,1) is a random Gaussian variable of mean 0 and variance 1. In other words, Xn+1 is a random Gaussian variable of mean Xn+g(Xn)Δt and variance |g(Xn)|Δt. Since this describes the quantity of proteins, we also have to consider the change in volume over time to recover the equation for the concentration of proteins, [X](t). Thus, let’s consider the volume of the cell V(t) at the discrete time points tn, which is given recursively by the equation:

(6) Vn+1=Vn+λVnΔt=Vn(1+λΔt)

Here, we define Vn=V(tn) and recover the λVn term from Equation 2. We assume that the volume time evolution is noiseless for simplicity. To recover, the differential equation describing the protein concentration, we evaluate the expression [X]n+1=Xn+1Vn+1. Thus, combining Equations 5 and 6, we get:

[X]n+1=Xn+g(Xn)Δt+|g(Xn)|ΔtN(0,1)Vn(1+λΔt)[X]n+1=11+λΔt([X]n+g([X]n)Δt+|g([X]n)|ΔtVnN(0,1))

where we identify 1Vng(Xn)=g([X]n) as the rate equation describing the concentration of the consitutively expressed protein [X](t) when the volume V is fixed as described in the deterministic case. Then, we compute the derivative as:

d[X]dt|tn[X]n+1[X]nΔt=1Δt(1+λΔt)([X]n(1+λΔt)[X]n+g([X]n)Δt+|g([X]n)|ΔtVnN(0,1))=11+λΔt(g([X]n)λ[X]n+|g([X]n)|VnΔtN(0,1))

Finally, we recover the approximate full differential equation for the protein concentration by expanding the prefactor in the last equation to the 0th order in Δt assuming it to be small to give:

(7) d[X]dt|tng([X]n)λ[X]n+|g([X]n)|VnΔtN(0,1)

So far, we have assumed that there is only a single biochemical reaction rate g([X](t)). We can however easily generalize our approach to include additional reaction rates by summing up the contribution of each rate to the total differential equation. Given M independent reaction rates gi([X](t)) and the properties of random Gaussian variables, we can easily generalize:

(8) d[X]dt|tni=1Mgi([X]n)λ[X]n+i=1M(|gi([X]n)|VnΔtNi(0,1))

Importantly, the Ni are Gaussian vectors accounting for the noise correlations associated with single reactions. For instance, imagine one protein Xk turns into another protein Xp, then the corresponding Gaussian vector for this interaction takes the form N(0,1)(xkxp) where xk is vector of length corresponding with the number of variables in the system whose k-th component is equal to 1 with 0s elsewhere. This indicates that the molecular fluctuation due to this reaction should have opposite signs for Xk and Xp as expected.

The term λ[X]n is a dilution term that corresponds to an effective degradation of protein concentration [X]n as seen in the deterministic case. Interestingly, we highlight the 1Vn dependency in the noise term. We can understand this dependency intuitively by considering a protein production process with a Poisson parameter θ. In this scenario, the mean and the variance of the protein quantity distribution is given by the parameter θ. Going back to concentration space, there are an infinite number of combinations of protein quantity and volume that can give the same concentration [X]=X/V. Thus, we need to specify both the protein number X and the volume V to correctly model the molecular noise contributing to fluctuations in concentrations.

Size control

Initial seed network

To guide the evolutionary process, we begin with an initial seed network. We base our first seed network on the phenomenology of the budding yeast S. cerevisiae’s cell cycle, where cell size primarily regulates the timing of the START transition in late G1 (Schmoller et al., 2015). This regulation allows small daughter cells to delay the G1/S transition allowing them to catch-up in size by extending the G1 phase. S/G2/M duration on the other hand is largely independent of cell size. We note that while budding yeast divide asymmetrically, our simulated cells divide symmetrically. In our simple initial seed network, the cell cycle consists of two phases, G1 and S/G2/M, which respectively denote the pre-G1/S and post-G1/S phases of the cell cycle. The transition between these two phases is controlled by the level of a transcription regulator we call I. Like the Whi5 protein in S. cerevisiae, I is an inhibitor of the G1/S transition such that the lower its level, the higher the chances of cell cycle progression. Since protein production rates were assumed to be dependent on volume (as described in section Mathematical implementation), we found that I’s concentration alone was largely independent of volume and could not trigger a size-dependent G1/S transition as Whi5 does in budding yeast. Thus, we chose the quantity of I defined as I(t)=[I](t)×V(t) as the control variable for this transition. We chose to model the probability of the G1/S transition occurring at the next time point of the simulation with a sigmoid-shaped curve given by Equation 9 that can be visualized in Appendix 1—figure 1. We maintain θ=0.8 and nθ=8 fixed throughout this project. We chose these values because they give a similar amount of noise in the G1/S transition as observed experimentally (Di Talia et al., 2007; Chandler-Brown et al., 2017).

(9) PG1/S(t)=21+exp((I(t)/θ)nθ)
Appendix 1—figure 1
Probability of the G1/S transition occurring at the next time step.

X-coordinate is the quantity of the transcriptional regulator I. Y-coordinate is the probability of the G1/S transition occurring at the next time step Δt. Parameters θ=0.8 and nθ=8.

In our seed model, we encode cell cycle state using a binary variable S/G2/M Switch, which is 0 in G1 and turns to 1 in S/G2/M once the G1/S transition takes place. Following S. cerevisiae’s cell cycle structure where S/G2/M duration is independent of cell size, we fix S/G2/M duration to be 50% of the doubling time τ=1λln2 with uniform noise unless stated otherwise. This way, cells can tune the length of their cell cycle by adjusting G1 length while being constrained by the incompressible length of the timer in S/G2/M, TS/G2/M.

Following S/G2/M, cells divide such that VBirth(n+1)=VDivision(n)/f with f the division fraction. The n exponent here is referencing the n-th generation in the cell lineage. We choose f=2 for all simulations performed in this study unless explicitly mentioned otherwise. We assume perfect partitioning of all proteins between the two daughter cells such that the proteins’ concentrations remain the same before and after division. After division, we follow one of the two daughter cells during their own subsequent cycle. If we simulate the cell lineage for a long time, ergodicity guarantees that all volume states will be visited given stable growth and we can extract population statistics from the lineage data itself. Here, cell growth is exponential on the single cell level since we were assessing size control mechanisms that take size as an input to cell cycle control. We are not exploring the very interesting case where growth deviates from the exponential. In that case, size homoeostasis would have a contribution from some cells in the population outcompeting others in terms of their growth and we would have to simulate the entire cell population and not disregard one of the daughter cells as we do here.

Inspired by the dynamics of Whi5, which is produced in S/G2/M and diluted in G1, we chose an initial seed network where S/G2/M Switch activates the transcription of the I inhibitor. This ensures that the concentration [I] is ‘reset’ to a higher value following S/G2/M and prevents cells from skipping entirely the G1 phase of the subsequent cycle which would quickly send the volume of the cell converging quickly towards 0. We note that this interaction systematically appeared anyway in our early evolution simulation so we chose to include it in the initial seed network to accelerate the evolutionary process.

Size control archetypes

Size control mechanisms are often compared to three well-characterized models or archetypes in order to quantify the strength of the size control mechanism under study. Specifically, there are timers, adders and sizers. In this subsection, we will define each archetype and show that we can summarize them via a control volume response curve T(VC) as described in the main text.

First, let’s consider the timescale of growth. Given V(t=0)=V0, the solution V(t) of the volume Equation 2 is V(t)=V0eλt. From this equation, we can easily recover the doubling time τ defined as the time required to double a cell’s volume (V(τ)=2V0) which yields:

(10) τ=1λln2

The doubling time τ is only a function of the growth rate λ. Since we are only considering symmetrical division events, fixed interdivision times shorter than the doubling time will yield progressively smaller daughter cells. Similarly, fixed interdivision times longer than the doubling time will yield progressively larger daughter cells. Thus, the absolute growth rate λ sets the timescale for cell cycle dynamics if we want to simulate a stable cell lineage.

With the quantity sensing of the inhibitor I at the G1/S transition (see Equation 9), we find that the instantaneous volume at G1/S sets the concentration of the biochemical species for the rest of the cycle and until the next G1/S transition in the daughter’s cell cycle. Consequently, it regulates the timing of the G1 phase of the daughter cell and thus creates a return map for the volume at G1/S. We find that the volume at G1/S at the (n+1)-th generation VG1/S(n+1) is given recursively by:

VG1/S(n+1)=VG1/S(n)eλ(TS/G2/M+TG1(VG1/S(n)))2=VG1/S(n)eλT(VG1/S(n))2

To study the mechanisms of cell size control, we choose to define a useful new variable: the control volume VC. This control variable is independent from the biochemical network and maintained fixed allowing us to break the size feedback, and distinguish its input, the volume V, from its output, the induced cycle period T. With this new variable, we can modify the control at G1/S by forcing the transitions to trigger once IC=[I]×VC is low enough. We can then extract the response curve of the system, that is, the cell cycle period induced from sensing this control volume at G1/S T(VC). We represent this process schematically in Appendix 1—figure 2A.

The control volume at which the response curve T(VC) is equal to the doubling time τ corresponds to the equilibrium volume, Veq, for this network. Indeed, if T(Veq)=τ, then the cell cycle length will ensure that this cell exactly doubles its volume during its cell cycle and returns to the same Veq at the next generation. This volume is a fixed point of the volume return map and can be either stable or unstable. Theoretically, there could be size control mechanisms with multiple fixed points of the response curve, but practically we have not seen this emerge from any of our evolution experiments and therefore assume that Veq is unique. In the main text, we have substituted Veq by the average value of the volume at the time where volume is sensed as both of these values are essentially identical. This corresponds to VG1/S for models with a sizer in G1 and a timer in S/G2/M and VDivision for networks with a timer in G1 and a sizer in S/G2/M.

Size variation naturally occurs in our models due to the precise timing of G1 and S/G2/M cycle phases which are both noisy, so the volume does not stay at Veq for very long. The stability of growth around this equilibrium volume however will depend on the sign of the local derivative with respect to control volume of the T(VC) response curve and we will consider the following three cases:

  • is strictly increasing with VC. In this case, small deviations around Veq are amplified over successive generations and the volume quickly shrinks to 0 or explodes to . In this case, we say that Veq is an unstable fixed point of the response curve.

  • T(VC) is constant with VC. In this special case and assuming exponential growth of the volume, the only stable mode of growth corresponds to the response curve T(VC)=τ. This corresponds to the only stable timer archetype. In this particular scenario, Veq is not well defined as there are an infinite number of volumes where the response curve intersects the doubling time.

  • T(VC) is strictly decreasing with VC. In this case, cells correct for size deviations over successive generations and perform size control. In this case, we say that Veq is a stable fixed point of the response curve.

Here, we found the T(VC) curves that were selected by the evolutionary algorithm were all decreasing with VC as expected for stable size control mechanisms. Thus, for the remainder of this document, we will assume that Veq is uniquely defined and corresponds to a stable fixed point of the control volume response curve.

Timer

The timer archetype describes mechanisms that monitor time rather than size. If the cycle duration of the timer is tuned precisely to the doubling time τ=1λln2, cells will double their mass over the course the cell cycle to ensure that newborn daughter cells have the same volume at birth as their mothers did when they were born. Consequently, this category of mechanisms is notoriously bad at correcting for size deviations given exponential cell volume. If growth was linear however, this mechanism would allow for size control to take place.

(11) TTimer(VC)=τ=1λln2
Adder

The adder archetype describes mechanisms where cells add a constant amount of cell volume during each cycle. We define Δ the increment of added volume between birth and division, that is Δ=VDivision-VBirth. For adder mechanisms, the added volume at each cycle is constant and does not depend on cell size. In this case, initial size deviations are reduced by a factor of 2 at each division such that the volume at birth geometrically converges to the added volume Δ over successive generations. To recover the adder response curve TAdder(VC), we consider that by definition, VDivision=VBirtheλTAdder. From this equation and the definition of the adder, we can recover the cycle period TAdder:

TAdder(VBirth)=1λln(ΔVBirth+1)

We would like write this equation as a function of the control volume VC at G1/S and the equilibrium volume Veq alone. Assuming that the G1/S transition is followed by a timer in S/G2/M, we can write VBirth=VCeλTS/G2/M/2. Similarly, since we know by definition that TAdder(Veq)=τ, we can recover that the added volume increment is Δ=VeqeλTS/G2/M/2. Finally, we can combine these two expressions with to recover the final expression of the adder response curve:

(12) TAdder(VC)=1λln(VeqVC+1)

We note here that if the control volume was measured at division or at birth, the response curve of the adder would be unchanged. The only difference would be that both VC and Veq would correspond to volumes at division or birth volumes instead of volume at the G1/S transition. Here for example, the volume increment Δ corresponds to the Veq at birth.

Sizer

The sizer archetype describes mechanisms that measure size directly and allow a cell to return to a target volume VTarget after a single generation, irrespective of how big or small a cell was initially. For sizers, VDivision=VTarget=VBirtheλTSizer by definition.

We can then extract:

TSizer(VBirth)=1λln(VTargetVBirth)

We can then write VBirth and VTarget as a function of the control volume at G1/S VC and the equilibrium volume Veq. Using the same definitions as before VBirth=VCeλTS/G2/M/2 and TSizer(Veq)=τ, we find that Vtarget=VeqeλTS/G2/M. The sizer response curve thus follows:

(13) TSizer(VC)=1λln(2VeqVC)

We note again here that if the volume was measured at division or at birth, the equation for the sizer would be identical with the only difference being that the control volume VC and the equilibrium volumes Veq would correspond to division or birth volumes respectively.

Appendix 1—figure 2
Control volume and archetype response curves.

(A) Schematic representation of the way we break the feedback in the system and impose a control volume VC (red arrow) at the G1/S transition in order to record the induced cell cycle period T(VC). (B) Response curve of the 3 size control archetypes. X-coordinate is the control volume at G1/S VC normalized by the equilibrium volume Veq. Y-coordinate is the response curve of the models T(VC) normalized by the doubling time τ. The dark blue curve is the response curve for the timer of length τ, the orange curve the response curve for the adder, and the red curve the response curve for the sizer. The dotted grey line indicates the equilibrium volume Veq. The shaded region corresponds to the region where growth is unstable and volume diverges over successive generations.

The three archetypes’ response curves are shown in Appendix 1—figure 2B. From these curves and given the particular cell cycle structures we examined, we can extract multiple relevant measures of size control such as the volume at birth, G1/S, and division from which we get the added volume during each phase of the cell cycle. It is noteworthy that the derivative of the added volumes ΔV with respect to the birth volume VBirth for the timers, adders, and sizers, are respectively 1, 0, and –1. Size control mechanisms are typically compared to the 3 archetypes by measuring the amount of added cell volumes over their cell cycles ΔVCycle and then fitting a linear model to these data points. The fitted slope of the linear model then informs what archetype this particular mechanism is more akin to. Some models evolved with added volume slopes lower than –1 and we call those super-sizers. Such mechanisms overcompensate for volume deviations about the equilibrium value which can increase variation in the size distribution instead of decreasing it.

Response curve of the seed networks

Appendix 1—figure 3
Response curves for the initial seed networks.

Columns indicate the size scaling assumption of the protein production rates as indicated above the figure. Rows indicate quantity or concentration sensing of inhibitor I at the G1/S transition assumption as indicated on the left side of the figure. In each panel, we first show the seed network’s response curve T(VC) as a function of control volume VC at G1/S. Sizer (red), adder (orange) and timer (dark blue) archetypes are shown for comparison. Second, we provide a schematic representation of the [I] trajectory in G1. Schematic trajectories are shown for low volumes (light pink) and high volumes (dark red). (A) Quantity sensing of I at G1/S with a size-dependent production rate in S/G2/M. Here, cells are born with a constant concentration [I]. Because of the quantity sensing at G1/S, the concentration of [I] at G1/S scales as 1/V. Thus, the time spent in G1 scales with V. This is the initial seed network we chose for most of our evolution experiments. (B) Quantity sensing of I at G1/S with a size-independent production rate in S/G2/M. Here, cells are born with a concentration [I] at birth that scales as 1/V. Because of the quantity sensing at G1/S, we again find that the concentration of [I] at G1/S scales as 1/V. Thus, the time spent in G1 is constant. (C) Concentration sensing of [I] at G1/S with a size-dependent production rate in S/G2/M. Here, cells are born with a constant concentration [I]. Because of the concentration sensing at G1/S, we find that the concentration of [I] at G1/S is constant. Thus, the time spent in G1 is constant. (D) Concentration sensing of [I] at G1/S with a size-independent production rate in S/G2/M. Here, cells are born with a concentration [I] that scales as 1/V. Because of the concentration sensing at G1/S, we find that the concentration of [I] at G1/S is constant. Thus, the time spent in G1 scales as 1/V.

In light of the control volume and response curve definitions from subsection 2.2, we can revisit the initial seed model and investigate how different assumptions alter the stability of growth and division in a cell lineage. Specifically, we investigate the size scaling assumption of the protein production rates and the concentration vs. quantity sensing of the transcriptional regulator I at G1/S as previously described in Section 1. We summarize our results in Appendix 1—figure 3.

First, let us consider the G1 trajectory of the transcriptional regulator [I] who is solely produced during the S/G2/M timer. The dynamics of [I](t) in G1 will be described by the following equation:

(14) [I](t)|G1=[I]0e(δ+λ)t

Here, the time variable t represents the time since birth, [I]0=[I](t=0), δ is the protein’s degradation rate and λ is the growth rate of the cell volume. This equation holds until the G1/S transition where the S/G2/M Switch is turned on again and G1 ends. Because the degradation of the inhibitor does not yet depend on volume in any way, the time spent in G1 will only be dependent on the ratio between: (1) the initial condition at birth [I]0; (2) the final condition at the G1/S transition, [I]G1/S.

We found that the size scaling assumption of the protein production rates influences the initial condition at birth [I]0. When production rates scale with size, we find that [I]0 is independent of volume. This is expected as this assumption was chosen specifically to model proteins whose concentrations are independent of size. Conversely, when we modify this assumption and consider that protein production rates are independent of size (e.g. like Whi5 in budding yeast), we find that the system produces a constant quantity of inhibitor instead of a constant concentration. This means that the initial concentration at birth [I]0 scales as 1/V.

Similarly, when imposing quantity sensing of I at G1/S, we found that the concentration of inhibitor at G1/S [I]G1/S scales as 1/V. Finally, when imposing concentration sensing of I at G1/S, we found a constant concentration of inhibitor at G1/S [I]G1/S as was expected by design.

Together, those assumptions alter the scaling of the duration of the G1 phase of the initial seed cycle. We summarize these results and present the models’ response curves T(VC) in Appendix 1—figure 3 where each row and column corresponds to a specific combination of assumptions. There, in Appendix 1—figure 3A, we see that for the combination of size-scaling production rate and quantity sensing at G1/S, we get a cell cycle period that is increasing with control volume VC. This is undesirable and leads to unstable growth of the cell lineage towards 0 or , but rewards the evolution of size control mechanisms that can prevent this unstable growth. We chose this initial seed model for most of our evolutionary simulations. In Appendix 1—figure 3B,C, we found that the two assumptions compensated each other to create size-independent timer models. The parameters of the network can be precisely fine-tuned to yield a response period of exactly τ as was done to produce the response curves shown here. Thus, it is technically possible to evolve a size control mechanism using these initial seed models, but we chose not to go down that path because we wanted to evolve an active size control mechanism. Finally, in Appendix 1—figure 3D, we see that if we assume that protein production rates do not scale with size and that the G1/S transition depends on the concentration of inhibitor [I], we get an initial seed model that already accomplishes size control as it displays a response curve that decreases with control volume VC. This simple model loosely corresponds to the Whi5 inhibitor dilution model of budding yeast (Schmoller et al., 2015) where a constant quantity of inhibitor Whi5 is present at birth (and thus a concentration [Whi5]1/V) and is passively diluted in G1 until it reaches concentration threshold that triggers the G1/S transition.

Evolutionary algorithm

Here we briefly describe the φ-evo evolutionary algorithm from Henry et al., 2018 that we used to evolve size control networks. We refer the reader to the original publication’s main text and supplementary material for a more thorough description of the algorithm. A schematic representation of the algorithm’s architecture is shown in Appendix 1—figure 4.

First, an initial seed network is selected by the user as the starting point of the evolution simulation. φ-evo then clones this first individual to create a population of networks. At each epoch, mutations are randomly applied to the networks of the population. Those mutations vary from topological changes to the network, where biochemical species or interactions can be added or removed, to non-topological changes, where the networks’ kinetic parameter values are modified. Following mutations, networks are ranked based on their performance at accomplishing the biological function we select for. This performance is encoded via a user-defined fitness function that is problem specific. We give details about the specific implementation of the fitness functions for cell size control in the following subsection. After ranking the networks, φ-evo proceeds to select the most fit half of the network population. The less fit half is then discarded and replaced by a copy of the most fit half to maintain a constant population size. With this, φ-evo completes the first epoch of the evolutionary process. We use the term epoch here rather than the term ’generation’ which we retain to describe a cell lineage. A predetermined number of epochs of mutation and selection are then performed after which a final population of networks is extracted.

Appendix 1—figure 4
Schematic representation of the φ-evo algorithm.

We begin with a user-defined initial seed network as starting point of the evolutionary process. The seed network is cloned to give a first population of networks. Individuals are then mutated randomly given the mutation parameters of the run. The dynamics and fitness scores of the networks are then computed and ranked. The best half of the population is selected and retained and the rest are discarded. The best half is then duplicated to maintain a constant population size N. We then repeat these instructions for a predefined number of epochs, after which a final population of networks is extracted and analyzed.

Fitness

In order to rank and select networks based on their performance at accomplishing a specific biological function, we design a specific objective function that we call fitness. Here, we chose a fitness function that could quantify a model’s ability to produce many viable descendants during a fixed time period of length t. We initially considered a simple fitness function to be minimized by φ-evo, f0=-NDiv. Here NDiv is the number of divisions or generations in a cell lineage produced during a total time period of length t. Noise at the G1/S transition and in the S/G2/M timer duration act as a source of variation in volume at each generation which needs to be controlled by the evolved networks in order to prevent the cell volumes from diverging. Thus, networks that perform size control display a high number of divisions NDiv. The cycle duration distribution of a size control network will be centered around the doubling time τ in order to promote stable growth. Thus, on average, we expect a fit network to exhibit a maximum number of max(NDiv)=t/τ divisions during a simulation of length t. Note that this number is mostly independent of the volume range selected by the evolutionary simulation as the doubling time τ is independent of the initial volume of the cell at the beginning of each cycle. There is a small effect on NDiv from the initial conditions chosen for the system of ODEs modelling the cell cycle, but this effect is mostly negligible as long as the total time period tτ.

In our first evolution experiments, we found that the single objective function f0=-NDiv was insufficient alone to evolve a cycling network. A possible reason for this is that the fitness landscape in parameter space defined by f0 is mostly flat far from the optimum and is difficult to navigate as it doesn’t incrementally guide the evolution process towards a proper size control phenotype. Indeed, a network that does not perform size control will display very few divisions before diverging towards sizes of 0 or . In contrast, a network that does performs some size control, even if performed badly, will display mostly stable growth with many divisions and the volume will not diverge over successive generations. There is thus an all or nothing effect with this fitness function. We found that optimization process would often get stuck in a local optima with a low number of divisions and could not find a path to the global optima of size control. Because of this, we chose to turn towards multi-objective Pareto optimization which aims to simultaneously optimize several fitness functions. The idea here is that an additional fitness function can guide the evolutionary process through a different path in parameter space and could allow the evolutionary procedure to escape local optima.

In the Pareto optimization framework, we consider N equally important objective functions f=(f0,f1,,fN). Assuming that fitness functions are to be maximized, we say that individual i (strictly) dominates individual j if and only if their fitness fi and fj are such that:

k{1,N},fkifkj and k{1,N}|fki>fkj

The algorithm then selects half of the population of the highest rank using a fitness sharing algorithm to maximize population diversity. We refer the reader to Warmflash et al., 2012 for more details on this procedure.

For this project, we chose to limit the optimization process to two fitness functions. We chose f0=-Ndiv as the first fitness function. We tested different measures of size control for the second objective function which are described in the following subsections. In most of the cases, we chose the coefficient of variation of the size distribution at birth CVBirth to be minimized as the second fitness function. In any case, we typically run 10 independent realizations of a network’s performance and compute the average fitness score over those runs to buffer variations in fitness scores.

Residuals

We first tested a least squared residuals fitness function to be minimized by φ-evo. This function yielded successful evolutionary runs but was abandoned due to requiring a user-defined target volume Vt. Indeed, we wanted to avoid the bias where we could select for biochemical networks matching a specific volume range. We used the following equation for the fitness function with VBirth(n) indicating the cell volume at birth at the n-th generation in a lineage.

(15) f1=1NDivVtn=1NDiv(VBirth(n)-Vt)2

We nevertheless present the result of a successful evolution run using Pareto optimization of NDiv and f1 in Appendix 1—figure 11 where we obtain a version of the feedback-based network topology of Model A1.

Coefficient of variation

We then considered the coefficient of variation of the size distribution at birth (CVBirth) to be minimized by φ-evo. The CVBirth is a measure of size control that normalizes the variance of the size distribution at birth with respect to its mean and is thus mostly insensitive to the absolute volume range of the cell. We chose this second fitness function in most of the Pareto optimization evolution experiments as described in the main text.

(16) f2=CVBirth=Var[VBirth]E[VBirth]
Added volume slope in G1

We also considered directly optimizing the fitted slope of the volume added in G1 as a function of volume at birth to reinforce the sizer behavior in G1. As described in the main text, given a series of volume values at birth and at G1/S, VBirth(i) and VG1/S(i) for i{1,NDiv}, the added volumes in G1 are defined as:

ΔVG1(i)=VG1/S(i)-VBirth(i)

The best fit slope m of a linear model ΔVG1=mVBirth+b has a closed-form equation which is given as a function of the lineage data directly:

(17) m=NDiv(i=1NDivVBirth(i)ΔVG1(i))-(i=1NDivVBirth(i))(i=1NDivΔVG1(i))NDivi=1NDiv(VBirth(i))2-(i=1NDivVBirth(i))2

A slope of –1 corresponds to a sizer, a slope of 0 to an adder and a slope of +1 corresponds to a timer. In order to directly optimize the size control mechanism in G1 and to bias towards sizer mechanisms and to keep the fitness values positive, we considered the following fitness function to be minimized by φ-evo:

(18) f3=m+1
Fitness penalties

In order to keep biochemical concentrations and cell volumes at reasonable levels, we chose to bound volume growth to a range V[Vmin,Vmax] with Vmin=0.1 and Vmax=100. To guide the evolution of size control and to have a well defined steady-state distribution of cell sizes, we chose to impose fitness penalties on networks that would see their volumes reach one of those bounds at some point during a run. This applies to all evolution experiments performed in this study, but we will describe the case of the NDiv and CVBirth fitness functions as they were used most of the time during this study.

Firstly, if the volume reached Vmin, we considered the cell too small and declared it dead. At that point, any further cycles would not contribute to fitness scores. With this penalty, networks are guided towards preventing or at least delaying the time at which the volume becomes too small for the cell to remain viable.

Similarly, when volume reaches Vmax, we penalize the fitness functions but in a different way. Because volume growth is restricted to the domain below Vmax, we set the growth rate λ=0 when V=Vmax. In this case, we sometimes see networks make use of this growth arrest and exploit this artificial feature. This phenomenon can sometimes be seen in computational evolution where digital mirages are often exploited by optimization processes (Lehman et al., 2020). Here, such exploitative networks have initially long cycle periods Tτ and can sometimes spend the majority of their cell cycle in this growth arrest phase. Over successive epochs, this period T gets shortened to get more and more divisions to take place during a run, shortening the time spent in growth arrest at each cycle. Eventually, this optimization leads to a particular type of model that has very sloppy size control but that scores highly with the fitness functions because of the artificial growth arrest. Such networks develop a timer with period Tτ and use the artificial growth arrest to buffer any volume variation incurred from late G1/S transition or noisy S/G2/M duration. These networks are optimal from an NDiv perspective since they exactly double their mass at each cycle and perform the same number of division cycles during a run as an actual size control network. They are also more than optimal from a CVBirth since their growth is always stopped at Vmax. Thus, without fail, their birth volume VBirth=Vmax/2 and there is no variation at all in the distribution. This size control illusion is a global optimum of the 2D-fitness space and must thus be heavily penalized to prevent the optimization from selecting this phenotype. Thus, when V=Vmax, we only count 70% of the divisions NDiv which is sufficient to distinguish this artificial phenotype from actual size control mechanisms. Additionally, we penalize the CVBirth score by adding to it a penalty of +10. This is significantly different from the usual range of coefficient of variations which lie between 0 and 0.5 typically, and prevents the optimization from selecting the artificial phenotype.

The introduction of these fitness penalties improved the convergence rate of the φ-evo algorithm significantly. Specifically, we believe the penalty on Vmax being somewhat less severe than the one for Vmin improved the convergence rate drastically. This is probably because the initial cell cycle model chosen for the evolution runs exhibits unstable growth and inevitably sees the cell volume grow to Vmax or shrink to Vmin rapidly as shown in Appendix 1—figure 3. Thus, the fitness landscape surrounding the initial cell cycle model is quite flat and is difficult to navigate from an optimization perspective. Gradual improvements to the NDiv fitness function at Vmax by spending less and less time in growth arrest seemed to have helped the optimization process. This guided the algorithm towards better size control models more often than if penalties were absent.

Overall, even if the penalties somewhat biased the evolutionary process into following a phenotypic trajectory, they improved the convergence rate of the evolutionary process dramatically to the point that we decided to keep them for all experiments.

Biochemical interactions

Many ‘inverse-approach’ approaches in systems biology have focused on purely transcriptional networks (Cotterell and Sharpe, 2010; François et al., 2007; Fujimoto et al., 2008; Ten Tusscher and Hogeweg, 2011;), because they are generic, easier to study and can efficiently describe many biological dynamics (Alon, 2007). In this project, we extend the biochemical interactions available for evolution: we not only model transcriptional activation and repression but also include complexation also known as protein-protein interaction (PPI), and assume there passive degradation. Adding PPIs is especially crucial because they are well known to lead to non-linear effects (Buchler and Cross, 2009; Buchler and Louis, 2008) allowing for the simple implementation of complex dynamics such as genetic oscillations (François and Hakim, 2005) observed e.g. in circadian clocks (François, 2005), and such non-linear effects indeed play crucial roles for control in our evolved model. The equations for these interactions are presented in this subsection.

We model transcriptional activation of arbitrary network species Y’s production from species X using the following Hill equation:

(19) ActivationXY([X],tX:Y,nX:Y)=([X]/tX:Y)nX:Y1+([X]/tX:Y)nX:Y

Similarly, repression of arbitrary species Y’s production from species X is modeled via the following Hill equation:

(20) RepressionXY([X],tX:Y,nX:Y)=11+([X]/tX:Y)nX:Y

In these equations, tX:Y is the threshold required for [X] to activate or repress [Y] at 50% of its capacity and nX:Y is the Hill coefficient. While both types of interactions are modeled using Hill equations, their combined effects on a network species’ dynamics is computed differently. Indeed, only the maximum of all the activations is accounted for whereas the inhibitions are multiplicative and all are accounted for. Additionally, we allow some species to be produced at a basal rate b independent of any activator which counts as an additional activation.

Altogether using an example, assuming multiple species X1,,Xq activate the production of species Z while multiples species Y1,,Yr repress it, and assuming that Z has a basal production rate bZ and a maximum production rate pZ, then the total contribution of these interactions to the ODE for the dynamics of [Z] is given by:

(21) d[Z]dt=max[pZmax[([X1]/tX1:Z)nX1:Z1+([X1]/tX1:Z)nX1:Z,,([Xq]/tXq:Z)nXq:Z1+([Xq]/tXq:Z)nXq:Z],bZ]i=1r11+([Yi]/tYi:Z)nYi:Z

We use the law of mass-action to model PPIs and passive degradation. Specifically, if arbitrary species X and Y interact together and form a complex Z given a forward rate kf and a backwards rate kb, then the contribution of these interactions to the ODE for the dynamics of the system will be given by:

(22) d[X]dt=d[Y]dt=-d[Z]dt=-kf[X][Y]+kb[Z]

Lastly, all network species are assumed to be degraded at a passive rate. Thus, if an arbitrary species X is solely degraded with a rate δX, then the dynamic equation for the dynamics of [X] will be given by:

(23) d[X]dt=-δX[X]

Analysis of sources of noise

As mentioned in the main text and in this Appendix, we chose a hierarchical way of introducing noise in the system, starting with the biggest contributing factor and incrementally adding additional sources of noise in subsequent analyses. We first included noise in the cell cycle phases, specifically in the timing of the G1/S transition and in the length of the S/G2/M phase. Then in the later parts, we introduced protein production noise modeled as Langevin noise.

In the simulations presented in the main text, we chose not to include noise in the growth rate and in the division ratio as the recorded noise level for in experiments for these measures is lower than that for the timing of the cell cycle and the protein concentration noise (Di Talia et al., 2007; Newman et al., 2006; Zatulovskiy et al., 2020). Nevertheless, those are crucial assumptions that we made that we chose to investigate in more details here.

In subsection S/G2/M noise, we investigate how the level of noise in S/G2/M affects the conclusions drawn in the main text. Then, in subsection Growth rate noise we do the same for noise in the growth rate and in subsection Division ratio noise for the division ratio.

Appendix 1—table 1 shows the values of CVBirth of the three models presented in the main text compared to the values reported in the literature for budding yeast (Di Talia et al., 2007), fission yeast (Sveiczer et al., 1996) and mouse epidermal stem cell grown in the animal (Xie and Skotheim, 2020).

Appendix 1—table 1
Coefficients of variation: models and experiments.
ModelsData
NameCVBirthCell typeTime of size measureCV
A10.098Haploid budding yeastBudding0.17
A20.095Haploid fission yeastFission0.06
B0.235Mouse epidermal stem cellBirth0.17

S/G2/M noise

Here, we perform similar evolution experiments to those reported in Figure 4 of the main text to examine the effect of modulating the noise in the S/G2/M timer. We thus perform three independent experiment where we set the CV in the timer period to 0%, 5%, and 8% corresponding to no, medium, and high noise respectively. For reference, the CV of the timer period in the control condition where Model A1 was evolved is 3%. Note that we maintain the average duration of the timer to be about half the time it takes to double the cell’s volume. Having specified the S/G2/M timer parameters and starting from the initial seed network of Model A1, we perform evolution and select networks as previously. We compare ensembles of 60 networks for each noise level, half of them evolved under the Pareto optimization of NDiv and CVBirth and the other half under the single objective optimization of NDiv. The results are shown in Appendix 1—figure 5.

Increasing the noise, progressively leads to a loss of the sizer signature and increases the CVBirth. This is likely because the fixed duration of S/G2/M allows the system to accurately reset protein concentrations for the subsequent cell cycle to promote accurate G1 control (Willis et al., 2020). Thus, an increasing level of noise becomes associated with a worse accuracy in the size control mechanism which leads to loss of the sizer signature and increased CVBirth as can be seen in Appendix 1—figure 5D,E. Other results from the main text remain unchanged.

Appendix 1—figure 5
S/G2/M noise analysis.

Summary statistics for evolutionary simulations each having 500 epochs. Model A1 was used as the initial seed network. 30 simulations were performed using Pareto optimization of the number of divisions (NDiv) and the CV of cell size at birth (CVBirth), are labeled Pareto and are shown in full colors. 30 more simulations were performed using only the number of divisions as the fitness function, are labeled NDiv and are shown in colored outlines only. Scatter plots show the coefficient of variation of the size distribution at birth (CVBirth, Y-coordinate) as a function of the fitted added volume slope over the whole cycle as a function of volume at birth (Slope ΔVCycle, X-coordinate) for the most fit models evolved during each of the 60 independent simulations. Horizontal box plots above the scatter plots in A-C display the distributions of the added volume slopes for the Pareto and NDiv simulations. Timer (dark blue), adder (orange) and sizer (red) slopes are shown respectively at 1, 0, and –1 for comparison. Vertical box plots on the right of the scatter plots in A-C show the distributions of CVBirth for the Pareto and NDiv simulations. Asterisks represent p-values for the Welch’s t-Test between the distributions. For reference, NS indicates p>0.05, * indicates p<0.05, ** indicates p<102, *** indicates p<103 and **** indicates p<104. The values of CVBirth and Slope ΔVCycle for the initial seed Model A1 are shown as a black square in the scatter plot or as a dashed black line in the box plots. Each panel explores different S/G2/M noise levels. (A) Evolution results for no noise in S/G2/M duration. (B) Evolution results for a noise level in S/G2/M duration equal to 5%. (C) Evolution results for a noise level in S/G2/M duration equal to 8%. (D) Evolved Slope ΔVCycle distributions as a function of noise level in S/G2/M. For reference, noise level for the Control experiment from Figure 4 corresponds to 3%. Box plots in D-E represent the distributions for both the Pareto and NDiv evolution experiments. Here, increased S/G2/M noise leads to loss of the sizer signature. (E) Evolved CVBirth distributions as a function of noise level in S/G2/M. Here, increased S/G2/M noise leads to increased variability in the cell size distributions at birth.

Growth rate noise

Here, we perform similar evolution experiments as we did for the noise in S/G2/M but this time by adding noise in the growth rate λ. Specifically, at each generation of a cell’s lineage, we sample a growth rate from a Gaussian distribution centered around 0.25, the initial value we used in the rest of this project. We perform three evolution experiments with coefficient of variations for the growth rate distributions set to 3%, 5% and 8% corresponding to low, medium and high noise respectively. We perform 30 independent evolution runs with the Pareto optimization framework for each λ noise level, each of them starting from the initial seed network of Model A1. The results are shown in Appendix 1—figure 6.

We note that on average, the growth rate will remain centered around the same value, thus not affecting the optimal fitness score networks are able to achieve. However, individual variations at each cycle perturb the ability of the system of accomplishing size control by always modifying the doubling time τ=ln(2)/λ. This leads to progressive loss of the sizer signature and also increases the CVBirth. This increased noise in the system sometimes sends the cell volume towards 0 or as volume is kicked outside of the control mechanism’s working range. Since this behavior is highly penalized in the fitness functions scoring, the evolution finds a way to prevent this from happening via different strategies. Interestingly, at higher noise levels, strong sizers can still evolve but are not the most common phenotype. Instead, evolution seems to favor timers that reliably ensure timely cell division generation after generation. There is however a trade-off and these models exhibit a higher CVBirth due to the lower amount of size control. Adders can also be evolved at all tested noise levels and provide good size control with reliably low CVBirth.

Appendix 1—figure 6
Growth rate noise analysis.

Summary statistics for evolutionary simulations each having 500 epochs. Model A1 was used as the initial seed network. Only 30 simulations were performed using Pareto optimization of the number of divisions (NDiv) and the CV of cell size at birth (CVBirth), are labeled Pareto. Scatter plots show the coefficient of variation of the size distribution at birth (CVBirth, Y-coordinate) as a function of the fitted added volume slope over the whole cycle as a function of volume at birth (Slope ΔVCycle, X-coordinate) for the most fit models evolved during each of the 30 independent simulations. Horizontal box plots above the scatter plots in A-C display the distributions of the added volume slopes. Timer (dark blue), adder (orange) and sizer (red) slopes are shown respectively at 1, 0, and –1 for comparison. Vertical box plots on the right of the scatter plots in A-C show the distributions of CVBirth. Asterisks represent p-values for the Welch’s t-Test between the distributions. For reference, NS indicates p>0.05, * indicates p<0.05, ** indicates p<102, *** indicates p<103 and **** indicates p<104. The values of CVBirth and Slope ΔVCycle for the initial seed Model A1 are shown as a black square in the scatter plot or as a dashed black line in the box plots. Each panel explores different growth rate noise levels. (A) Evolution results for low noise in growth rate with associated coefficient of variation at 3%. (B) Evolution results for medium noise in growth rate with associated coefficient of variation at 5%. (C) Evolution results for high noise in growth rate with associated coefficient of variation at 8% (D) Evolved Slope ΔVCycle distributions as a function of noise level in the growth rate. For reference, noise level for the Control experiment from Figure 4 corresponds to no noise. Here, increased growth rate noise leads to rapid loss of the sizer signature. (E) Evolved CVBirth distributions as a function of noise level in the growth rate. Here, increased growth rate noise leads to increased variability in the cell size distributions at birth.

Division ratio noise

Here, we perform similar evolution experiments as we did for the noise in S/G2/M and in λ, but this time by adding noise in the division fraction f. Specifically, at each generation of a cell’s lineage, we sample f from a Gaussian distribution centered around 2, the initial value we used in the rest of this project for symmetrical divisions. We perform three evolution experiments with coefficient of variations for the growth rate distributions set to 2%, 4% and 8% corresponding to low, medium and high noise respectively. We perform 30 independent evolution runs with the Pareto optimization framework for each f noise level, each of them starting from the initial seed network of Model A1. The results are shown in Appendix 1—figure 7.

The results of this experiment are very similar to those for noise in the growth rate λ described in the previous subsection. Indeed, changing the division fraction f does not change the doubling time τ directly like for λ. Instead, it changes the cycle period around which cells see their volume shrink or grow over successive generations. Indeed, for a division fraction f, this equilibrium time between shrinking and growth becomes τf=ln(f)/λ. Intuitively, if f is bigger than 2, then cells need to spend a little bit more time in their cell cycles for growth to occur in order to compensate for this increased division fraction. This increased noise in f leads to progressive loss of the sizer signature as seen before and also increases the CVBirth. The division fraction noise has a big effect on premature cell death. Indeed, at the highest noise level tested where CVf=8%, we saw many runs end with premature cell death to the point where three evolution runs were completely unable to find a mechanism able to prevent this. As a result of this strong pressure to avoid premature cell death, evolution turns once again to timers instead of sizers in the noisier regime. As before, adders seem to be the most reliable size control phenotype exhibiting low CVBirth.

Appendix 1—figure 7
Division fraction noise analysis.

Summary statistics for evolutionary simulations each having 500 epochs. Model A1 was used as the initial seed network. Only 30 simulations were performed using Pareto optimization of the number of divisions (NDiv) and the CV of cell size at birth (CVBirth), are labeled Pareto. Scatter plots show the coefficient of variation of the size distribution at birth (CVBirth, Y-coordinate) as a function of the fitted added volume slope over the whole cycle as a function of volume at birth (Slope ΔVCycle, X-coordinate) for the most fit models evolved during each of the 30 independent simulations. Horizontal box plots above the scatter plots in A-C display the distributions of the added volume slopes. Timer (dark blue), adder (orange) and sizer (red) slopes are shown respectively at 1, 0, and –1 for comparison. Vertical box plots on the right of the scatter plots in A-C show the distributions of CVBirth. Asterisks represent p-values for the Welch’s t-Test between the distributions. For reference, NS indicates p>0.05, * indicates p<0.05, ** indicates p<102, *** indicates p<103 and **** indicates p<104. The values of CVBirth and Slope ΔVCycle for the initial seed Model A1 are shown as a black square in the scatter plot or as a dashed black line in the box plots. Each panel explores different division fraction f noise levels. (A) Evolution results for low noise in division fraction with associated coefficient of variation at 2%. (B) Evolution results for medium noise in division fraction with associated coefficient of variation at 4%. (C) Evolution results for high noise in division fraction with associated coefficient of variation at 8%. 27/30 evolution runs succeeded and are shown here. (D) Evolved Slope ΔVCycle distributions as a function of noise level in the division fraction. For reference, noise level for the Control experiment from Figure 4 corresponds to no noise. Here, increased noise leads to rapid loss of the sizer signature. (E) Evolved CVBirth distributions as a function of noise level in the division fraction. Here, increased noise leads to increased variability in the cell size distributions at birth.

Model descriptions

In this section, we give the full set of equations and parameter values of the models from the main text. We remind the reader that we scale all our variables so that a concentration of one arbitrary unit corresponds roughly to 1000 proteins in a 100fL cell (Milo et al., 2010). Additionally, we scale the time variable such that 1 arbitrary time unit corresponds roughly to 30 min (Di Talia et al., 2007).

Model A1

The parameter values of the Model A1 are shown in Appendix 1—table 2 along with the corresponding differential equations in Equation 24.

(24) d[I]dt=p2(S/G2/M Switch/t1:2)n1:21+(S/G2/M Switch/t1:2)n1:2-kf1[I][R]+kb1[S4]-kf2[I]2+kb2[S5]-(δ2+λ)[I]d[R]dt=p311+([I]/t2:3)n2:3-kf1[I][R]+kb1[S4]-(δ3+λ)[I]d[S4]dt=kf1[I][R]-kb1[S4]-(δ4+λ)[S4]d[S5]dt=kf2[I]2-kb2[S5]-(δ5+λ)[S5]
Appendix 1—table 2
Model A1 parameter values.
ParameterValueParameterValue
p20.369408δ30.4574764
t1:21.104619kf16.783001
n1:23.215727kb10.195168
δ20.083827δ44.244007
p30.703658kf20.021174
t2:30.044938kb21.075146
n2:33.266732δ51.1010876

Model A2

The parameter values of the Model A2 are shown in Appendix 1—table 3 along with the corresponding differential equations in Equation 25.

(25) d[I]dt=p2(S/G2/MSwitch/t1:2)n1:21+(S/G2/MSwitch/t1:2)n1:211+([R]/t3:2)n3:2kf1[I][R]+kb1[S4]kf2[I]2+kb2[S5](δ2+λ)[I]d[R]dt=p311+([I]/t2:3)n2:311+(S/G2/MSwitch/t1:3)n1:3kf1[I][R]+kb1[S4](δ3+λ)[I]d[S4]dt=kf1[I][R]kb1[S4](δ4+λ)[S4]
Appendix 1—table 3
Model A2 parameter values.
ParameterValueParameterValue
p21.915601n1:32.751652
t1:20.17872t2:30.441106
n1:22.09054n2:32.780297
t3:20.962612δ30.045051
t3:22.144666kf10.839879
δ20.019495kb10.381067
p31.944803δ40.913992
t1:30.4229390

Model B

The parameter values of the fluctuation-sensing Model B are shown in Appendix 1—table 4 along with the corresponding stochastic differential equations in Equation 26. We note that the C0 parameter in the noise terms of the equation is a concentration scaling factor.

(26) d[I]dt=p2max((S/G2/M Switch/t1:2)n1:21+(S/G2/M Switch/t1:2)n1:2,([A]/t3:2)n3:21+([A]/t3:2)n3:2)(δ2+λ)[I]+N1(0,1)p2max((S/G2/M Switch/t1:2)n1:21+(S/G2/M Switch/t1:2)n1:2,([A]/t3:2)n3:21+([A]/t3:2)n3:2)+δ2[I]VC0Δtd[A]dt=p3([A]/t3:3)n3:31+([A]/t3:3)n3:3kf[A][S4]+kb[S5](δ3+λ)[I]+N2(0,1)p3([A]/t3:3)n3:31+([A]/t3:3)n3:3+δ3[A]VC0ΔtN3(0,1)kf[A][S4]VC0Δt+N4(0,1)kb[S5]VC0Δtd[S4]dt=p4kf[A][S4]+kb[S5](δ4+λ)[S4]N3(0,1)kf[A][S4]VC0Δt+N4(0,1)kb[S5]VC0Δt+N5(0,1)p4+δ4[S4]VC0Δtd[S5]dt=kf[A][S4]kb[S5](δ5+λ)[S5]+N3(0,1)kf[A][S4]VC0ΔtN4(0,1)kb[S5]VC0Δt+N6(0,1)δ5[S5]VC0Δt
Appendix 1—table 4
Model B parameter values.
ParameterValueParameterValue
p24.968896n3:33.189190
t1:22.569361δ30.246561
n1:20.952178p44.674459
t3:20.131057δ41.010978
n3:29.219961kf3.194116
δ20.521437kb4.634009
p30.113586δ51.165439
t3:32.183079C01

Additional models

Here we present additional models that were evolved but not discussed in the main text to provide more examples of size control mechanisms.

Model A3

Model A3 is similar to Model A1, but lacks the homodimerization interaction of I (see Equation 24). This specific model results in a weak adder/timer that displays a non-linear size control response curve. This is similar to Model A2 where we see a sizer/adder behavior in the low control volume regime and adder/timer behavior in the high control volume regime.

This model’s behavior is summarized in Appendix 1—figure 8. The parameter values of the model are shown in Appendix 1—table 5 along with the corresponding differential equations in Equation 27. This model was evolved with the Pareto fitness optimization framework to maximize NDiv and minimize CVBirth. The initial seed model topology for this evolutionary simulation was the quantity sensing oscillator shown in Appendix 1—figure 3A and was optimized during 3000 epochs.

(27) d[I]dt=max(p2(S/G2/M Switch/t1:2)n1:21+(S/G2/M Switch/t1:2)n1:2,b2)-kf[I][S3]+kb[S4]-(δ2+λ)[I]d[S3]dt=p31+(I/t2:3)n2:3-kf[I][S3]+kb[S4]-(δ3+λ)[S3]d[S4]dt=kf[I][S3]-kb[S4]-(δ4+λ)[S4]
Appendix 1—figure 8
Model A3’s behavior.

(A) Network topology of the evolved Model A3. S3 is as a size sensor and titrates I in a size-dependent manner. (B) Size distributions at birth (red), G1/S (orange), and division (purple). (C) Added volumes in G1 (red), S/G2/M (blue) and over the whole cycle (purple) as a function of the volume at the beginning of those phases. (D) Temporal dynamics of the model, colors correspond to the variables in A. (E) Response curve T(VC) as a function of control volume VC at the G1/S transition.

Appendix 1—table 5
Model A3 parameter values.
ParameterValueParameterValue
p25.606156t2:31.410420
t1:20.066136n2:33.434213
n1:28.177965kf1.930909
b20.911279kb3.871610
δ20.257920δ30.013294
p35.887584δ41.803926

Model A4

Model A4 is similar in essence to Model A1, albeit more unstable. In this model, S4 is the size sensor. Instead of using a PPI to titrate I in a size-dependent manner in G1, this model leverages a transcriptional repression to modulate the production of I directly rather than its effective degradation. The model’s behavior is summarized in Appendix 1—figure 9. The parameter values of the model are shown in Appendix 1—table 6 along with the corresponding differential equations in Equation 28. Notably, Appendix 1—figure 9E shows that the model’s response curve in the low control volume regime is ill-defined. In this model specifically, when the low volume regime is reached, the concentration of inhibitor I at G1/S increases due to quantity sensing [I]1/V. Then, because of the homodimerization of I into S3, we see the concentration [S3] also rise. We can see both of these curves spike up momentarily in the trajectories of Appendix 1—figure 9D around t136. The problem arises if this increase is too strong. Then, S3 activates the production of additional I, kick-starting a positive feedback loop that creates more and more inhibitor I, effectively interrupting the oscillator and making the period ill-defined. In practice, due to intrinsic noise at the G1/S transition and in the S/G2/M timer length, we’ve seen a cell lineage terminate prematurely before it can reach the maximum number of NDiv allowed in a simulation because of this problem. One could say this model is thus less fit than those presented in the main text, although it displays an added volume slope over its cycle of –0.41 and is close to a sizer.

This model was evolved using a Pareto fitness optimization that maximizes NDiv and minimizes CVBirth. The initial model topology was the quantity sensing oscillator of Appendix 1—figure 3A which was optimized over 4000 epochs.

(28) d[I]dt=max(p2max((S/G2/M Switch/t1:2)n1:21+(S/G2/M Switch/t1:2)n1:2,([S3]/t3:2)n3:21+([S3]/t3:2)n3:2),b2)11+(S4/t4:2)n4:2-kf[I]2+kb[S3]-(δ2+λ)[I]d[S3]dt=kf[I]2-kb[S3]-(δ3+λ)[S3]d[S4]dt=p41+(I/t2:4)n2:4-(δ4+λ)[S4]
Appendix 1—figure 9
Model A4’s behavior.

(A) Network topology of the evolved Model A4. S4 is the size sensor and represses the production of I in a size-dependent manner. (B) Size distributions at birth (red), G1/S (orange) and division (purple). (C) Added volumes in G1 (red), S/G2/M (blue) and over the whole cycle (purple) as a function of initial volume at the start of those phases. (D) Temporal dynamics of the model, colors correspond to the variables in A. (E) Response curve T(VC) as a function of control volume VC at the G1/S transition.

Appendix 1—table 6
Model A4 parameter values.
ParameterValueParameterValue
p25.6604334δ20.001059
t1:20.408208kf0.843408
n1:21.769233kb1.680321
t3:20.887392δ30.825150
n3:26.437683p44.674459
t4:20.197495t2:40.744119
n4:23.238633n2:43.451469
b20.800255δ41.929584

Model A5

Model A5 is another variation on Model A1 but here within the S. pombe cell cycle framework where I controls the timing of division directly and the Switch is turned on in G1 instead of S/G2/M. Here, S3 directly senses size and does so via the PPI linking I, S3, and S7. Indeed, since I is inversely proportional to the volume at G1/S due to quantity sensing and since S7 is solely produced via complex formation of I with S3, S7 is also inversely proportional to volume. Consequently, S3, which is almost solely produced via the dissociation of S7, becomes a direct sensor of the size of the cell.

Then, instead of using a PPI to titrate I in a size-dependent manner as is done in Models A1 and A2, Model A5 leverages a transcriptional repression mechanism to modulate the production of I directly instead of its degradation, precisely like in Model A4. The model’s behavior is summarized in Appendix 1—figure 10. The parameter values of the model are shown in Appendix 1—table 7 along with the corresponding differential equations in Equation 29.

This model was evolved using a Pareto fitness optimization that maximizes N Div and minimizes CVBirth. The seed model topology was the quantity sensing oscillator of Appendix 1—figure 3A which was optimized over 2500 epochs.

(29) d[I]dt=p2max((G1 Switch/t1:2)n1:21+(G1 Switch/t1:2)n1:2,b2)11+(S3/t3:2)n3:211+(S4/t4:2)n4:2-kf1[I][S4]+kb1[S6]-kf2[I]2+kb2[S5]-kf3[I][S3]+kb3[S7]-(δ2+λ)[I]d[S3]dt=b3-kf3[I][S3]+kb3[S7]-(δ3+λ)[S3]d[S4]dt=b4-kf1[I][S4]+kb1[S6]-(δ4+λ)[S4]d[S5]dt=kf2[I]2-kb2[S5]-(δ5+λ)[S5]d[S6]dt=kf1[I][S4]-kb1[S6]-(δ6+λ)[S6]d[S7]dt=kf3[I][S3]-kb3[S7]-(δ7+λ)[S7]
Appendix 1—figure 10
Model A5’s behavior.

(A) Network topology of the evolved Model A5. S3 is the size sensor and represses the production of I in a size-dependent manner. (B) Size distributions at birth (red), G1/S (orange) and division (purple). (C) Added volumes in G1 (red), S/G2/M (blue) and over the whole cycle (purple) as a function of initial volume at the start of those phases. (D) Temporal dynamics of the model, colors correspond to the variables in A. (E) Response curve T(VC) as a function of control volume VC at division.

Appendix 1—table 7
Model A5 parameter values.
ParameterValueParameterValue
p22.379635b42.624939
t1:20.142770δ41.499653
n1:20.383384δ50.006318
t3:20.714386δ60.983140
n3:28.642875δ70.481592
t4:21.754776kf10.226150
n4:27.783080kb13.793388
b20.195855kf21.608325
δ20.576282kb23.322565
b30.582449kf32.999118
δ30.318466kb30.906997

Model A6

Model A6 is yet another version of Model A1 evolved with slightly different fitness functions. This model was evolved using a Pareto fitness optimization that maximizes NDiv and minimizes the sum of squared residuals from a target volume at birth as described in Equation 15. For reference, the target volume at birth chosen for this simulation was Vt=10. The initial model topology was the quantity sensing oscillator of Appendix 1—figure 3A which was optimized over 3000 epochs.

The model’s behavior is summarized in Appendix 1—figure 11. The parameter values of the model are shown in Appendix 1—table 8 along with the corresponding differential equations in Equation 30.

(30) d[I]dt=p2(S/G2/M Switch/t1:2)n1:21+(S/G2/M Switch/t1:2)n1:2-kf1[I]2+kb1[S4]-kf2[I][S3]+kb2[S5]-(δ2+λ)[I]d[S3]dt=p311+(([I]/t2:3)n2:3)11+(([S3]/t3:3)n3:3)-kf2[I][S3]+kb2[S5]-(δ3+λ)[S3]d[S4]dt=kf1[I]2-kb1[S4]-(δ4+λ)[S4]d[S5]dt=kf2[I][S3]-kb2[S5]-(δ5+λ)[S5]
Appendix 1—figure 11
Model A6’s behavior.

(A) Network topology of the evolved Model A6. S3 is the size sensor and represses the production of I in a size-dependent manner. (B) Size distributions at birth (red), G1/S (orange) and division (purple). (C) Added volumes in G1 (red), S/G2/M (blue) and over the whole cycle (purple) as a function of initial volume at the start of those phases. (D) Temporal dynamics of the model, colors correspond to the variables in A. (E) Response curve T(VC) as a function of control volume VC at the G1/S transition.

Appendix 1—table 8
Model A6 parameter values.
ParameterValueParameterValue
p20.369408δ20.093159
t1:20.748731kf10.533415
n1:23.850889kb13.141514
p33.868044δ30.626507
t2:30.082081kf23.911233
n2:33.638939kb21.885120
t3:30.599375δ40.821609
n3:36.300643δ51.882484

Data availability

This is a theory paper, so there is no experimental data, and all results were generated by the code. The code used is freely available at https://github.com/FelixPG/PhiEvo_SizeControl, (copy archived at swh:1:rev:afa7f16a2f8a9d793aa3685116c2436faae100dd). Reference to the code has been added in the text.

References

    1. François P
    2. Hakim V
    (2005) Core genetic module: the mixed feedback loop
    Physical Review. E, Statistical, Nonlinear, and Soft Matter Physics 72:031908.
    https://doi.org/10.1103/PhysRevE.72.031908
  1. Book
    1. Lynch M
    (2007)
    The Origins of Genome Architecture
    Sinauer Associates.
  2. Book
    1. Morgan DO
    (2007)
    The Cell Cycle: Principles of Control
    OUP/New Science Press.

Decision letter

  1. Sandeep Krishna
    Reviewing Editor; National Centre for Biological Sciences­‐Tata Institute of Fundamental Research, India
  2. Aleksandra M Walczak
    Senior Editor; CNRS LPENS, France
  3. Kabir Husain
    Reviewer; University of Chicago, United States
  4. Shiladitya Banerjee
    Reviewer; Carnegie Mellon University, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Evolution of cell size control is canalized towards adders or sizers by cell cycle structure and selective pressures" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individuals involved in the review of your submission have agreed to reveal their identity: Kabir Husain (Reviewer #1); Shiladitya Banerjee (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

The reviewers liked many aspects of the manuscript but have suggested a number of revisions, which include showing the robustness of the results, a better justification of the assumptions/methods, a comparison with existing data and mathematical approaches, including clarification of how the evolutionary approach differs from other approaches, and a more careful rewording of the conclusions. Please address all of the reviewer comments (see below).

Reviewer #1 (Recommendations for the authors):

Figures and text:

I think the paper would be greatly strengthened by less dense Figures. As it stands, I found it difficult to figure out what the main points are. As an example: the main point of Figure 5, as I understand it, might be best served by a time-series plot of Slope \Δ-V_{cycle}, or some other summary statistic that shows the transient evolution of a sizer. Otherwise, this point is buried in the legend of panels B, E, and H.

On feedback control mechanisms:

There are no statistical summaries of the simulations described in Figures 2 and 3 of the main text -- Only Figure 4 (whose simulations are initialised with the evolved Model A1) contains statistics on evolved networks. If I understand the text, all the simulations resulted in topologically similar networks -- is this the case. What is the range of CV_births, and does the achieved CV_birth depend on the molecular implementation of the feedback control (PPI, dimerisation, transcriptional control), or do these molecular details affect the \Δ-V_{slope}?

On adders vs sizers (Figures4, 5, and lines 557 to 561 in the discussion):

Figure 4, 5, and the text around it, suggest that the CV_{birth} for adders is lower than that of sizers, in contrast to the common view that sizers are better at controlling cell size. I wonder if the distinction is between the *steady-state CV_{birth}* and the time taken to return to equilibrium from a large 'perturbation' of the cell size?

If it is true that adders are better at the former, but sizers are better at the latter, then would this also explain why cell size control late in the cycle tends to favour sizers? The intuition perhaps being that size control later in the cycle needs to deal with perturbations that occurred earlier in the cell cycle (as suggested by Lines 557 to 561)?

Reviewer #2 (Recommendations for the authors):

We have the following specific comments and recommendations for the authors.

1. Lines 74-77: This is the one piece of the introduction where readers unfamiliar with the eukaryotic cell cycle would be confused. Including background information about G1/S and S/G2/M phases would help expand the target audience of the paper since the techniques discussed therein are widely applicable.

2. Lines 101-106: The "poisson rate that corresponds to the deterministic rate" is unclear. These two sentences could be elaborated on further since significant prior knowledge of the reader is currently assumed.

3. Line 115: "While size-dependent growth mechanisms exist and do support size homeostasis" – This assertion should be backed up with relevant citations.

4. Lines 209-210: "Upon passing the G1/S transition, we assume cells are committed to division and there is a fixed time delay before they divide thus modeling S/G2/M as a timer." – What motivates this assumption? A citation or further discussion is warranted.

5. Line 285: This would be a good place for a citation to direct readers to sources discussing concentration vs quantity sensing. Note that in the bacterial size control literature, quantity sensing of division initiators has been shown to regulate adder behavior (Si et al., Curr Biol 2019).

6. Figure 1B: The dashed grey line is defined afterwards in 1C. It would be good to include it in the defined interactions here instead. In addition, the clarity of this schematic would be better with a line showing how the final network becomes the initial network in the next epoch.

7. Lines 301-302: "… its production is completely shut down in S/G2/M" does not come through clearly in the associated figure. You should clearly describe the chosen dynamics for the inhibitor protein in different phases of the cell cycle, and justify why the choices are different from known inhibitors such as Whi5.

8. Figure 2A: The message of this figure is not presented clearly. There is clustering with high CV volume and low N, another with negligible CV volume and widespread N division, and then the circled optimum. However, the trajectory of how a network evolves is not clear in this picture. Do all of them converge to the optimum eventually? Do they move to low CV before high N division or at the same time? How many epochs does it take to cross the large gap between the clustered networks and the optimum? Recommend somehow indicating sample evolutionary trajectories in addition to the aforementioned clarifications to remedy this issue. Additionally, why are there no numerical values in the axes? This makes it very difficult to assess the degree to which original values have changed.

9. Lines 309-322: This paragraph is somewhat confusing, in particular lines 314-316. The motivation for the control volume is unclear, especially in the physical sense of why a cell would use a non-physical volume to control a transition. While the idea makes sense later in the supplementary material, it needs to be clear from the very start that the goal here is that the control volume is a tool to examine how size at G1/S affects the cycle time.

10. Figures 3A,3C: While intuitive, specifying the role of the dashed red line would improve clarity.

11. Figure 3D: The predictions for the cell cycle scatter appear much stronger than the scatter itself. Can you comment on this?

12. Figure 3B, D: Compare how the model predictions compare with binned means for the scatter.

13. Lines 357-361: The numbers of 120 simulations and 500 epochs appear to be chosen arbitrarily. Why did you choose these initial conditions, and are the results of your paper robust with respect to higher/lower values? If so, including that point here would strengthen the argument, especially with a brief discussion on the lower limit. In addition, roughly how long in time is an epoch? The speed at which evolution is occurring would be of interest to many readers.

14. Figure 4: 4B is created to resemble S. pombe. Do the other panels have real-life analogies or are they arbitrarily chosen for qualitative representations of the discussed effects in the main text?

15. Figures 2F, 2I, 5A, 5D, 5G, 6C: The second zoomed-in panel of 6C is essential to understand the inner-generational dynamics of your modeling. The first many-generation panel shows stability in V but fails to address the other variables and the multi-phase dynamics. The other figures (2F, 2I, 5A, 5D, 5G) would benefit greatly from either a similar treatment or just fewer generations. The stability can be shown with significantly fewer divisions than are currently used.

16. Figure 5B: The ΔV scatters for S/G2/M and the whole cycle could be grouped into two – a positive correlation and a negative correlation. A best fit to the entire scatter is misleading therefore and does not describe the correlation trend.

17. Lines 539-542: Why is a one-step implementation of size control discarded? Surely a simpler control mechanism could be preferred naturally despite being a lesser theoretical interest in evolution simulations.

18. Lines 793-794: In this paper, you consider parallels to organisms that divide asymmetrically, such as budding yeast. Have you run simulations considering asymmetric division? Surely that would impact cell size distribution and variability.

19. Lines 794-795: Wouldn't disregarding one of the two daughter cells add a bias against faster dividing cells? I.e. if the number of divisions is one of your fitness functions, doesn't this method eliminate the natural advantage of a relatively larger population size for multi-cell level exponential growth? Also an issue at S130-132.

20. Lines 801-802: How is the extraction of the nullcline performed?

21. 829: Recommend providing the conversion from the arbitrary units used to physical values, here and all other figures.

22. 838-839: Model A2 does not receive an explanation comparable to A1 in this figure; either move to supplemental materials or explain it clearly as well.

23. Figures 6H, 6I, 6J: These subfigures are not very clear, together with captions for 6I and 6J that do not sufficiently explain to the reader how they are read. Why is the Burst Amplitude axis extended so far beyond the heatmap?

24. S120-S121: How do these theta and n theta values come to be? Currently, it seems like they are chosen with no supporting reasoning or explanation.

25. S239: V target missing a capital T.

26. S328-S329: Need to use a left apostrophe rather than two right apostrophes for epoch and generation.

27. Eq. (S18): Why do you minimize m+1 rather than just m?

28. S406-S408: Why do you choose these interactions to include? How much of all biochemical interactions do they encompass together? Are there others that you are aware of that you are choosing to neglect, and why? Can you provide citations and/or an argument to motivate this choice? This is another crucial ansatz for your modeling that needs to be discussed more carefully.

29. Supplementary Section 4: Translating the arbitrary units into physical values (when possible) would be immensely useful/helpful here.

30. Figure S6E: Why does cut off just below 1 (here, and not in any other plots)?

31. Since the model is generally applicable to any organism, comparisons to size control in bacterial cells (even qualitative) would be useful to widen the appeal. For example, could you predict why almost all bacterial cells (even evolutionary divergent ones) behave as adders? It has been shown that adder is regulated by threshold accumulation of an initiator protein that is produced at a rate proportional to cell volume, which your model could perhaps capture. Furthermore, many bacterial cells also exhibit biphasic size regulation during the cell cycle. It has been shown that Bacillus subtilis behave as sizers during the first phase, followed by a timer phase till division (DOI:10.1016/j.cub.2020.04.030). By contrast, Caulobacter crescentus cells implement a timer first, followed by an adder phase of size control (DOI:10.1038/nmicrobiol.2017.116). Both these organisms behave as approximate adders overall.

https://doi.org/10.7554/eLife.79919.sa1

Author response

Reviewer #1 (Recommendations for the authors):

Figures and text:

I think the paper would be greatly strengthened by less dense Figures. As it stands, I found it difficult to figure out what the main points are. As an example: the main point of Figure 5, as I understand it, might be best served by a time-series plot of Slope \Δ-V_{cycle}, or some other summary statistic that shows the transient evolution of a sizer. Otherwise, this point is buried in the legend of panels B, E, and H.

We have clarified the main text and the captions to better explain Figure 5 but otherwise have kept the figure mostly unchanged as we find useful to see the actual temporal dynamics of the networks throughout evolution to illustrate the sloppy sizer evolving into a weak adder in order to reduce the CV of the size distribution at birth. This key point has now been added to the discussion as described above. To re-emphasize its importance, and clarify the purpose of Figure 5, we have modified the top of the figure caption, which now includes the following sentence: “Evolutionary dynamics continually reduce the selected for CVextBirth and proceed through a noisy sizer to a less noisy adder”

On feedback control mechanisms:

There are no statistical summaries of the simulations described in Figures 2 and 3 of the main text -- Only Figure 4 (whose simulations are initialised with the evolved Model A1) contains statistics on evolved networks. If I understand the text, all the simulations resulted in topologically similar networks -- is this the case.

It is the case as all successful evolution simulations ended with some version of Model A. In Figure 2, we present two versions of Model A that have slightly different network topologies but perform the feedback mechanism in the same way. In the Supplement, we also give additional examples of evolved models in Figures S8, S9 and S10. We rewrote the second paragraph of the Results section, which now begins as: “Evolution simulations are in part reproducible and most often lead to similar network topologies. The evolution trajectory leading to Model A1 is a typical example in which all simulations produced similar networks (Figure 2B).”

What is the range of CV_births, and does the achieved CV_birth depend on the molecular implementation of the feedback control (PPI, dimerisation, transcriptional control), or do these molecular details affect the \Δ-V_{slope}?

We have now included a Table in the Supplement where we show the ranges of CV_births obtained by our evolved models (see new Table S1). We’ve seen different combinations of biochemical interactions evolve and perform feedback in different ways. Yet, they generally result in similar CV_Birth as seen in the 5 different versions of Model A shown in Figures 2, 3, S8, S9, and S10. Our understanding is that the resulting \Δ V_Slope is a direct consequence of the strength of the feedback mechanism which is in itself dictated by the molecular interactions of the network. But, the feedback control can be implemented in different ways with similar results. We now describe these results in the last paragraph of the section ‘Evolution of quantity-based size control mechanisms’ which reads as: “We give additional examples of similarly evolved networks in Figures S8-S10 where we can see the sensing and the feedback mechanism being implemented in different ways. Yet, despite these mechanistic differences in feedback regulation the resulting function of the evolved networks were similar as indicated by their CVBirth (Table S1).”

On adders vs sizers (Figures4, 5, and lines 557 to 561 in the discussion):

Figure 4, 5, and the text around it, suggest that the CV_{birth} for adders is lower than that of sizers, in contrast to the common view that sizers are better at controlling cell size. I wonder if the distinction is between the *steady-state CV_{birth}* and the time taken to return to equilibrium from a large 'perturbation' of the cell size?

If it is true that adders are better at the former, but sizers are better at the latter, then would this also explain why cell size control late in the cycle tends to favour sizers? The intuition perhaps being that size control later in the cycle needs to deal with perturbations that occurred earlier in the cell cycle (as suggested by Lines 557 to 561)?

As we discussed above, a noisy sizer can have a higher CV at steady state than an accurate adder. However, in line with the reviewers intuition, expect sizers to be better at reducing the effect of large fluctuations where the system is far from steady state. This is because a sizer will return the system to the steady state in one cell cycle. We now emphasize this point in the discussion, where we write:

‘However, we anticipate even noisy sizers will be better than adders at controlling cell size in response to large deviations away from the steady state distribution. This is because sizers will always return the cell size to be within the steady state distribution within a cell cycle.’

The argument here is agnostic to the distribution of cell size control in the different phases of the cell cycle and independent of the fact that G2 size control promotes sizers (which is discussed at length in the text already).

Reviewer #2 (Recommendations for the authors):

We have the following specific comments and recommendations for the authors.

1. Lines 74-77: This is the one piece of the introduction where readers unfamiliar with the eukaryotic cell cycle would be confused. Including background information about G1/S and S/G2/M phases would help expand the target audience of the paper since the techniques discussed therein are widely applicable.

We have added some basic information on the cell cycle and cited the main book of the field by David Morgan. The introduction now contains the following lines of text: “Cell size is regulated through the cell cycle control network that governs transitions from one phase of the cell cycle to the next. The division cycle can be broken up into distinct phases that are characterized by different molecular activities (D. O. Morgan, 2007). While it is typically considered that there are 4 phases of the cell cycle (G1, S, G2, and M), we here consider a two phase model based on a G1 phase and a composite S/G2/M phase. This is because size control in general has been associated with either the G1/S transition or mitosis at the end of the cell cycle.”

2. Lines 101-106: The "poisson rate that corresponds to the deterministic rate" is unclear. These two sentences could be elaborated on further since significant prior knowledge of the reader is currently assumed.

We now write: “Thus, each biochemical reaction takes place with a rate that corresponds to the deterministic rate, to which we add one white Gaussian noise with a variance equal to that rate. For example, given a deterministic biochemical rate k and a time interval of size extΔt, we consider a tau-leaping change of kΔt+N(0,kΔt) where N(0,kΔt) is a random gaussian variable of mean 0 and variance extkΔt.”

3. Line 115: "While size-dependent growth mechanisms exist and do support size homeostasis" – This assertion should be backed up with relevant citations.

We now cite Miettinen and Bjorklund (2016; https://doi.org/10.1016/j.devcel.2016.09.004) and Ginzburg et al. (2018; https://doi.org/10.7554/eLife.26957).

4. Lines 209-210: "Upon passing the G1/S transition, we assume cells are committed to division and there is a fixed time delay before they divide thus modeling S/G2/M as a timer." – What motivates this assumption? A citation or further discussion is warranted.

This assumption is justified by what happens in S. cerevisiae’s cell cycle. Upon passing the G1/S transition, during an event called Start, the cell becomes irreversibly committed to division. This means that at that point, even when treated with drugs that would disable the late cell cycle machinery, cells will divide no matter what. We now cite Doncic et al. (2011) as a reference for the commitment point, and Chandler-Brown et al. (2017) as a reference for the time nature of S/G2/M phase.

5. Line 285: This would be a good place for a citation to direct readers to sources discussing concentration vs quantity sensing. Note that in the bacterial size control literature, quantity sensing of division initiators has been shown to regulate adder behavior (Si et al., Curr Biol 2019).

As suggested, we now refer the reader to references discussing mechanisms to sense protein quantities rather than concentrations. We now discuss quantity sensing when it first is introduced in the methods section. The subsection on initial cell cycle model now contains the following sentences: “One way the amount rather than the concentration of a molecule could be sensed is through its titration against a fixed cellular quantity such as the genome, which is part of a general class of titration-based cell size sensing mechanisms (Amodeo et al., 2015; Heldt et al., 2018; Si et al., 2019; Wang et al., 2009).”

6. Figure 1B: The dashed grey line is defined afterwards in 1C. It would be good to include it in the defined interactions here instead. In addition, the clarity of this schematic would be better with a line showing how the final network becomes the initial network in the next epoch.

The dashed grey line is a special interaction that cannot be evolved or mutated by the PhiEvo algorithm. This grey dashed line always connects the inhibitor I to the S/G2/M Switch to indicate that I is in fact an inhibitor of the S/G2/M cell cycle phase. Because of this distinction, we want to separate this interaction from the evolvable interactions (transcriptional activation, transcriptional repression and complexation) that are shown in Figure 1B. To avoid confusion, we have removed the dashed grey line from the cartoon networks of Figure 1B as the focus of this panel is on the general evolution algorithm itself and not the specific implementation used in this study.

As for a line showing how the final network becomes the initial network in the next epoch, we have included a figure in the Supplement showing this process in more detail (Figure S4).

7. Lines 301-302: "… its production is completely shut down in S/G2/M" does not come through clearly in the associated figure. You should clearly describe the chosen dynamics for the inhibitor protein in different phases of the cell cycle, and justify why the choices are different from known inhibitors such as Whi5.

Here in these lines, we were talking about the R protein and not the I inhibitor that plays the role of Whi5 in our system. We apologize for this confusion and have clarified this point in the text, which now reads as: “For example, Model A2 contains extra interactions for the volume sensing gene [R], where [R] is repressed by the S/G2/M Switch (meaning its production is completely shut down in S/G2/M leading to sawtooth-like dynamics).”

8. Figure 2A: The message of this figure is not presented clearly. There is clustering with high CV volume and low N, another with negligible CV volume and widespread N division, and then the circled optimum. However, the trajectory of how a network evolves is not clear in this picture. Do all of them converge to the optimum eventually? Do they move to low CV before high N division or at the same time? How many epochs does it take to cross the large gap between the clustered networks and the optimum? Recommend somehow indicating sample evolutionary trajectories in addition to the aforementioned clarifications to remedy this issue. Additionally, why are there no numerical values in the axes? This makes it very difficult to assess the degree to which original values have changed.

We have completely remade this figure panel to render it more transparent (shown below) and have updated the caption accordingly. We have also added additional details about the evolutionary trajectory in the subsection S3A – Fitness of the Supplement.

Evolution happens in several stages. First, there are several epochs without any size control; networks then cluster in two regions of the Pareto front, essentially corresponding to volume going to 0 ([ii] in the new panel) and volume going to maximum volume ([i] in the new panel) both cases which are highly penalized in their fitness score (as is now explained in the Supplement). In Figure S3A, we show that the initial cell cycle network for the evolution process leads to unstable growth which is why we inevitably see the cell volume crash to 0 or maximum volume and cluster in [i] or [ii]. Evolution goes back and forth between those two clusters with a slow increase of the number of divisions (the fitness on the x axis). Eventually, some volume control evolves, most of the time corresponding to the Model 1 architecture described in the main text. Then, both the number of divisions N_div and the CV_Birth of those networks are optimized considerably, which is what we described in the main text as an “all or nothing” fitness score. Finally, CV_Birth keeps decreasing slowly until it reaches a plateau of 6-9% at the last generation. Not all networks go to optimum: in fact our Pareto evolution favours population diversity so makes sure that some networks are always relatively far from optimum. Also, notice that all simulations are different: those are stochastic simulations. Some simulations converge while others get lost on their way to the optimum as is expected. However, we implement the equivalent of a Drake’s rule: on average each network is mutated once per epoch. So the number of epochs before an evolutionary jump is a proxy of the typical number of mutations needed to select for one “good” mutation (changing Pareto front). We have added more details on all of this.

9. Lines 309-322: This paragraph is somewhat confusing, in particular lines 314-316. The motivation for the control volume is unclear, especially in the physical sense of why a cell would use a non-physical volume to control a transition. While the idea makes sense later in the supplementary material, it needs to be clear from the very start that the goal here is that the control volume is a tool to examine how size at G1/S affects the cycle time.

The reviewer is absolutely correct and phrases this statement better than we did. We have now redone Figure 3 rewritten this part of the text, which now reads as: “We then numerically integrate the differential equations of the model and measure the period T(VC) of the simulated cell-cycle for this control volume. Use of the control volume allows us to break the size feedback system and distinguish its input, VC, from its output, the induced cycle period T (Angeli et al., 2004).”

10. Figures 3A,3C: While intuitive, specifying the role of the dashed red line would improve clarity.

We clarified this and have updated the figure accordingly.

11. Figure 3D: The predictions for the cell cycle scatter appear much stronger than the scatter itself. Can you comment on this?

We think there is some confusion here. We note that we are predicting the average response using the control volume framework described above. We do not predict the scatter, which results from the stochastic simulations in steady state. We have clarified this in the caption of Figure 3B which now reads: “Added volumes extΔV for different phases of the cell cycles for simulations of Model A1. Individual dots correspond to different cell cycles for a simulation at steady-state. The full line corresponds to the extrapolation from the T(VC) curve shown in A for a restricted range of VC relevant to the scatter. The black cross, star and square indicate the average added volumes corresponding to when the system senses a volume corresponding to VG1/Sangle at the G1/S transition. We see that the model is predicted to follow an adder over a large range of volumes.”

12. Figure 3B, D: Compare how the model predictions compare with binned means for the scatter.

In evolutionary simulations, fitness typically improves rapidly in the beginning but then plateaus and only very gradually increases after a couple of hundred epochs. This is a general property of optimization simulations, as is also generally seen in machine learning. Based on our extensive previous experience with PhiEvo simulations, 500 epochs typically is sufficient to find the fitness plateau, but not so much that the network extensively explores the neutral mutations around the plateau. An epoch has no length per se, it just is a cycle of mutation/selection, however, the evolutionary algorithm adjust its mutation rates to have on average 1 mutation per generation (this is in fact an experimental fact called Drake’s law, and practically it prevents the so called “code-bloat” that can be observed in evolutionary simulations). This means that if a network evolves within, say, 100 generations, it is at most 100 mutations away from the initial state, and practically much less than that because most mutations are either neutral or deleterious (and in that case are not kept). To clarify this point in the text, we now write: “We chose to use 500 epochs in our simulations because in our previous experience this was sufficient for networks to evolve to be near the optimum, but not so much that they were forced to extensively explore the effects of neutral mutations near the optimum.”

13. Lines 357-361: The numbers of 120 simulations and 500 epochs appear to be chosen arbitrarily. Why did you choose these initial conditions, and are the results of your paper robust with respect to higher/lower values? If so, including that point here would strengthen the argument, especially with a brief discussion on the lower limit. In addition, roughly how long in time is an epoch? The speed at which evolution is occurring would be of interest to many readers.

For clarification, as suggested by the reviewer, we have redone the figure panels on the dynamics with fewer generations to allow for better visualization of the multi-generational protein and volume dynamics.

14. Figure 4: 4B is created to resemble S. pombe. Do the other panels have real-life analogies or are they arbitrarily chosen for qualitative representations of the discussed effects in the main text?

Panel A was designed to resemble S. cerevisiae, at least qualitatively, which was discussed in the text as having a size control at the G1/S transition.

15. Figures 2F, 2I, 5A, 5D, 5G, 6C: The second zoomed-in panel of 6C is essential to understand the inner-generational dynamics of your modeling. The first many-generation panel shows stability in V but fails to address the other variables and the multi-phase dynamics. The other figures (2F, 2I, 5A, 5D, 5G) would benefit greatly from either a similar treatment or just fewer generations. The stability can be shown with significantly fewer divisions than are currently used.

For clarification, as suggested by the reviewer, we have redone the figure panels on the dynamics with fewer generations to allow for better visualization of the multi-generational protein and volume dynamics.

16. Figure 5B: The ΔV scatters for S/G2/M and the whole cycle could be grouped into two – a positive correlation and a negative correlation. A best fit to the entire scatter is misleading therefore and does not describe the correlation trend.

This is a good point, we have done this as there are indeed two distinct behaviors as pointed out by the referee depending on whether or not the cell is small or large. The figure panel now looks as follows and we have adjusted the caption accordingly:

17. Lines 539-542: Why is a one-step implementation of size control discarded? Surely a simpler control mechanism could be preferred naturally despite being a lesser theoretical interest in evolution simulations.

As explained above in response to the reviewers major comment #4, we chose to do the analysis the way we did because we want to explore how cell size control can be done by a network with multiple feedbacks rather than just the concentration of a single protein, such as the budding yeast Whi5 protein, that has a special dedicated synthesis mechanism to make its concentration directly reflect cell size.

18. Lines 793-794: In this paper, you consider parallels to organisms that divide asymmetrically, such as budding yeast. Have you run simulations considering asymmetric division? Surely that would impact cell size distribution and variability.

It is definitively of interest to explore the effects of asymmetric divisions, but this is outside the scope of this already dense manuscript and will be the subject of future investigations.

19. Lines 794-795: Wouldn't disregarding one of the two daughter cells add a bias against faster dividing cells? I.e. if the number of divisions is one of your fitness functions, doesn't this method eliminate the natural advantage of a relatively larger population size for multi-cell level exponential growth? Also an issue at S130-132.

Here, cell growth is exponential on the single cell level since we were assessing size control mechanisms that take size as an input to cell cycle control. We are not exploring the very interesting case where growth deviates from the exponential. In that case, what the reviewer says is absolutely essential because then size homoeostasis would have a contribution from some cells in the population outcompeting others in terms of their growth. We intend to explore this possibility in future work. We have now added this text to the supplementary material which reads: “Here, cell growth is exponential on the single cell level since we were assessing size control mechanisms that take size as an input to cell cycle control. We are not exploring the very interesting case where growth deviates from the exponential. In that case, size homoeostasis would have a contribution from some cells in the population outcompeting others in terms of their growth and we would have to simulate the entire cell population and not disregard one of the daughter cells as we do here.”

20. Lines 801-802: How is the extraction of the nullcline performed?

We clarify in the text how the fictitious nullcline is extracted. We write:

“An intermediate fictitious nullcline is shown as a line that connects the average concentration [I] at G1/S and at division.”

21. 829: Recommend providing the conversion from the arbitrary units used to physical values, here and all other figures.

Our arbitrary units are described in the methods where it states that: “Note that we scale all our variables so that a concentration of one arbitrary unit corresponds roughly to 1000 proteins in a 100fL cell (Milo et al., 2010). Additionally, we scale the time variable so that 1 arbitrary time unit corresponds roughly to 30 min (Di Talia et al., 2007).” However, we prefer to have the figures in AU since this leads to a numerically simpler presentation and corresponds directly to what we have evolved.

22. 838-839: Model A2 does not receive an explanation comparable to A1 in this figure; either move to supplemental materials or explain it clearly as well.

We think the shorter explanation of A2 is fine for the figure caption because there is a longer explanation in the main text to explain this model. We now refer the reader to see the text in the figure caption.

23. Figures 6H, 6I, 6J: These subfigures are not very clear, together with captions for 6I and 6J that do not sufficiently explain to the reader how they are read. Why is the Burst Amplitude axis extended so far beyond the heatmap?

To make these figure panels more clear, we truncated the Burst Amplitude axis as suggested by the reviewer, and moved the inset to be its own panel. We also modified the caption to reflect these changes.

24. S120-S121: How do these theta and n theta values come to be? Currently, it seems like they are chosen with no supporting reasoning or explanation.

This is an important point. We chose these values because they give a similar amount of noise in the G1/S transition as observed experimentally (e.g., Di Talia et al. 2007; Chandler-Brown et al. 2017). We have added this explanation in the supporting information.

25. S239: V target missing a capital T.

This correction has been made.

26. S328-S329: Need to use a left apostrophe rather than two right apostrophes for epoch and generation.

This correction has been made.

27. Eq. (S18): Why do you minimize m+1 rather than just m?

In principle, it is the same. We chose to optimize m+1 rather than m in order to keep the fitness function positive. This helps when imposing multiplicative penalties on the fitness when cell volume is too low or too high and does not affect the optimization process.

28. S406-S408: Why do you choose these interactions to include? How much of all biochemical interactions do they encompass together? Are there others that you are aware of that you are choosing to neglect, and why? Can you provide citations and/or an argument to motivate this choice? This is another crucial ansatz for your modeling that needs to be discussed more carefully.

We chose these interactions because they are commonly used in the systems biology literature. Moreover, we have found in our previous evolutionary simulations that a combination of transcriptional interactions with protein-protein interactions is sufficient to account for many standard mechanisms in systems biology (e.g., switches, oscillators, biochemical adaptation). We note that other work in this area often restricts itself to purely transcriptional networks because they are easier to study and understand, and we added references to those works. Our approach is more flexible and generic.

We now write in the Supplement: “Many “inverse-approach” approaches in systems biology have focused on purely transcriptional networks (Francois et al., 2007, Fujimoto et al., 2008, Cotterell et al., 2010, Ten Tusscher et al., 2011), because they are generic, easier to study and can efficiently describe many biological dynamics (Alon, 2007).

In this project, we extend the biochemical interactions available for evolution: we not only model transcriptional activation and repression but also include complexation also known as protein-protein interaction (PPI), and assume there passive degradation. Adding PPIs is especially crucial because they are well known to lead to non-linear effects (Buchler et al., 2008, Buchler et al., 2009) allowing for the simple implementation of complex dynamics such as genetic oscillations (François et al., 2005) observed e.g. in circadian clocks (François, 2005), and such non-linear effects indeed play crucial roles for control in our evolved model. The equations for these interactions are presented in this subsection.”

29. Supplementary Section 4: Translating the arbitrary units into physical values (when possible) would be immensely useful/helpful here.

We have made the requested modification.

30. Figure S6E: Why does cut off just below 1 (here, and not in any other plots)?

For smaller Vc than the cutoff, Model A4 does not produce an oscillation, but instead reaches a steady state so there is no defined T/tau. We now note this fact in the new S9E (old S6E) caption and include a more detailed explanation for this model in the Supplement.

We write: “Notably, Figure S9E shows that the model's response curve in the low control volume regime is ill-defined. In this model specifically, when the low volume regime is reached, the concentration of inhibitor I at G1/S increases due to quantity sensing [I]  1/V. Then, because of the homodimerization of I into S3, we see the concentration [S3] also rise. We can see both of these curves spike up momentarily in the trajectories of Figure S9D around t  136. The problem arises if this increase is too strong. Then, S activates the production of additional I, kick-starting a positive feedback loop that creates more and more inhibitor I, effectively interrupting the oscillator and making the period ill-defined.”

31. Since the model is generally applicable to any organism, comparisons to size control in bacterial cells (even qualitative) would be useful to widen the appeal. For example, could you predict why almost all bacterial cells (even evolutionary divergent ones) behave as adders? It has been shown that adder is regulated by threshold accumulation of an initiator protein that is produced at a rate proportional to cell volume, which your model could perhaps capture. Furthermore, many bacterial cells also exhibit biphasic size regulation during the cell cycle. It has been shown that Bacillus subtilis behave as sizers during the first phase, followed by a timer phase till division (DOI:10.1016/j.cub.2020.04.030). By contrast, Caulobacter crescentus cells implement a timer first, followed by an adder phase of size control (DOI:10.1038/nmicrobiol.2017.116). Both these organisms behave as approximate adders overall.

We thank the reviewer for indicating those references, that we now cite in the manuscript. We have a generic argument for why adders might arise, ie, if size control takes place early in the division cycle, but do not have any specific things to say about bacteria in comparison to eukaryotic cells.

https://doi.org/10.7554/eLife.79919.sa2

Article and author information

Author details

  1. Felix Proulx-Giraldeau

    Department of Physics, McGill University, Montreal, Canada
    Contribution
    Conceptualization, Software, Validation, Investigation, Visualization, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0238-5410
  2. Jan M Skotheim

    1. Department of Biology, Stanford University, Stanford, United States
    2. Chan Zuckerberg Biohub, San Francisco, San Francisco, United States
    Contribution
    Conceptualization, Supervision, Methodology, Writing - original draft, Project administration, Writing - review and editing
    Competing interests
    No competing interests declared
  3. Paul François

    Department of Physics, McGill University, Montreal, Canada
    Contribution
    Conceptualization, Software, Supervision, Funding acquisition, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    paul.francois2@mcgill.ca
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2223-839X

Funding

Natural Sciences and Engineering Research Council of Canada (Discovery Grant)

  • Paul François

Natural Sciences and Engineering Research Council of Canada (Alexander Graham Bell Canada Graduate Scholarship)

  • Felix Proulx-Giraldeau

Fonds de recherche du Québec – Nature et technologies (Doctoral research scholarship)

  • Felix Proulx-Giraldeau

National Institutes of Health (NIH R35 GM134858)

  • Jan M Skotheim

Chan Zuckerberg Initiative (Biohub Investigator Award)

  • Jan M Skotheim

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Rodrigo Reyes-Lamothe, Nicolas E Buchler, and Lucas Fuentes Valenzuela for helpful comments on the manuscript. JS was supported by the NIH (R35 GM134858), PF was supported by Natural Sciences and Engineering Research Council of Canada (NSERC), Discovery Grant Program, FPG was supported by a Fonds de Recherche du Québec Nature et Technologies (FRQNT) Doctoral scholarship (B2X) and by a Natural Sciences and Engineering Research Council of Canada (NSERC) Doctoral Canada Graduate Scholarship (CGS-D).

Senior Editor

  1. Aleksandra M Walczak, CNRS LPENS, France

Reviewing Editor

  1. Sandeep Krishna, National Centre for Biological Sciences­‐Tata Institute of Fundamental Research, India

Reviewers

  1. Kabir Husain, University of Chicago, United States
  2. Shiladitya Banerjee, Carnegie Mellon University, United States

Publication history

  1. Preprint posted: April 13, 2022 (view preprint)
  2. Received: May 2, 2022
  3. Accepted: September 25, 2022
  4. Accepted Manuscript published: September 30, 2022 (version 1)
  5. Version of Record published: November 11, 2022 (version 2)

Copyright

© 2022, Proulx-Giraldeau et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 633
    Page views
  • 231
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Felix Proulx-Giraldeau
  2. Jan M Skotheim
  3. Paul François
(2022)
Evolution of cell size control is canalized towards adders or sizers by cell cycle structure and selective pressures
eLife 11:e79919.
https://doi.org/10.7554/eLife.79919

Further reading

    1. Computational and Systems Biology
    2. Neuroscience
    Sergio Oscar Verduzco-Flores, Erik De Schutter
    Research Article Updated

    How dynamic interactions between nervous system regions in mammals performs online motor control remains an unsolved problem. In this paper, we show that feedback control is a simple, yet powerful way to understand the neural dynamics of sensorimotor control. We make our case using a minimal model comprising spinal cord, sensory and motor cortex, coupled by long connections that are plastic. It succeeds in learning how to perform reaching movements of a planar arm with 6 muscles in several directions from scratch. The model satisfies biological plausibility constraints, like neural implementation, transmission delays, local synaptic learning and continuous online learning. Using differential Hebbian plasticity the model can go from motor babbling to reaching arbitrary targets in less than 10 min of in silico time. Moreover, independently of the learning mechanism, properly configured feedback control has many emergent properties: neural populations in motor cortex show directional tuning and oscillatory dynamics, the spinal cord creates convergent force fields that add linearly, and movements are ataxic (as in a motor system without a cerebellum).

    1. Computational and Systems Biology
    2. Immunology and Inflammation
    Mingyao Pan, Bo Li
    Short Report Updated

    T cells are potent at eliminating pathogens and playing a crucial role in the adaptive immune response. T cell receptor (TCR) convergence describes T cells that share identical TCRs with the same amino acid sequences but have different DNA sequences due to codon degeneracy. We conducted a systematic investigation of TCR convergence using single-cell immune profiling and bulk TCRβ-sequence (TCR-seq) data obtained from both mouse and human samples and uncovered a strong link between antigen-specificity and convergence. This association was stronger than T cell expansion, a putative indicator of antigen-specific T cells. By using flow-sorted tetramer+ single T cell data, we discovered that convergent T cells were enriched for a neoantigen-specific CD8+ effector phenotype in the tumor microenvironment. Moreover, TCR convergence demonstrated better prediction accuracy for immunotherapy response than the existing TCR repertoire indexes. In conclusion, convergent T cells are likely to be antigen-specific and might be a novel prognostic biomarker for anti-cancer immunotherapy.