PDE models of adder mechanisms in cellular proliferation
aa r X i v : . [ q - b i o . CB ] M a r PDE MODELS OF ADDER MECHANISMS IN CELLULAR PROLIFERATION
MINGTAO XIA ∗ , CHRIS D. GREENMAN † , AND
TOM CHOU ‡ Abstract.
Cell division is a process that involves many biochemical steps and complex biophysical mechanisms. To simplify theunderstanding of what triggers cell division, three basic models that subsume more microscopic cellular processes associated withcell division have been proposed. Cells can divide based on the time elapsed since their birth, their size, and/or the volume addedsince their birth – the timer, sizer, and adder models, respectively. Here, we propose unified adder-sizer models and investigatesome of the properties of different adder processes arising in cellular proliferation. Although the adder-sizer model provides a directway to model cell population structure, we illustrate how it is mathematically related to the well-known model in which cell divisiondepends on age and size. Existence and uniqueness of weak solutions to our 2+1-dimensional PDE model are proved, leading to theconvergence of the discretized numerical solutions and allowing us to numerically compute the dynamics of cell population densities.We then generalize our PDE model to incorporate recent experimental findings of a system exhibiting mother-daughter correlationsin cellular growth rates. Numerical experiments illustrating possible average cell volume blowup and the dynamical behavior ofcell populations with mother-daughter correlated growth rates are carried out. Finally, motivated by new experimental findings,we extend our adder model cases where the controlling variable is the added size between DNA replication initiation points in thecell cycle.
Key words.
PDE, structured populations, cell size control
AMS subject classifications.
1. Introduction.
How cells regulate and maintain their sizes, as well as sizes of their appendages is alongstanding research topic in cell biology. Besides growth of an individual cell, the size distributions within apopulation of cells are also a quantity of interest. When considering proliferating cell populations, individual cellgrowth is interrupted by cell division events that generate smaller daughter cells. The biological mechanisms thatcontrol when and how a cell divides are complex and involve many steps such as metabolism, gene expression,protein production, DNA replication, chromosomal separation (for eukaryotic cells), and fission or cell wallformation [27, 13, 4, 3, 6]. These processes are regulated and may involve intricate biochemical signaling.Despite the complexity of cell growth and the cell cycle, three simple hypotheses for the underlying mech-anisms of cell division have been proposed. Cell division can be governed by cell age a , cell volume x [26], oradded volume since birth y [29, 28]. The division mechanism employed by a type of cell may be interrogated bytracking the volumes x , added volumes y , and ages a during division events. Volume growth of an individual cellcan be straightforwardly measured and can be modeled by an effective empirical law such as ˙ x = g ( a, x, y, t ). Acommonly used approximation that is supported by observations is the exponential growth law g ( x ) = λx [24].To describe population-level distributions, PDE approaches have been developed. For example, the timermodel, in which the cell division rate depends only on age of the cell is described by the classic McKendrick equa-tion for n ( a, t ), the expected density of cells at age a and time t [19, 9]. The McKendrick “transport” equationfor the cell density takes the form ∂ t n ( a, t ) + ∂ a n ( a, t ) = − ( µ ( a ) + β ( a )) n ( a, t ), in which β ( a ) and µ ( a ) are age-dependent birth and death rates, respectively. The associated boundary condition n ( t,
0) = 2 R a β ( s ) n ( s, t )d s describes the birth of zero-age cells. Fully demographically stochastic versions of the timer model have alsobeen recently developed [12, 5, 11].The timer (or age-dependent) model does not explicitly track cell sizes, but PDE models incorporatingsizer mechanisms have been developed [22, 8, 23]. In these studies size-dependent birth rates β ( x ) are pertinent.Depending on the form of β ( x ), cells can diverge in size x in the absence of death [16]. Existence and uniquenessof weak solutions to timer and sizer models have been proved for certain boundary and initial conditions. Thesetypes of structured population equations can be partially solved using the method of characteristics but theboundary conditions can only be reduced to a Volterra-type integral equation [22, 5].Much like a general growth law g ( a, x, y, t ) that can depend on age, size, added size, and time, the threedistinct mechanisms of cell division need not be mutually exclusive. In this paper, we mainly focus, at thecell population level, on the cell division mechanism that incorporates the added volume, or the so-called the“adder.” This mechanism, in which the cell seems to use added size as the factor controlling its division, hasbeen indicated in many recent experimental studies. Specifically, apart from the sizer and the timer models,the adder mechanism has been recently shown to be consistent with E. coli division [27, 28, 29] and can bemotivated by an initiator accumulation mechanism distinct from those used to justify sizers or timers [28, 4]. ∗ Department of Mathematics, UCLA, Los Angeles, CA, USA † School of Computing Sciences, University of East Anglia, Norwich, UK ‡ Depts. of Computational Medicine and Mathematics, UCLA, Los Angeles, CA, USA1 e will introduce the PDE model that describes cell population structure under the adder mechanism,which we describe as the “adder-sizer” PDE model, and show its connection to the classical “timer-sizer” PDEmodel that involves cell age and size as controlling parameters. The proof of the existence and uniqueness of aweak solution to the proposed three-dimensional “adder-sizer” PDE turns out to be more complex than the prooffor the timer and/or sizer counterparts [22]. Our proof leads to the convergence of the numerical solutions to theadder-sizer PDE, allowing us to numerically evaluate the corresponding structured cell populations, facilitatingfurther analysis, exploration of possible “blowup” behavior, and generalizations of the model. Stochastic Monte-Carlo simulations of the corresponding stochastic process are also generated and compared with numerical resultsfor n ( x, y, t ) and division-event densities.Next, we propose an extension to the adder-sizer model that incorporates cellular growth rates that arecorrelated across successive generations. Changes in growth rates at the single-cell level have been exploredusing stochastic mapping methods [7, 18]. By numerically solving the PDE, we found out that the population-averaged growth rate are larger when correlations between mother and daughter cell growth rates are larger.Finally, we generalized the adder model to include a different two-phase PDE system which could describe thelatest “initiation adder” mechanism, which states that the added mechanism takes effect on the cell’s size atinitiation instead of division in [25]. In contrast to the single-PDE division adder model, a model describing theinitiation adder mechanism requires two coupled PDEs.To model cell size control, stochastic maps that relate daughter cell sizes to mother cell sizes have beendeveloped [17, 21]. These models describe how cell sizes evolve with generation and can interpolate among timer,sizer, and adder mechanisms. Kessler and Burov [17] assumed stochastic growth which lead to a stochastic mapwith multiplicative noise. They found that an adder mechanism can admit “blow-up” in which the expected cellsizes can increase without bound with increasing generation observed experimentally in filamentous bacter. Modi et al. [21] assume additive noise and do not find blow-up in an adder model. Stochastic maps of generationalcell size do not describe population-level distributions in size or age.
2. Adder-sizer PDE models.
Here, we introduce adder-sizer PDE models and generalize them to de-scribe recently observed characteristics of population-level bacterial cell division. An adder-sizer model is onethat incorporates a cell division rate β ( x, y, t ) and a single-cell growth rate g ( x, y, t ) that, instead of depend-ing on a cell’s age a , are functions of cell size x and a cell’s volume added since birth y . Such an adder-sizerPDE model can be developed by defining n ( x, y, t )d x d y as the mean number of cells with size in [ x, x + d x ]and added volume in [ y, y + d y ]. As cells have finite size and their added volume must be less than total size, n ( x ≤ , y, t ) = n ( x, y ≥ x, t ) = 0. A derivation similar to that given in [20] for the sizer model yields a transportequation of the form(2.1) ∂n ( x, y, t ) ∂t + ∂ [ g ( x, y, t ) n ( x, y, t )] ∂x + ∂ [ g ( x, y, t ) n ( x, y, t )] ∂y = − β ( x, y, t ) n ( x, y, t )for the adder-sizer PDE. Here, we have neglected the effects of death, which can be simply added to theright-hand-side of Eq. 2.1.To explicitly outline our general derivation, consider the total population flux into and out of the sizeand added size domain Ω shown in Fig. 1(a) and define ˜ β ( x ′ , y ′ , z ′ , t )d z ′ as the rate of fission of cells of size x ′ and added size y ′ to divide into two cells, one with size in [ z ′ , z ′ + d z ′ ] and the other with size within[ x ′ − z ′ , x ′ − ( z ′ + d z ′ )]. For binary fission, conservation of daughter cell volumes requires ˜ β ( x ′ , y ′ , z ′ , t ) ≡ ˜ β ( x ′ , y ′ , x ′ − z ′ , t ). This differential division function allows mother cells to divide into two daughter cells ofdiffering sizes (asymmetric division), a process that has been observed in numerous contexts [14, 13, 2]. We alsoassume that daughter cells must have positive size so ˜ β ( x ′ , y ′ , z ′ = 0 , t ) = ˜ β ( x ′ , y ′ , z ′ = x ′ , t ) = 0.The change in the number of cells in Ω due to fission can arise in a number of ways. First, if a cell in Ωdivides, it can only produce two cells with size less than x . Thus, such fission events lead to a net change of+1 in the number of cells with y = 0 and size in [0 , x ]. If a cell with size within [0 , x ] but with added size > y divides, it creates two cells with added size y = 0 and size within [0 , x ], leading to a net change of +2 cells.For cells with any added size y ′ > x ′ > x , we have two subcases. If the dividing cell has size x < x ′ < x , it will produce one daughter cell in Ω if a daughter cell has size 0 < z ′ < x ′ − x or x < z ′ < x ′ asshown in Fig. 1(b). If x ′ − x < z ′ < x , both daughter cells have size < x . Finally, if the dividing cell has size x ′ > x , at most one daughter will have size x ′ < x (see Fig. 1(b)). Upon simplifying the above birth termsby using R x ′ d z ′ = R x d z ′ + R x ′ x d z ′ for x ′ > x and the symmetry ˜ β ( x ′ , y ′ , z ′ , t ) = ˜ β ( x ′ , y ′ , x ′ − z ′ , t ), we combineterms to balance proliferation with transport and find ig. 1: The size and added-size state space for cell populations. The expected total number of cells at time t with added size within [0 , y ] and volume (or “size”) within [0 , x ] is defined as N ( x, y, t ). Over an increment intime d t , the domain Ω = [0 , y ] × [0 , x ] infinitesimally distorts Ω → Ω + dΩ through the growth increment g d t .The total population within this distorted domain changes only due to birth and death. Cells within Ω thatdivide always give rise to two daughters within Ω, leading to a net change of +1 cell. (b) The z ′ and x ′ domainsof the differential birth rate function ˜ β ( x ′ , y ′ , z ′ , t ). Cells outside of Ω can contribute a net +1 or +2 cells in Ωdepending on the division patterns defined in the depicted regions. Z x d x ′ Z y d y ′ ∂n ( x ′ , y ′ , t ) ∂t + Z x d x ′ g ( x ′ , y, t ) n ( x ′ , y, t ) + Z y d y ′ g ( x, y ′ , t ) n ( x, y ′ , t )= Z ∞ d y ′ Z x d x ′ Z x ′ d z ′ ˜ β ( x ′ , y ′ z ′ , t ) n ( x ′ , y ′ t )+ Z ∞ y d y ′ Z x d x ′ Z x ′ d z ′ ˜ β ( x ′ , y, z ′ , t ) n ( x ′ , y ′ , t )+ 2 Z ∞ d y ′ Z ∞ x d x ′ Z x d z ′ ˜ β ( x ′ , y ′ , z ′ , t ) n ( x ′ , y ′ , t ) . (2.2)Upon taking the derivatives ∂ ∂x∂y , we find the PDE given in Eq. 2.1 where the total division rate is defined by β ( x, y, t ) := R x ˜ β ( x, y, z, t )d z . For the boundary condition at y = 0, we take the derivative ∂/∂x and set y → + to find(2.3) g ( x, y = 0 , t ) n ( x, y = 0 , t ) = 2 Z ∞ x d x ′ Z x ′ d y ′ ˜ β ( x ′ , y ′ , z = x, t ) n ( x ′ , y ′ , t ) . The other boundary condition defined by construction is n ( x, x, t ) = 0.In the special restricted case of symmetric cell division, ˜ β ( x, y, z, t ) = β ( x, y, t ) δ ( z − x/ g ( x, y = 0 , t ) n ( x, y = 0 , t ) = 4 Z x β (2 x, y ′ , t ) n (2 x, y ′ , t )d y ′ . The above derivation provides an explicit boundary condition representing newly born cells that may be asym-metric in birth size. Quantities such as the total cell population N ( t ) and the mean total biomass M ( t ) (thetotal volume over all cells) can be easily constructed from the density n ( x, y, t ):(2.5) N ( t ) = Z ∞ d x Z x d y n ( x, y, t ) , M ( t ) = Z ∞ d x Z x d y xn ( x, y, t ) . Higher moments of the total volume can also be analogously defined. By applying these operations to Eq. 2.1and using the boundary condition (Eq. 2.3), we find the dynamics of the total population and biomass N ( t )d t = Z ∞ d x Z x d y β ( x, y, t ) n ( x, y, t ) , d M ( t )d t = Z ∞ d x Z x d y g ( x, y, t ) n ( x, y, t ) . Finally, we also define the distribution of division events over the size and added size variables, accumulatedover a time T :(2.7) ρ d ( x, y, T ) = Z T β ( x, y, t ) n ( x, y, t )d t Z T d t Z ∞ d x ′ Z x d y ′ β ( x ′ , y ′ , t ) n ( x ′ , y ′ , t ) . In general, the birth rate functions˜ β ( x, y, z, t ) and β ( x, y, t ) associated with adder-sizer models can take many forms that make biological sense.However, some classes of β ( x, y, t ) may allow the adder-sizer model to be transformed into the well-known “sizer-timer” structured population model [26]. To illustrate the relationship, we consider a division rate function β which depends explicitly only on age a and see how it could be converted to a function of size and added size.For a cell born at time t , the probability that the cell splits within time [ a, a + d a ] is defined by γ ( a ; ¯ a )d a .In the absence of death, to ensure that any single cell will eventually split, R ∞ γ ( a ; ¯ a )d a = 1. Reasonablechoices for γ ( a ; ¯ a ) are Gamma, lognormal, or normal distributions. Without loss of generality, we propose asimple gamma distribution for γ ( a ; ¯ a ):(2.8) γ ( a ; ¯ a ) = 1 a Γ((¯ a/σ a ) ) exp " − a ¯ aσ a + (cid:18) ¯ aσ a (cid:19) ln (cid:18) a ¯ aσ a (cid:19) , where ¯ a is the mean division age and σ a is the variance. This type of distribution can be derived from the sumof independent, exponentially distributed ages.For determinisitic exponential growth g = λx , age a and the parameter ¯ a can be explicitly expressed interms of x, y and possibly other fixed parameters:(2.9) a ( x, y ) = 1 λ ln (cid:18) xx − y (cid:19) , ¯ a ( x, y ) = 1 λ ln (cid:18) x − y + ∆ x − y (cid:19) , in which ∆ is the fixed added size parameter that represents the adder mechanism.With a ( x, y ) and ¯ a ( x, y ) defined in Eqs. 2.9, the division rate function β ( x, y ) can be expressed in terms of x and y by using the splitting probability γ ( a ( x, y ); ¯ a ( x, y )):(2.10) β ( x, y, t ) = γ ( a ( x, y ); ¯ a ( x, y ))1 − R a ( x,y )0 d a ′ γ ( a ′ ; ¯ a ( x, y )) . Assuming this “hazard function” form of a growth law, cells born at small initial size x (0) = x = x − y take longer time to divide, while cells born with large size split sooner. Using the gamma distribution, we finda division rate of the form(2.11) β ( x, y ) = Γ (cid:16) ¯ a ( x,y ) σ a (cid:17) γ ( a ( x, y ); ¯ a ( x, y ))Γ (cid:16) ¯ a ( x,y ) σ a , a ( x,y )¯ a ( x,y ) σ a (cid:17) , where Γ( · , · ) is the upper incomplete gamma function. We plot two examples of the time-independent rate β ( x, y ) in Fig. 2.With β ( x, y, t ) defined, we still need to construct the full fission rate ˜ β , which we will assume is a productof the overall division rate β ( x, y, t ) and a differential division probability. The simplest model is to assume ig. 2: The size and added-size dependent rate β ( x, y ) constructed using a gamma distribution for the splittingprobability γ (Eq. 2.8) and Eq. 2.10. We show projections at fixed values of x . In (a) the parameters are σ a = 0 .
2, while in (b) σ a = 1. Note the difference in scale and that γ ( a ) with a higher standard deviationleads to a lower overall cell division rate β . When x is large, ¯ a defined in 2.9 is small, a nonzero division rate β ( x, y → > β as a hazard function. Modifying birth rate at smallvalues of y so that β ( x, y = 0) → y .that the differential division probability h ( r ) is a function of only the ratio r between the size of the daughtercell and that of the mother cell, and independent of the cell size just before division. Thus,(2.12) ˜ β ( x, y, z, t ) = β ( x, y, t ) h ( z/x ) /x, where r ≡ z/x ∈ [0 , g ( x, , t ) n ( x, , t ) = 2 Z ∞ x d x ′ Z d s β ( x ′ , sx ′ , t ) h ( x/x ′ ) n ( x ′ , sx ′ , t ) . A reasonable model for h ( r = x/x ′ ) is a lognormal form that is symmetric about r = 1 / h ( r ) = h ( r ) + h (1 − r ) Z ( σ r , δ ) , h ( r ) = e − ( − δ +ln r )22 σ r e − ln2(1 − r )2 σ r , where the parameters δ and σ r determine the bias and spread of the daughter cell size distribution, and thenormalization constant is Z ( σ r , δ ) = R ( h ( r ) + h (1 − r ))d r . With the differential birth ratefunction ˜ β defined, we can now consider the implementation of numerical solutions to Eqs. 2.1 and 2.3 as wellas event-based simulations of the underlying corresponding stochastic process. Since a typical initial conditionmay not be smooth, a classical solution to Eqs. 2.1 and 2.3 may not exist. Thus, we provide a proof of existenceand uniqueness of the weak solution to Eqs. 2.1 and 2.3 in Appendix A. We show convergence of a discreteapproximation to our problem, allowing us to confidently numerically approximate the weak solution.The numerical approximation to the weak solution will be based on an upwind finite difference scheme inwhich both x and y are discretized with step size h . We define locally averaged functions by(2.15) f i + ,j + := 1 h Z ( i +1) hih d x Z ( j +1) hjh d y f ( x, y, t ) , here f ( x, y, t ) can represent n ( x, y, t ), g ( x, y, t ), or β ( x, y, t ). Similarly,(2.16) ˜ β i + ,j + (( s + 12 ) h, t ) = h − Z ( i +1) hih d x Z ( j +1) hjh d y Z ( k +1) hkh d z ˜ β ( x, y, z, t )in the domain i, j ≥ j, k < i . The discretization of the transport equation can be expressed as(2.17) n i + 12 ,j + 12 ( t +∆ t ) − n i + 12 ,j + 12 ( t )∆ t + g i +1 ,j + 12 ˜ n i +1 ,j + 12 − g i,j ˜ n i,j + 12 h + g i + 12 ,j +1 ˜ n i + 12 ,j +1 − g i + 12 ,j ˜ n i + 12 ,j h = − β i + ,j + n i + ,j + ( t ) , for 1 ≤ i, j ≤ L , where Lh is the maximum size which we take sufficiently large such that n i,j>K ( t = 0) =0 , n i ≤ j = 0. We also set g i + ,i = 0 to prevent density flux out of the y < x domain. In Eq. 2.17, g i +1 ,j + ( t )can be taken as g (( i + 1) h, ( j + ) h, t ) while ˜ n i +1 ,j + ( t ) = R ( j +1) hjh d y n (( i + ) h, y, t ) is a finite-volume numericalapproximation to R ( j +1) hjh d y n (( i + 1) h, y, t ). The discretized version of the boundary condition (Eq. 2.3) can beexpressed as(2.18) g i + , n i + , ( t ) = 2 h L X k = i +1 k − X j =0 ˜ β k + ,j + (( i + 12 ) h, t ) n k + ,j + ( t ) . The full explicit discretization scheme for the numerical calculation is provided in Appendix B.Direct Monte-Carlo simulations of the birth process are also performed and compared with our numericallycomputed deterministic distributions (see Appendix C). We construct a list of cells and their associated sizesand their sizes at birth. This list is updated at every time step ∆ t . The cell sizes grow according to g ( x, y, t ).If a cell divides, the initial sizes of the daughter cells are randomly chosen according to the distribution h ( z/x ).The daughter cells then replace the mother cell in the list. Simulations of the underlying stochastic processresults in, at any given time, a collection of cells, each with a specific size and added size. This collection of cellsrepresents a realization of the population that should be approximated by the distributions that are solutionsto Eqs. 2.1 and 2.3.
3. Analysis and Extensions.
In this section, we numerically investigate the adder-sizer model and plotvarious cell population densities and birth event distributions under different parameter regimes. We also showthe consistency of numerical solutions of the adder-sizer PDE with results from direct Monte-Carlo simulations ofthe corresponding stochastic process, which demonstrates that numerical solutions of the linear PDE model forcell population is in agreement with single-cell level stochastic models. After investigating birth rate parametersthat can lead to blow-up of population-averaged cell sizes, we extend the basic adder model to include mother-daughter growth rate correlations and processes that measure added size from different points in the cell cycle, i.e. , an initation-adder model.
We evaluated our adder-sizer PDE model by using the divisionrate given in Eq. 2.10 and first assuming the simple and well-accepted growth function g ( x, y, t ) = λx . Fig. 3shows the numerical results for the density ¯ n ( x, y, t ) = n ( x, y, t ) /N ( t ) at successive times t = 1 , ,
12, respectively.Stochastic simulations of the underlying process yield cells populations consistent with the deterministic densitiesderived from the PDE model. In Fig. 4, we compare the cell densities ¯ n ( x, y, t ) the division event densities ρ d ( x, y, T ) for two different differential division functions h ( r ). As before, the more asymmetric the division thebroader the cell and event densities. At the single-cell level, a stochastic map model by Kessler and Burov as-sumed a multiplicative noise and predicted that cell sizes can eventually grow without bound, in agreement withwhat was experimentally observed for filamentous bacteria [17]. However, stochastic maps of generational cellsize do not capture population-level distributions in size or age. In this subsection, we will numerically explorehow a possible ”blowup” in the population-averaged cell volumes. Within PDE models that describe populationdistributions, timer and sizer mechanisms have been shown to exhibit blow-up depending on properties of thebirth rate β ( a, x ) [1, 7, 16]. Analysis of the conditions on full differential division rate ˜ β ( x, y, z, t ) that would - -- - - Fig. 3: Numerically computed densities ¯ n ( x, y, t ) = n ( x, y, t ) /N ( t ) using g ( x, y, t ) = λx and ˜ β ( x, y, z, t ) definedby Eqs. 2.10, 2.8, and 2.14. For all plots, we use σ a = 0 . γ ( a ) (Eq. 2.8) and rescale size in units of ∆. In(a-c), we use the sharp, single-peaked differential division function h ( r ) shown in the inset ( σ r = 0 . , δ = 0)and plot ¯ n ( x, y, , ¯ n ( x, y, n ( x, y, h ( r ) with parameters σ r = 0 . , δ = 0 .
7. In all calculations, weassumed an initial condition corresponding to a single newly born ( y = 0) cell with size x = 1. For moreasymmetric cell division in (d-f), the density spreads faster. In these cases, the densities closely approach asteady-state distribution by about t = 12. Also shown in each plot are realizations of Monte-Carlo simulationsof the discrete process. Individual cells are represented by blue dots which accurately sample the normalizedcontinuous densities ¯ n ( x, y, t ).result in blow-up in the “adder-sizer” PDE model is more involved. Here, we provide only a heuristic argumentfor sufficient conditions for blow-up.First, we characterize the shape of the densities in the adder-sizer model. In the analogous McKendrickequation [15] one can investigate the age profile defined by dividing the number density by the total populationsize. The long term age profile may be stable even when the total population size continuously increases. Wetake a similar approach here by analyzing ¯ n ( x, y, t ) = n ( x, y, t ) /N ( t ) where N ( t ) is given by Eq. 2.5. Writingthe adder-sizer PDE in terms of ¯ n , we find(3.1) ∂ ¯ n∂t + ¯ nN d N d t + ∂ ( g ¯ n ) ∂x + ∂ ( g ¯ n ) ∂y = − β ¯ n. Integrating this equation over x, y leads to ˙
N/N = R ∞ d x R x d y β ¯ n , which can be substituted into the first termin Eq. 3.1 to yield the nonlinear PDE(3.2) ∂ ¯ n∂t + ∂ ( g ¯ n ) ∂x + ∂ ( g ¯ n ) ∂y = − (cid:18) β + Z Ω β ¯ n (cid:19) ¯ n. A number of standard approaches may be applied to analyze Eq. 3.2. For example, in [15], solutions areattempted by controlling the analogous non-linear integral term. In the adder-sizer problem, we can define ig. 4: Comparison of cell densities ¯ n ( x, y, t ) and cell division event densities ρ d ( x, y, T ) (Eq. 2.7). The standarddeviation σ a = 0 . n ( x, y, t = 12) and ρ d ( x, y, T ) using σ r = 0 . , δ = 0 while in (c) an (d) we used a broader differential division function in which σ r = 0 . , δ = 0 . T = 12. h β ( t ) i = R Ω β ¯ n in the above expression to find a self-consistent condition on h β ( t ) i . One can also assess thesteady-state ¯ n ss by setting ∂ ¯ n ss ∂t = 0 and establishing convergence.One indication of blow-up is a diverging mean cell size h x ( t ) i = M ( t ) /N ( t ). By multiplying the Eq. 3.1 by x and integrating (using the boundary condition and symmetry of the ˜ β distribution) we find(3.3) d h x ( t ) i d t + h β ( t ) ih x ( t ) i = q ( t ) , in which q ( t ) := R Ω g ¯ n . If β , g , and ¯ n = ¯ n ss are time-independent and a steady state mean cell size exists, weexpect it to obey h x ( ∞ ) i = q ( ∞ ) / h β ( ∞ ) i . For the special case of deterministic exponential growth g ( x ) = λx ,we can write the time evolution of the mean size as(3.4) d h x ( t ) i d t = [ λ − h β ( t ) i ] h x ( t ) i , h β ( t ) i ≡ Z ∞ d x Z x d y β ( x, y, t )¯ n ( x, y, t ) . If β ( ∞ ) is bounded above by λ , then we expect blow-up. For β ( ∞ ) that is not bounded, as in our example(Eq. 2.10), one cannot determine if blow-up occurs without a more detailed and difficult analysis. Since theprecise conditions on β leading to cell volume explosion are difficult to find, we will explore this possiblephenomena using numerical experiments. We numerically examine the density n ( x, y, t → ∞ ) and the meancell size h x ( t ) i using the β, ˜ β defined in 2.10, 2.8, and 2.14.In Fig. 5(a) and (b) we plot the marginal distribution ¯ n ( x, t ) := R ∞ x d y n ( x, y, t ) / R ∞ d x R ∞ x d y n ( x, y, t ) fordifferent values of the division rate variability σ a at different times. The associated division rates correspond tothose plotted in Fig. 2(a) and (b). In Fig. 5(c) we plot the mean cell sizes h x ( t ) i = M ( t ) /N ( t ) corresponding tothe distributions in (a) and (b). For sufficiently broad division probabilities γ ( a ) (large σ a ), the division rates β are small, and h x ( t ) i fails to saturate and diverges. Fig. 5: (a) Size distributions ¯ n ( x, t ) for σ a = 0 . t = 1 , , ,
10. (b) ¯ n ( x, t = 1 , , ,
10) for σ a = 1, σ r = 0 .
1, and δ = 0. (c) The corresponding mean cell sizes h x ( t ) i . The curve associated with the σ a = 0 . σ a = 1 exhibits blow-up. However, the blowup is suppressed if a deathterm ( µ = ln 2) is included. Recent experiments indicate that the growth rate of amother cell is “remembered” by its daughter cells. For growth rates of the form g ( x, y, t ) = λx , the exponentialgrowth parameter λ between successive generations i, i + 1 have been proposed to evolve [18, 7]. In [18],fluctuations in λ have been discussed at the single-cell level to explore their effects on the population-averagedgrowth rate while in [7], changes in growth rates across two consecutive generations are modeled as a Markovprocess in order to estimate a division rate function β . In this subsection, we first introduce a generalizedadder-sizer PDE incorporating variability in λ and then explore the mother-daughter growth rate correlationaffects the population dynamics.A mother-daughter growth rate correlation between two consecutive generations can be described by(3.5) λ i +1 = ( λ i − ¯ λ ) R + ¯ λ + ξ, where ξ is a random variable, 0 ≤ R < λ is the meanlong-term, or preferred growth rate. Given a growth rate λ i of a mother cell, Eq. 3.5 describes the predictedgrowth rate λ i +1 of its daughter cells. We assume that the random variable has mean zero and is distributedaccording to some probability density P ( ξ ), which vanishes for ξ ≤ (1 − R )¯ λ to ensure that the growth ratesremain positive.To incorporate the memory of growth rates between successive generations in the adder-sizer PDE model,we extend the cell density in the growth rate variable λ . Thus, n ( x, y, t, λ ) is the density of cells with volume x ,added volume y , and growth rate λ . The growth function g ( x, y, t, λ ) is now explicitly a function of the growthrate λ . We propose the extended PDE model(3.6) ∂n ( x, y, t, λ ) ∂t + ∂ ( gn ) ∂x + ∂ ( gn ) ∂y = − β ( x, y, t ) n ( x, y, t, λ ) ,g ( x, , t, λ ) n ( x, , t, λ ) =2 Z ∞ d λ ′ Z ∞ x d x ′ Z x ′ d y ˜ β ( x ′ , y, x, t ) n ( x ′ , y, t, λ ′ ) P ( ξ = λ − Rλ ′ − (1 − R )¯ λ ) , ˜ β ( x, y, x ′ , t ) = ˜ β ( x, y, x − x ′ , t ) ,n ( x, y, , λ ) = n ( x, y, λ ) , A possible symmetric mean zero distribution that vanishes at − (1 − R )¯ λ takes on a log-normal form:(3.7) P ( ξ ) ∝ exp " − ln ( ξ + (1 − R )¯ λ )2 σ ξ − ln ((1 − R )¯ λ − ξ )2 σ ξ . If we start with one newly born daughter cell at size x and growth rate λ , the initial condition in our PDEmodel would be n ( x, y, λ ) = δ ( x − x ) δ ( y ) δ ( λ − λ ). .4 0.5 0.6 Fig. 6: Population-level evolution of cellular growth rate. Parameters used are ¯ λ = ln 2 , σ a = 0 . , σ r = 0 . , δ = 0.(a-b) The marginalized density ¯ n ( λ, t ) as a function of growth rate λ for no correlation ( R = 0) and initial growthrate λ = 0 .
55. The peak in the distribution broadens as the mean evolves towards the preferred mean value¯ λ = ln 2. (c) The evolution of the mean h λ ( t ) i for different values of correlation R . Note that the steady-statevalues h λ ( ∞ ) i depend on the correlation R .Numerical solutions of Eqs. 3.6 shown in Fig. 6 indicate that although ¯ λ is the same for two different cases, R = 0 and R = 0 .
4, their corresponding mean growth rates h λ ( t ) i converge to different values. For largercorrelation R , the daughter cells’ growth rates do not deviate much from those of their mothers’ growth rates.This means that the offspring of faster growing cells tend to grow faster and the offspring of slower growing cellstend to grow slower. Because it takes shorter time for faster cells to divide, they will produce more generationsof faster-growing cells, leading to a larger average growth rate defined as(3.8) h λ ( t ) i = R ∞ d x R x d y R ∞ d λ λn ( t, x, y, λ ) R ∞ d x R x d y R ∞ d λn ( t, x, y, λ ) . On the other hand, for a fixed mother growth rate λ i , smaller correlations R lead to mean daughter cellgrowth rates h λ i +1 i that are closer to ¯ λ . Since cells with growth rates less than ¯ λ will live longer before division,these cells persist in the population longer than those with larger λ , pushing the average growth rate h λ ( t ) i tovalues smaller than ¯ λ . Fig. 6(c) explicitly shows that when R = 0, the mean growth rate approaches a valuesmaller than ¯ λ = ln 2. Recent experiments suggest a new type of adder mechanism for bacterialcell size control [25]. Rather than a fixed volume added between birth and division as the primary controlparameter, new experimental evidence suggests that the control parameter in
E. coli is the added volumebetween successive initiations of DNA replication. Initiation occurs when the ori sites in a cell’s genome areseparated, leading to DNA replication and segregation. The number of ori sites depend on cell type and species,typically one in prokaryotic cells and more than one in eukaryotic cells. The initiation-adder model assumesthat a cell’s volume per initiation site (the ori site in the genome) tends to add a fixed volume between twoconsecutive initiations.If the number of ori sites in a cell is q , initiation increases the number to 2 q . Immediately after divisionand DNA separation, the number of ori s decreases back to q in each daughter cell.In this subsection, we generalize the adder PDE model to describe this new initiation-adder mechanism.We classify all cells into two subpopulations: cells that have not yet undergone initiation and cells that haveinitiated DNA replication but that have not yet divided. We define n ( x, y, t )d x d y as the expected number ofpre-initiation cells in with volume in [ x, x + d x ] and with added volume y < x in [ y, y + d y ]. Mean post-initiationcell numbers with volume in [ x, x + d x ] and added volume in [ y, y + d y ] are described by n ( x, y, t )d x d y . In thegeneral initiation-adder process, when a pre-initiation cell commences DNA replication (initiates) can dependon the volume or added volume. Thus, we describe transitions from a pre-initiation cell transitions into apost-initiation cell by the rate k i ( x, y, t ). After initiation, the number of ori sites doubles and the added volumeis reset to zero in the newly formed post-initiation cell. In analogy with the differential division rate in Eq. 2.1, yy n n y y Fig. 7: Schematic for the initiation adder process. DNA replication is initiated (indicated by the red dot) beforecopied DNA is segregated and cell division. In this example, q = 1 and y is and added volume per originationsite for two origination sites. The density of cells with q = 1 copy of DNA (before DNA replication initiation)is denoted n ( x, y, t ) while the density of cells post-initiation is denoted n ( x, y, t ), where y denotes the volumeadded after initiation. The factor that controls y + y in the initiation-adder model is the volume ∆ addedbetween successive initiation events, rather than between successive cell divisions. Thus, the controlled variable(added volume in this case) spans the pre-initiation and post-initiation states.we define β ( x, y, t ) as the rate of division of post-initiation cells. Under a general asymmetric division event, weassume that the added volume is divided proportionally to the volume of the daughter cells, i.e. , if the mothercell’s volume is x with added volume y since initiation, and if one daughter cell’s volume is z < x and theother daughter cell’s volume is x − z , the added volume since division for the first daughter will be set to yz/x while the added volume for the second daughter will be y ( x − z ) /x . The resulting PDE model now involves twocoupled densities n and n :(3.9) ∂n ( x, y, t ) ∂t + ∂ [ g n ] ∂x + ∂ [ g n ] ∂y = − k i ( x, y, t ) n + 2 Z ∞ x zx n ( z, yz/x, t ) ˜ β ( z, x, yz/x, t )d z,∂n ( x, y, t ) ∂t + ∂ [ g n ] ∂x + ∂ [ g n ] ∂y = − β ( x, y, t ) n ,n ( x, , t ) = 0 , g n ( x, , t ) = Z x k i ( x, y, t ) n ( x, y, t )d y,β ( x, y, t ) = Z x ˜ β ( x, z, y, t )d z, in which we have allowed for different growth rates in the different cell phases. Both n and n are defined in thedomain { R +2 ∩ { y < x }} × R + . These coupled PDEs are different from the PDE associated with the standard“division-adder” described in Eqs. 2.1 and 2.3. Here, the added volume is reset to zero not after division, butafter initiation.In [30], a strong size control acting on initiation initiation was proposed where all cells will have inititatedDNA replication before reaching some fixed volume x i . This hypothesis can be implemented in our initiation-adder model by setting k i ( x → x i , y, t ) → ∞ . The probability that a cell born at time t has not yet initiated, e − R tt k i ( x ( s ) ,y ( s ) ,s ) d s , always vanishes for all ( t , x t , y t ) before some finite time t and x ( t ) < x i . Thus, n ( x, , t )is nonzero only in [0 , x i ] for all t . If there exists a constant τ such that lim τ → τ e − R t τt k i ( x ( s ) ,y ( s ) ,s ) d s = 0 for all t , then the largest volume that any cell can attain will be e λτ x i , leading to strict size control and no blowup.Fig. 8 shows numerical solutions to Eq. 3.9 using the same birth rate function as that used in Fig. 3(d-f). Note that due to cell size control affecting the pre-initiation stage, initial daughter cell sizes stay small atinitiation and n ( x, y, t ) is more peaked near y ≈ x .If one takes k i sufficiently large, both daughter cells will nearly instantly initiate DNA replication afterdivision. We have checked numerically that for constant k i = 10 , the densities n ( x, y, t ) are negligible while n ( x, y, t ) approaches the density of the division adder shown in Fig. 3 (for the same differential division functions - -- -- Fig. 8: Normalized densities of pre-initiation cell populations ¯ n and post-initiation cell populations ¯ n for atvarious fixed times t = 1 , ,
12. Here, we used k i ( x ) = p ( x ) / (cid:2) − R x p ( x ′ )d x ′ (cid:3) with p ( x ) ∼ N (1 , .
1) and thesame ˜ β ( x, y, z, t ) as that used in Fig. 3(d-f). (a-c) shows the normalized densities ¯ n ( x, y, t ) ≡ n ( x, y, t ) /N ( t )where N ( t ) = R d y R d x ( n + n ). (d-f) shows the normalized post-initiation density ¯ n ( x, y, t ). For the k i used in this example, the pre-initiation densities span larger volume and added volumes. The densities areindistinguishable from those at steady state after about t = 2.˜ β ). Thus, the initiation adder model converges to the standard division adder model when k i → ∞ . This canbe seen from the first of Eqs. 3.9 where n can be neglected and is dominated by the two terms on the right-hand-side. Substituting k i ( x, y, t ) n ≈ R ∞ x d zx n ( z, yz/x, t ) into the integral terms in the second equation, wefind Eq. 2.1 for n ( x, y, t ).
4. Summary and Conclusions.
In this paper, we used PDE models to describe population dynamicalbehavior under the adder division mechanism. Under certain conditions, this PDE for the adder mechanism canalso be converted to the well-known size- and age-structured PDE. In the absence of death, we motivated modelsfor the differential birth rate function ˜ β ( x, y, z, t ) that are consistent with normalized division probabilities InAppendix A we showed existence and uniqueness of a weak solution to the PDE model within a time interval[0 , T ] during which the solution’s support can be bounded. One can prove similar results when both time andspace are unbounded as this problem is related to other first order PDE models that have been studied in moredetail.With a weak solution justified, we explored the “adder-sizer” PDE via numerical experiments and Monte-Carlo simulations of the underlying stochastic process. Our results show that event-based Monte-Carlo simu-lations of the discrete process generate sample configurations. The observed configurations are are consistentwith samples from the cell densities numerically computed from our PDE model.When broader differential division rates are used (when cell division is more asymmetric), we find, underthe same initial conditions, a broader cell density n ( x, y, t ) and a broader event density ρ D ( x, y, T ). We alsodemonstrate numerically, the divergence of the mean cell size h x ( t ) i = M ( t ) /N ( t ). We showed that divisionprobabilities that are broader in the age or added size (and smaller in magnitude) more likely lead to mean cellsizes that explode with time.We then incorporated growth rate correlation between cells of successive generations [17] into our “adder-sizer” PDE model. By extending the dimension of the density function to include growth rates and allowing for ariability in growth rate as new cells are born, we developed a PDE model that incorporated the stochasticnature of growth rate inheritance and that describes evolution of the growth rate distribution of cells. We foundthat the steady-state value of the mean growth rate depends on the correlation of growth rates between motherand daughter cells. This dependence arises from a subtle interaction between the shape of the growth ratedistribution and the distribution of variations in the growth rate from one generation to the next.Finally, we proposed a coupled partial integro-differential equation (PIDE) to model two-phase cell popu-lation dynamics under a new initiation-adder mechanism suggested by recent experimental results. In the limitthat the initiation rate k i of DNA replication is significantly faster than all other time scales in the problem, thenumerical solutions of the initiation adder model (Eq. 3.9) converge to those of division adder model (Eqs. 2.1and 2.3). Under proper assumptions that come from experimental findings, we found that the initiation adderwould also lead to effective cell size control [30].There are new cellular processes and size control mechanisms that have been recently discovered and thatcan be mathematically modeled. Thus, there is likely general mathematical topics that remain to be exploredwithin PDE and PIDE models of structured populations. For example, a recent experimental study indicatesthat an adder mechanism may be the result of several consecutive processes in the cell division cycle, suggestingthat a much more complicated coupled system of PDEs/PIDEs would be required. Acknowledgments.
This work was supported in part by grants from the NSF (DMS-1814364), the Na-tional Institutes of Health (R01HL146552) and the Army Research Office (W911NF-18-1-0345).
Appendix A. Existence and uniqueness of a weak solution for the adder-sizer model.
In this section we show the existence and uniqueness of the solution to the “adder-sizer” model PDE. Thefull problem is defined as(A.1) ∂n∂t + ∂ ( ng ) ∂x + ∂ ( ng ) ∂y = − β ( x, y, t ) n ( x, y, t ) ,g ( x, , t ) n ( x, , t ) = 2 Z ∞ x d x ′ Z x ′ d y ˜ β ( x ′ , y, x, t ) n ( x ′ , y, t ) ,β ( x, y, t ) := Z x ˜ β ( x, y, z, t )d z, ˜ β ( x, y, z ′ , t ) = ˜ β ( x, y, z − z ′ , t ) , ˜ β ( x, y, , t ) = 0 , n ( x, x, t ) = 0 ,n ( x, y, t = 0) := n ( x, y ) . where the independent variables ( x, y, t ) ∈ R ∩ { y < x } × R + .First, we assume that(A.2) 0 < g min ≤ g ∈ C ( { ( R + ) ∩ { y ≤ x }} × R + ) ,n ( x, y ) ∈ L ∩ L ∞ ∩ C ( R + ∩ { y < x } ) , ≤ ˜ β ∈ L ∞ ∩ L ∩ C ( { ( R + ) ∩ { y < x, z < x }} × R + ) ,β ( x, y, t ) ∈ L ∞ ∩ L ∩ C ( { ( R + ) ∩ { y ≤ x }} × R + ) , and nondimensionalize the size and added size by ∆, the added size parameter defined in Eq. 2.9. We alsoimpose an additional assumption on g :(A.3) | g ( x, y, t ) | < K ( t + x + 1) , K < ∞ . We also assume the initial distribution n ( x, y ) is compactly supported in (0 , Ω) × [0 , Ω) , Ω < ∞ . From thisassumption and A.3, the closure of n ( x, y, T )’s support is compact for any finite time T since n = 0 onlywhen y < x and x ( s ) ≤ Ce Ks − (1 + T ) from Gr¨onwall’s Inequality, where C < T + Ω is given by theinitial condition. At any finite time T , the support of n ( x, y, T ) is bounded and we assume it is contained in[0 , Ω( T )) × [0 , Ω( T )). Furthermore, by setting g, β, ˜ β = 0 at the given time T when ( x, y ) is out of the supportof n , we can assume the closure of g, β, ˜ β ’s support to be compact. One can generalize the definition of theweak solution n to [( R + ) ∩ { y < x } ] × [0 , ∞ ) as in [22]. efinition A.1 Given time
T < ∞ and assuming A.2, for a function n ∈ L ((([0 , Ω( T )]) ∩ { y < x } ) × [0 , T ]) , Ω( T ) < ∞ with n ( x, y, t ) = 0 in [0 , Ω( T )) × [0 , Ω( T )) , y < x, t ∈ [0 , T ], n is said to satisfy the adder-sizerPDE in the weak sense in time [0 , T ], if(A.4) − Z T d t Z ∞ d x Z x d y n ( x, y, t ) (cid:20) ∂ Ψ ∂t + g ( x, y, t ) ∂ Ψ ∂x + g ( x, y, t ) ∂ Ψ ∂y − β ( x, y, t )Ψ( x, y, t ) (cid:21) = Z ∞ d x Z x d y n ( x, y )Ψ ( x, y ) + Z T d t Z ∞ d x Ψ( x, , t ) n ( x, , t ) g ( x, , t )holds for all test functions Ψ ∈ C (([0 , Ω( T )]) ∩ { y ≤ x } ) × [0 , T ]) satisfying Ψ( x, y, T ) ≡ , Ψ(Ω( T ) , y, t ) = 0and Ψ( x, x, t ) = 0, where we set g, ˜ β, β = 0 for x ≥ Ω( T ) , x ≤ y or x ≤ z . Upon using the boundary conditionin A.1, the right-hand-side becomes Z ∞ d x Z x d y n ( x, y )Ψ ( x, y ) + 2 Z T d t Z ∞ d x Z x d y Z x d z Ψ( z, , t ) ˜ β ( x, y, z, t ) n ( x, y, t ) . Note that if n ∈ C ((( R + ) ∩ { y < x } ) × R + ) is a classical solution to the PDE (Eq. A.1), then it mustalso satisfy Eq. A.4 in any time interval [0 , T ]. We refer to [22] for a proof of the existence and uniqueness of aweak solution of a related, simpler renewal equation. However, our adder-sizer PDE is more complicated. Theproof of uniqueness requires very different techniques from the sizer PDE; yet the proof of existence is similarto the proof in [22]. A.1. Uniqueness.
First, we prove uniqueness of the solution to A.4. Assume there are two weak solutions n (0) and n (1) for the adder-sizer PDE satisfying A.4 with the same initial condition n (0)0 ( x, y ) = n (1)0 ( x, y ). Takingthe difference between using these purported solutions, we obtain(A.5) − Z T d t Z ∞ d x Z x d y ∆ n ( x, y, t ) (cid:20) ∂ Ψ ∂t + g ( x, y, t ) ∂ Ψ ∂x + g ( x, y, t ) ∂ Ψ ∂y − β ( x, y, t )Ψ( x, y, t ) (cid:21) = 2 Z T d t Z ∞ d x Z x d y Z x d z Ψ( z, , t ) ˜ β ( x, y, z, t )∆ n ( x, y, t ) , where ∆ n = n (1) − n (0) . A.1.1. Adjoint Problem.
We consider the adjoint problem for Ψ in the given time interval [0 , T ] andwith a source term S ( x, y, t ): (A.6) ∂ Ψ ∂t + g ( x, y, t ) ∂ Ψ ∂x + g ( x, y, t ) ∂ Ψ ∂x − β ( x, y, t )Ψ( x, y, t ) = − Z x Ψ( z, , t ) ˜ β ( x, y, z, t )d z − S ( x, y, t ) , ≤ y < x Ψ( x, y, T ) = 0 , Ψ(Ω( T ) , y, t ) = 0 , Ψ( x, x, t ) = 0 . Theorem A.1
Assume A.2, and S ∈ C ([0 , Ω( T )] × [0 , T ]) , S (Ω( T ) , y, t ) = 0, and S = 0 when x ≤ y .Then there exists a unique C solution to the adjoint problem. Proof.
We can transform the above equation into an ODE along the characteristic line and use contractionmapping, which is a standard practice in functional analysis to prove existence and uniqueness of the solutionto a PDE problem. On the left-hand-side of Eq. A.6, we apply the characteristic line method. Setting X ( c, t ) =( x ( c, t ) , y ( c, t )) on the characteristic lines leads to ∂X ( c, s ) ∂s = ( g ( x, y, s ) , g ( x, y, s )) , t ≤ s ≤ T,X ( c, t ) = ( x t , y t ) , ≤ y t < x t , x t − y t = c. Since we have x ( s ) − y ( s ) = x t − y t , the above equation can be simplified to ∂X ( c, s ) ∂s = ˜ g ( X ( c, s ) , s ) , x ( c, t ) = x t , y ( c, t ) = x t − c here ˜ g ( X ( c, s ) , s ) = ( g ( x ( c, s ) , x ( c, s ) − c, s ) , g ( x ( c, s ) , x ( c, s ) − c, s )). Once c is fixed and x t is given, the aboveequation becomes an ordinary differential equation. Given x t , we define ( ˜Ψ( c, s ) := Ψ( X ( c, s ) , s ) e − R st β ( X ( c,v ) ,v ) dv ,U ( c, z, s ) := 2 ˜ β ( X ( c, s ) , z, s ) e − R st β ( X ( c,v ) ,v ) dv , ˜ S ( c, s ) := S ( X ( c, s ) , s ) e − R st β ( X ( c,v ) ,v ) dv . Thus, along the characteristic line we can write A.6 as(A.7) ∂∂s ˜Ψ( c, s ) = − Z x ( c,s )0 Ψ( z, , s ) U ( c, z, s )d z − ˜ S ( c, s ) . Since ˜Ψ( c, T ) = 0 and ˜Ψ( c, t ) = Ψ( x t , x t − c, t ),(A.8) Ψ( x t , x t − c, t ) = Z Tt ˜ S ( c, s )d s + Z Tt d s Z x ( c,s )0 d z Ψ( z, , s ) U ( c, z, s ) , < c ≤ x t . We can see that if x ≤ y or x t ≥ Ω( T ), Ψ( t, x t , x t − c ) = Ψ( t, x, x ) = 0 since U, ˜ S = 0 for c ≤ x t > Ω( T ).Using c = x t , Eq. A.8 becomes(A.9) Ψ( x t , , t ) = Z Tt ˜ S ( x t , s )d s + Z Tt d s Z x ( x t ,s )0 d z Ψ( z, , s ) U ( x t , z, s ) . From condition A.3 we obtain x ( s ) ≤ ( x t + 1 + T ) e K ( s − t ) − (1 + T ). From condition A.3, we define ˜ B = 2 k ˜ β k ∞ < ∞ . Next, we choose s = max { T − K ln(1+
12 ˜ B (1+ T ) ) , T − K ln 2 , T − } such that e K ( T − t ) ≤
12 ˜ B (1+ T ) , s ≤ t ≤ T ,and choose x s small enough such that x s < min { ,
18 ˜ B ( T − s ) } . We denote a mapping T defined on the functionalspace as T (Ψ)( x t , , t ) = Z Tt ˜ S ( x t , s )d s + Z Tt d s Z x ( s,x t )0 d z Ψ( z, , s ) U ( x t , z, s ) , t ∈ [ s, T ] , x t ∈ [0 , x s ] . It is easy to verify that T is a contraction mapping for Ψ( x t , , t ) and thus there exists a unique solution Ψ satisfying A.6 in D defined as D = { ( x, t ) | s ≤ t ≤ T, ≤ x ≤ x ( x s , t ) } . We then let x s > x s and define D = { ( x, t ) | s ≤ t ≤ T, ≤ x ≤ x ( x s , t ) } such that the difference of the area between regions D and D is lessthan ˜ B − . Next, define a second mapping T by T (Ψ)( x t , , t ) = Z Tt d s Z x ( x t ,s ) x ( x s ,s ) d z Ψ( z, , s ) U ( x t , z, s ) + I ( x s , t ) , t ∈ [ s, T ] , x t ∈ [ x ( t, x s ) , x s ] ,I ( x s , t ) = Z Tt d s ˜ S ( x t , s ) + Z Tt d s Z x ( x s ,s )0 d z Ψ ( z, , s ) U ( x t , z, s ) .T is also a contraction mapping and we can obtain a Ψ on D such that T (Ψ ) = Ψ . Denote(A.10) Ψ( x, , t ) = ( Ψ ( x, , t ) , ( x, t ) ∈ D , Ψ ( x, , t ) , ( x, t ) ∈ D , and it is easy to verify that Ψ is C continuous on D ∩ D by first proving it is continuous and then taking thepartial derivatives, and Ψ satisfy A.6 in the region D ∪ D .Following the same procedure, we can extend Ψ to satisfy A.6 in the region t ∈ [ s, T ]. Then, for [0 , s ], wechoose a ˜ s close enough to s and use the same strategy by defining T as A.11) T (Ψ)( x t , , t ) = Z st d r ˜ S ( x t , r ) + Z st d r Z x ( x t ,r )0 d z Ψ( z, , r ) U ( x t , z, r ) + ˜ I ( t, x s ) , t ∈ [˜ s, s ] , ˜ I ( x s , t ) = Z Ts d r ˜ S ( x t , r ) + Z Ts d r Z x ( x t ,r )0 d z Ψ( z, , r ) U ( x t , z, r ) . We finally obtain a unique function Ψ satisfying A.6 in [0 , T ] × [0 , ∞ ).From A.8, the value of Ψ is determined by ˜ S, Ψ( x, , t ) , U and we conclude that there exists a unique C solution for A.6. A.1.2. Uniqueness of weak solution for the adder-sizer model.
From Section A.1.1 we obtain theexistence and uniqueness of Ψ of the adjoint problem. Given any time T and S ( x, y, t ) ∈ C ( R + × ( R + ) )satisfying the condition in Theorem A.1, since we can set g, β, ˜ β ’s support to be compact in [0 , T ], we can finda unique C continuous Ψ satisfying A.6. By substituting A.6 into A.5, we obtain(A.12) Z T d t Z Ω( T )0 d x Z x d y ∆ n ( x, y, t ) S ( x, y, t ) = 0for any S ( x, y, t ) ∈ C ( R + × R +2 ) satisfying S ( x ≤ y, t ) = S ( x ≥ Ω( T ) , y, t ) = 0, which implies n ≡ y < x ≤ Ω( T ). So at any given time T the weak solution, if exists, is unique.One can also set the condition for ˜ β, g weaker even when we define the weak solution in unbounded region[0 , ∞ ) × ( R + ) ∩{ y < x } . In [22] such work is done for the renewal equation. We do not discuss this generalizationin detail here. A.2. Existence of the weak solution.
We construct a series of functions { n i } with a limit n for thisseries satisfying A.6 for all test functions Ψ. We use a semi-discrete approximation to discretize the PDE andobtain piecewise solutions. As the mesh size becomes smaller, we expect the piecewise solution to converge toa function n satisfying A.4. The idea of constructing a series of piecewise constant solutions and proving theirconvergence to a weak solution is similar to that in [22]. A.2.1. Semi-discrete approximation for the PDE.
We choose a uniform grid with mesh size h > x and y axis and let time t be continuous. We denote(A.13) ( x i , y j ) = ( ih, jh ) , ( x i + , y j + ) = (( i + 12 ) h, ( j + 12 ) h ) , j < i ∈ N ,β i + ,j + ( t ) = 1 h Z ( i +1) hih d y Z ( j +1) hjh d x β ( x, y, t ) , j < i ∈ N , ˜ β i + ,j + (( s + 12 ) h, t ) = 1 h Z ( i +1) hih d z Z ( j +1) hjh d y Z ( s +1) hsh d x ˜ β ( x, y, z, t ) , s ≤ i,g i,j ( t ) = g ( ih, jh, t ) , j < i ∈ N . Here, β i + ,j + ( t ) = h P is =0 ˜ β i + ,j + (( s + ) h, t ). Given a fixed time T , we wish to find a solution ofpointwise function n h ( t ), which takes values on the grid points ( x i + , y j + ). Then n h can be seen as a vectorfunction. According to our assumption there exists Ω such that the initial value n is nonzero within the region { ( x, y ) | y < x, x < Ω } , and from our previous calculation there exists Ω( T ) < ∞ such that n is nonzero withinthe region { ( x, y ) | y < x, x < Ω( T ) } . Eventually, we will set h ( k ) = Ω( T ) /k and let the mesh size h → k → ∞ .By discretizing Eqs. A.1, we expect the vector function n h ( t ) to satisfy the below equations for t ∈ [0 , T ]and 0 < j < i < L ( L is the number of discretization grid points along one direction): A.14) h d n i + ,j + ( t )d t + h ( g i +1 ,j + ( t ) n i + ,j + ( t ) − g i,j + ( t ) n i − ,j + ( t ))+ h ( g i + ,j +1 ( t ) n i + ,j + ( t ) − g i + ,j ( t ) n i + ,j − ( t )) + h β i + ,j + ( t ) n i + ,j + ( t ) = 0 , ≤ j < i − h d n i + ,j + ( t )d t + hg i +1 ,j + ( t ) n i +1 ,j + ( t ) − hg i + ,j ( t ) n i + ,j − ( t ) + h β i + ,j + ( t ) n i + ,j + ( t ) = 0 , ≤ j = i − g i + , ( t ) n i + , − ( t ) = 2 h L − X ℓ = i ℓ − X j =0 ˜ β ℓ + ,j + (( i + 12 ) h, t ) n ℓ + ,j + ( t ) ,n i + ,j + (0) = 1 h Z x i +1 x i d y Z y j +1 y j d x n ( x, y ) , n i + ,i + ( t ) = 0 , where we henceforth omit the h superscript in the proof. In the two-dimensional upwind scheme, derivatives inone direction are neglected on neighboring sites in the other direction: n i,j ± = n i − ,j ± , n i ± ,j = n i ± ,j − .The boundary condition n ( x, x, t ) = 0 is implemented by n i + ,i + ( t ) = 0 for any t and i .We will obtain a uniform bound irrelevant of h for n . All coefficients in the above ODE equations are C continuous, which means that there exists a unique solution in time [0 , T ] , T < ∞ . Theorem A.2
For t ∈ [0 , T ] and assuming A.2 holds, we find the bound(A.15) L − X i =1 i X j =0 | n i + ,j + ( t ) | ≤ e Mt L − X i =1 i X j =0 | n i + ,j + (0) | , where ˜ B = 2 k ˜ β k ∞ , M = 2 B − b, B = k β k ∞ , and b = min t min i,j β i + ,j + ( t ). The L ∞ bound is given by k n h ( t ) k ∞ ≤ max { g min ˜ Be MT k n (0) k , k n h (0) k ∞ } e g ′ t where ˜ g ′ is the L ∞ bound of ∂g/∂x, ∂g/∂y . Proof
For the summation of n over all grid points, we multiply the first equation in A.14 by sign( n i + ,j + )for each i, j ≤ i ,(A.16) h dd t | n i + ,j + ( t ) | + hg i +1 ,j + ( t ) | n i + ,j + ( t ) | + hg i + ,j +1 ( t ) | n i + ,j + ( t ) | + h β i + ,j + ( t ) | n i + ,j + ( t ) | ≤ hg i,j + ( t ) | n i − ,j + ( t ) | + hg i + ,j ( t ) | n i + ,j − ( t ) | By multiplying the second equation in A.14 by sign( n i + ,j + ) for each i, j ≤ i pair and summing over index P L − i =1 P i − j =0 , h L − X i =1 i − X j =0 | n i + ,j + ( t ) | + h i − X j =0 g L,j + ( t ) | n L − ,j + ( t ) | + h L − X i =1 i − X j =0 β i + ,j + ( t ) | n i + ,j + ( t ) | ≤ h L − X i =0 g i + , ( t ) | n i + , − ( t ) | . We can simplify the above expression to h ddt L − X i =1 i − X j =0 | n i + ,j + ( t ) | + h L − X i =1 i − X j =0 β i + ,j + | n i + ,j + ( t ) | h L − X i =0 (cid:12)(cid:12) L − X ℓ = i ℓ − X j =0 ˜ β ℓ + ,j + (( i + 1 / h, t ) n ℓ + ,j + ( t ) (cid:12)(cid:12) ≤ h L − X ℓ =1 ℓ − X j =0 | β ℓ + ,j + ( t ) || n ℓ + ,j + ( t ) | . We then have dd t L − X i =1 i − X j =0 | n i + ,j + ( t ) | ≤ (2 B − b ) L − X i =1 i − X j =0 | n i + ,j + ( t ) | , which yields(A.17) L − X i =1 i − X j =0 | n i + ,j + ( t ) | ≤ e Mt L − X i =1 i − X j =0 | n i + ,j + (0) | . A.17 states that the l norm of all the values on the grid points are uniformly bounded and independent of h .Next, we estimate the L ∞ bound of n h . First, we consider j = 0 and assume S ( t ) = max ≤ i ≤ L − | n i + , ( t ) | e − ˜ g ′ t for t ∈ [0 , T ]. For the maximum value of S at some index i , we find h d | n i + 12 , ( t ) | d t + h ( g i +1 , ( t ) | n i + , ( t ) | − g i, ( t ) | n i − , ( t ) | ) + h ( g i + , ( t ) | n i + , ( t ) | − g i + , ( t ) | n i + , − ( t ) | ) ≤ ,h d | n i + 12 , ( t ) | d t + hg i +1 , ( t ) | n i + , ( t ) | − g i + , ( t ) | n i + , − ( t ) | ≤ , i = 1 . and d( | n i + , ( t ) | e − ˜ g ′ t )d t + h − g i + , ( t ) | n i + , ( t ) | e − ˜ g ′ t ≤ h − g i + , ( t ) | n i + , − ( t ) | e − ˜ g ′ t , By the assumption that g ( x, y, t ) ≥ g min ( t ) ≥ g min > g < K ( T + 1 + Ω( T )), we have(A.18) d( | n i + , ( t ) | e − ˜ g ′ t )d t + h − g min ( t ) | n i + , ( t ) | e − ˜ g ′ t ≤ h − (cid:18) g min ( t ) g min (cid:19) max ≤ i ≤ L − | g i + , ( t ) n i + , − ( t ) | . Finally, defining G ( t ) = h − R t g min ( s )d s yieldsd( | n i + , ( t ) | e − ˜ g ′ t e G ( t ) )d t ≤ h (cid:18) g min ( t ) g min (cid:19) max ≤ i ≤ L − | g i + , ( t ) n i + , − ( t ) | e G ( t ) . From the L bound, we can deducemax t max ≤ i ≤ L − | g i + ( t ) n i + , − ( t ) | ≤ h ˜ Be MT k n h (0) k ≤ ˜ Be MT k n (0) k , t > S ( t ) e G ( t ) (A.19) S ( t ) e G ( t ) ≤ S (0) + 1 g min ˜ Be MT k n (0) k ( e G ( t ) − , and S ( t ) ≤ max ≤ i ≤ L − { n i + , (0) , g min ˜ Be MT k n (0) k } , which then gives the L ∞ bound for the pointwise solution n h when j = 0. ow, we estimate | n i + ,j + ( t ) | by first defining P ( t ) ≡ max ≤ i ≤ L − , ≤ j ≤ i − {| n i + ,j + ( t ) | e − g ′ t } . At a fixedtime t , specific values of i and j define P ( t ). If the maximum occurs at j = 0, P ( t ) = S ( t ) e − ˜ g ′ t . If the maximumoccurs at i − > j >
0, we have h dd t ( | n i + ,j + ( t ) | e − g ′ t ) = − h g i,j + ( t ) | n i + ,j + ( t ) | − g i +1 ,j + ( t ) | n i + ,j + ( t ) | + g i + ,j ( t ) | n i + ,j + ( t ) | − g i + ,j +1 ( t ) | n i + ,j + ( t ) | + 2 h ˜ g ′ | n i + ,j + ( t ) | i e − g ′ t ≤ , (A.20)while if the maximum occurs at j = i − >
0, we have(A.21) dd t ( | n i + ,j + ( t ) | e − g ′ t ) ≤ h h − (cid:16) g i + ,j ( t ) − g i +1 ,j + ( t ) (cid:17) | n i + , ( t ) | − g ′ | n i + ,j + ( t ) | i e − g ′ t ≤ . In Eqs. A.20 and A.21, i, j are the maximizing indices that define P ( t ).For any t ∈ (0 , T ] we can find a minimum ˜ t < t such that P ( v ) > S ( v ) e − ˜ g ′ v for v ∈ (˜ t, t ]. If ˜ t = 0, and since P ( t ) is nonincreasing from Eq. A.21, P ( t ) ≤ P (0) = k n h (0) k ∞ . If t > ˜ t > P ( t ) ≤ P (˜ t ) ≤ S (˜ t ) ≤ max ≤ t ≤ T S ( t ),while if ˜ t = t , P ( t ) = S ( t ) ≤ max ≤ t ≤ T S ( t ). Thus, P ( t ) = k n h ( t ) k ∞ e − g ′ t ≤ max { max ≤ t ≤ T { S ( t ) } , k n h (0) k ∞ } and(A.22) k n h ( t ) k ∞ ≤ max { max ≤ t ≤ T { S ( t ) } , k n h (0) k ∞ } e g ′ t , giving the second conclusion in Theorem A.2 that the L ∞ bound is uniform and independent of h . A.2.2. Existence of the weak solution.
For a given time
T < ∞ , we can take the grid size h ( k ) =Ω( T ) /k → k → ∞ . Spatially piecewise constant functions can then be defined based onthe sequence of vector functions { n h ( k ) } . By setting n hi + ,i + ( t ) = 0, we define n h ( x, y, t ), β h , and ˜ β h as n h ( x, y, t ) = k − X i =0 i − X j =0 n hi + ,j + ( t ) ( ih ≤ x < ( i + 1) h, jh ≤ y < ( j + 1) h ) ,β h ( x, y, t ) = k − X i =0 i X j =0 β i + ,j + ( t ) ( ih ≤ x < ( i + 1) h, jh ≤ y < ( j + 1) h ) , ˜ β h ( x, y, z, t ) = k − X i =0 i − X j =0 i − X ℓ =0 ˜ β i + ,j + (( ℓ + 12 ) h, t ) ( ih ≤ x < ( i + 1) h, jh ≤ y < ( j + 1) h, ℓh ≤ z < ( ℓ + 1) h ) ,n h ( x, , t ) = n hi + , − ( t ) , ih ≤ x < ( i + 1) h, where above, h = h ( k ) and is the indicator function. Since there is an upper bound for both β and ˜ β , andboth β, ˜ β are continuous, we have the following resultlim k →∞ β h ( k ) ( x, y, t ) → β ( x, y, t ) a.e. ≤ β h ( k ) ≤ k β k ∞ < ∞ , lim k →∞ ˜ β h ( k ) ( x, y, z, t ) → β ( x, y, z, t ) a.e. ≤ ˜ β h ( k ) ≤ k ˜ β k ∞ < ∞ , lim k →∞ n h ( k ) ( x, y, → n ( x, y, a.e.. Then, we can apply Theorem A.2 to the piecewise constant solutions n h ( k ) of Eqs. A.14. Corollary A.3
Under the conditions of Theorem A.2, for any t ∈ [0 , T ] and any h ,(A.23) Z Ω( T )0 d y Z Ω( T )0 d x | n h ( x, y, t ) | ≤ e Mt Z Ω(0)0 d y Z Ω(0)0 d x | n h ( x, y, | nd(A.24) k n h ( t ) k ∞ ≤ max {| n (0) | ∞ , Be MT | n (0) | } e g ′ t , where B, M, ˜ g ′ are defined in Theorem A.2. The proof is the direct consequence of Theorem A.2.The sequence of piecewise constant functions { n h ( k ) } is uniformly bounded and n h ( k ) ∈ L ∩ L ∞ ([0 , Ω( T )] ∩{ y < x } × [0 , T )), so n h ( k ) are all L functions. There exists a function n ∈ L ([0 , Ω( T )] ∩ { y < x } × [0 , T )) anda subsequence k i → ∞ that satisfies n h ( k i ) ⇀ n . Since L ([0 , Ω( T )] ∩ { y < x } × [0 , T )) implies L integrability,we can deduce that n is an L function as desired.To prove n h ( k i ) ⇀ n , we need only to verify that there exists a subsequence n h ( k i ) such that for all testfunctions f ∈ L , R T d t R Ω( T )0 d x R x d y n h ( k i ) f → R T d t R Ω( T )0 d x R x d y nf . Since L space is separable, we have acountable set of basis function { b i ( x, y, t ) } for the space L ([0 , Ω( T )] ∩ { y < x } × [0 , T )). Thus, every n h ( k ) can be decomposed as n h ( k ) = P ∞ i =1 α ki b i . The sequence { n h ( k ) } is uniformly L ∞ bounded, so P α k are alluniformly bounded. We can then select a subsequence { n h ( k i ) } from { n h ( k ) } satisfying lim i →∞ α k i j = α j so that P ∞ i =1 α j < ∞ . If we set n = P ∞ i =1 α i b i , then, by decomposing any test function Ψ ∈ L ([0 , Ω( T )] ∩ { y 0) and − R T d t R ∞ d x R x d y R x d z Ψ( z, , t ) ˜ β ( x, y, z, t ) n ( x, y, t ), respectively.The third term on the RHS of A.26 h R T d t P L − i =1 g i + ,i ( t ) n hi + ,i − ( t )Ψ i + ,i − ( t ) tends to 0 since Ψ is C continuous and is 0 on the boundary x = y . Since Ψ is continuous and is 0 at x = Ω( T ), h R T P L − j =0 g L,j + ( t ) n hL − ,j + ( t )Ψ L − ,j + ( t )d t → h → 0. Finally, the last term on the RHS of of A.26 h R T P L − i =1 P i − j =0 h β i + ,j + ( t ) n hi + ,j + ( t )Ψ i + ,j + ( t )d t → R T d t R Ω( T )0 d x R x d y β ( x, y, t ) n ( x, y, t )Ψ( x, y, t ).By passing to the limit h → 0, we conclude that n exactly satisfies the condition of a weak solution in A.4.Since the numerical solution obtained by the scheme in Appendix B is a discretization in time for the ODEsystem A.14 it is an approximation to the solution of A.14. Provided h, ∆ t → k g k ∞ ∆ t < h , conclude that at least a subsequence of the numerical solutions converge to the unique weaksolution of A.1. Furthermore, recently, the existence to an eigenpair of the adder-sizer PDE A.1 under specificsmooth conditions satisfied by the coefficients g, β, ˜ β has been proved in [10], allowing for studying asymptoticbehavior of the solution. Appendix B. Numerical Scheme. We denote u ( t ) = { n ( t ) , n ( t ) , . . . , n L − ( t ) } T where n j ( t ) = { n ,j − , n ,j − , . . . , n L − ,j − } and n i ≤ j = 0. Equations 2.17 and 2.18 can then be written in the form u ( t + ∆ t ) = A ( t ) u ( t ), where(B.1) A ( t ) = B + C C C C · · · C L − C L − D B · · · D B · · · · · · B L − 00 0 0 0 · · · D L − B L − , is made up of the following L − L × L matrices i = i × i ) 0 ( i × ( L − i ))0 (( L − i ) × i ) b i , C i = × i ) 0 (1 × ( L − i ))0 (( L − × i ) c i , and D i = i × i ) 0 ( i × ( L − i ))0 (( L − i ) × i ) d i , in which b i is a lower bidiagonal matrix with diagonal(B.2) diag( b i ) = 1 − h g j +1 ,i − ( t )d t − h g j + ,i ( t )d t − β j + ,i − ( t )d t, j = i, i + 1 , . . . , L − , and lower off diagonal ( b i ) − = g j,i − ( t ) d th , j = i + 1 , . . . , L − , (B.3) ( c i ) sj = (cid:26) ˜ β i − + j,i − (( s + ) h, t )d t, i + j − s − > , i + j ≤ L . and d i is a diagonal matrix diag( d i ) = g j + ,i − ( t ) d th , j = i, i + 1 , . . . , L − . Appendix C. Monte-Carlo Simulations. In this section we describe the implementation of our Monte-Carlo simulations of the process underlying the adder-sizer mechanism. Suppose we have a list of cells at time t denoted by S ( t ) = { c ( x i , y i , t, b ) , ..., c i ( x i , y i , t, b i ) } , where x i is cell c i ’s volume and y i is its added volume.The cell’s division factor b i is determined at birth, which is drawn from a uniform distribution U (0 , β of the form 2.10 and ˜ β of the form 2.12. We set the maximum allowable time step to∆ t = 0 . 01 and determine the next state of the system at time t ′ by the following • Step 1: For each cell i , calculate its age a i at time t by the exponential growth law d x d t = λx . We requirethat G i = R a i γ ( a ′ )d a ′ < b i at the beginning of each step for every i . • Step 2: For each cell, calculate G i = R a i +∆ t γ ( a ′ )d a ′ . If G i ≥ b i , then we numerical calculate a ∆ t i such that R a i +∆ t i γ ( a ′ )d a ′ ≈ b i . • Step 3: Choose the smallest ∆ t i among all possible ∆ t i s as the new time step, set time t ′ = t + ∆ t i andlet all cells gain an extra volume λx i ∆ t i . If there is no such ∆ t i , which means G i < b i for every i , goto step 5. • Step 4: Remove cell i from S ( t ′ ), record its volume x at t ′ , and generate the random numbers r from thedistribution h ( r ) and b m , b m +1 from U (0 , S ( t ′ ) labeled by c m ( rx, , t, b m )and c m +1 ( x − rx, , t, b m +1 ). • Step 5: If G i < b i for all i , set t = t ′ and let all cells gain an extra volume λx i ∆ t i . • Step 6: Return to step 1 until t ′ > t max , the maximum time of the simulation.Here, we set the initial added volume of all cells to zero so the condition in step 1 above is automaticallysatisfied at t = 0. For our runs, we used 10 cells of initial volume 0 . t max = T is the same as the maximumtime for the numerical PDE experiments. We can also generalize the model to incorporate the mother-daughtergrowth coefficient correlation by including a new label λ i to each cell. REFERENCES[1] E. Bernard, M. Doumic, and P. Gabriel , Cyclic asymptotic behaviour of a population reproducing by fission into two equalparts , (2016). 222] J. Betschinger and J. A. Knoblich , Dare to be different: asymmetric cell division in drosophila, C. elegans and vertebrates ,Current Biology, 14 (2004), pp. R674–R685.[3] C. Cadart, S. Monnier, J. Grilli, P. J. S´aez, N. Srivastava, R. Attia, E. Terriac, B. Baum, M. Cosentino-Lagomarsino, and M. Piel , Size control in mammalian cells involves modulation of both growth rate and cell cycleduration , Nature Communications, 9 (2018), p. 3275.[4] D. Chandler-Brown, K. M. Schmoller, Y. Winetraub, and J. M. Skotheim , The Adder Phenomenon Emerges fromIndependent Control of Pre- and Post-Start Phases of the Budding Yeast Cell Cycle , Current Biology, 27 (2017), pp. 2774–2783.[5] T. Chou and C. D. Greenman , A Hierarchical Kinetic Theory of Birth, Death and Fission in Age-Structured InteractingPopulations , Journal of Statistical Physics, 164 (2016. PMCID: PMC3894939), pp. 49–76.[6] M. Delarue, D. Weissman, and O. Hallatschek , A simple molecular mechanism explains multiple patterns of cell-sizeregulation , PLoS ONE, 12 (2017), p. e0182633.[7] M. Doumic, M. Hoffmann, N. Krell, and L. Robert , Statistical estimation of a growth-fragmentation model observed ona genealogical tree , Bernoulli, 21 (2015), pp. 1760–1799.[8] M. Doumic, B. Perthame, and J. P. Zubelli , Numerical solution of an inverse problem in size-structured populationdynamics , Inverse Problems, 25 (2009), p. 045008.[9] H. V. Foerster , Some remarks on changing populations , Kinetics of Cellular Proliferation Grune and Stratton, (1959),pp. 382–407.[10] P. Gabriel and H. Martin , Steady distribution of the incremental model for bacteria proliferation , Networks & HeterogeneousMedia, 14 (2019), p. 149.[11] C. D. Greenman , A path integral approach to age dependent branching processes , Journal of Statistical Mechanics: Theoryand Experiment, 2017 (2017), p. 033101.[12] C. D. Greenman and T. Chou , Kinetic theory of age-structured stochastic birth-death processes , Physical Review E, 93(2016), p. 012112.[13] M. Guo, L. Y. Jan, and Y. N. Jan , Control of Daughter Cell Fates during Asymmetric Division: Interaction of Numb andNotch , Neuron, 17 (1996), pp. 27–41.[14] H. R. Horvitz and I. Herskowitz , Mechanisms of asymmetric cell division: two Bs or not two Bs, that is the question ,Cell, 68 (1992), pp. 237–255.[15] M. Iannelli , Mathematical theory of age-structured population dynamics , Giardini editori e stampatori in Pisa, (1995).[16] M. D. Jauffret and P. Gabriel , Eigenelements of a general aggeregation-fragmentation model , Mathematical Models andMethods in Applied Sciences, 20 (2010), pp. 757–783.[17] D. A. Kessler and S. Burov , Effective Potential for Cellular Size Control , arXiv:1701.01725, (2017).[18] J. Lin and A. Amir , The effects of stochasticity at the single-cell level and cell size control on the population growth , CellSyst, 5 (2017).[19] A. G. McKendrick , Applications of mathematics to medical problems , Proc. Edinburgh Math. Soc., 44 (1926), pp. 98–130.[20] J. A. J. Metz and O. Diekmann , The Dynamics of Physiologically Structured Populations , Springer, 1986.[21] S. Modi, C. A. Vargas-Garcia, K. R. Ghusinga, , and A. Singh , Analysis of Noise Mechanisms in Cell-Size Control ,Biophysical Journal, 112 (2017), pp. 2408–2418.[22] B. Perthame , Introduction to structured equations in biology , 2008.[23] L. Robert, M. Hoffmann, N. Krell, S. Aymerich, J. Robert, , and M. Doumic , Division in Escherichia coli is triggeredby a size-sensing rather than a timing mechanism , BMC Biology, 12 (2014), p. 17.[24] M. Schaechter, O. Maaloe, and N. O. Kjeldgaard , Dependency on medium and temperature of cell size and chemicalcomposition during balanced grown of Salmonella typhimurium. , Journal of General Microbiology, 19 (1958), pp. 592–606.[25] F. Si, G. Le Treut, J. T. Sauls, S. Vadia, P. A. Levin, and S. Jun , Mechanistic origin of cell-size control and homeostasisin bacteria , Current Biology, 29 (2019), pp. 1760–1770.[26] J. W. Sinko and W. Streifer , A New Model for Age-size Structure of a Population , Ecology, 48 (1967), pp. 910–918.[27] L. Sompayrac and O. Maaloe , Autorepressor model for control of DNA replication , Nat New Biol, 241 (1973), pp. 133–135.[28] S. Taheri-Araghi, S. Bradde, J. T. Sauls, N. S. Hill, P. A. Levin, J. Paulsson, M. Vergassola, and S. Jun , Cell-sizecontrol and homeostasis in bacteria , Current Biology, 25 (2015), pp. 385–391.[29] W. J. Voorn, L. J. Koppes, and N. B. Frover , Mathematics of cell division in Escherichia coli cell division: comparisonbetween sloppy-size and incremental-size kinetics , Current Topics in Mol. Genet., 1 (1993), pp. 187–194.[30]