{'fraction_non_alphanumeric': 0.06372623574144487, 'fraction_numerical': 0.013840304182509505, 'mean_word_length': 3.9241467850647265, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
PDE models of adder mechanisms in cellular proliferation

Cell division is a process that involves many biochemical steps and complex biophysical mechanisms. To simplify the understanding of what triggers cell division, three basic models that subsume more microscopic cellular processes associated with cell division have been proposed. Cells can divide based on the time elapsed since their birth, their size, and/or the volume added since their birth -- the timer, sizer, and adder models, respectively. Here, we propose unified adder-sizer models and investigate some of the properties of different adder processes arising in cellular proliferation. Although the adder-sizer model provides a direct way to model cell population structure, we illustrate how it is mathematically related to the well-known model in which cell division depends on age and size. Existence and uniqueness of weak solutions to our 2+1-dimensional PDE model are proved, leading to the convergence 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 correlations in cellular growth rates. Numerical experiments illustrating possible average cell volume blowup and the dynamical behavior of cell 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 the cell cycle.

1. Introduction. How cells regulate and maintain their sizes, as well as sizes of their appendages is a longstanding research topic in cell biology. Besides growth of an individual cell, the size distributions within a population of cells are also a quantity of interest. When considering proliferating cell populations, individual cell growth is interrupted by cell division events that generate smaller daughter cells. The biological mechanisms that control 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 wall formation [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 mechanisms of cell division have been proposed. Cell division can be governed by cell age a, cell volume x [26], or added volume since birth y [29,28]. The division mechanism employed by a type of cell may be interrogated by tracking the volumes x, added volumes y, and ages a during division events. Volume growth of an individual cell can be straightforwardly measured and can be modeled by an effective empirical law such asẋ = g(a, x, y, t). A commonly 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 timer model, in which the cell division rate depends only on age of the cell is described by the classic McKendrick equation for n(a, t), the expected density of cells at age a and time t [19,9]. The McKendrick "transport" equation for 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 agedependent birth and death rates, respectively. The associated boundary condition n(t, 0) = 2 We 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" PDE model that involves cell age and size as controlling parameters. The proof of the existence and uniqueness of a weak solution to the proposed three-dimensional "adder-sizer" PDE turns out to be more complex than the proof for the timer and/or sizer counterparts [22]. Our proof leads to the convergence of the numerical solutions to the adder-sizer PDE, allowing us to numerically evaluate the corresponding structured cell populations, facilitating further 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 results for n(x, y, t) and division-event densities.
Next, we propose an extension to the adder-sizer model that incorporates cellular growth rates that are correlated across successive generations. Changes in growth rates at the single-cell level have been explored using stochastic mapping methods [7,18]. By numerically solving the PDE, we found out that the populationaveraged 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 the latest "initiation adder" mechanism, which states that the added mechanism takes effect on the cell's size at initiation instead of division in [25]. In contrast to the single-PDE division adder model, a model describing the initiation adder mechanism requires two coupled PDEs.
To model cell size control, stochastic maps that relate daughter cell sizes to mother cell sizes have been developed [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 map with multiplicative noise. They found that an adder mechanism can admit "blow-up" in which the expected cell sizes 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 generational cell 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 describe recently observed characteristics of population-level bacterial cell division. An adder-sizer model is one that incorporates a cell division rate β(x, y, t) and a single-cell growth rate g(x, y, t) that, instead of depending on a cell's age a, are functions of cell size x and a cell's volume added since birth y. Such an adder-sizer PDE model can be developed by defining n(x, y, t)dxdy as the mean number of cells with size in [x, x + dx] and added volume in [y, y + dy]. As cells have finite size and their added volume must be less than total size, n(x ≤ 0, y, t) = n(x, y ≥ x, t) = 0. A derivation similar to that given in [20] for the sizer model yields a transport equation 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 the right-hand-side of Eq. 2.1.
To explicitly outline our general derivation, consider the total population flux into and out of the size and added size domain Ω shown in Fig. 1(a) and defineβ(x ′ , y ′ , z ′ , t)dz ′ 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 ′ + dz ′ ] and the other with size within . This differential division function allows mother cells to divide into two daughter cells of differing sizes (asymmetric division), a process that has been observed in numerous contexts [14,13,2]. We also assume that daughter cells must have positive size soβ( 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 ′ > 0 but with size x ′ > x, we have two subcases. If the dividing cell has size x < x ′ < 2x, it will produce one daughter cell in Ω if a daughter cell has size 0 < z ′ < x ′ − x or x < z ′ < x ′ as shown in Fig. 1(b). If x ′ − x < z ′ < x, both daughter cells have size < x. Finally, if the dividing cell has size x ′ > 2x, at most one daughter will have size x ′ < x (see Fig. 1(b)). Upon simplifying the above birth terms by using x dz ′ for x ′ > x and the symmetryβ(x ′ , y ′ , z ′ , t) =β(x ′ , y ′ , x ′ − z ′ , t), we combine terms to balance proliferation with transport and find Fig. 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 in time dt, the domain Ω = [0, y] × [0, x] infinitesimally distorts Ω → Ω + dΩ through the growth increment gdt. The total population within this distorted domain changes only due to birth and death. Cells within Ω that divide always give rise to two daughters within Ω, leading to a net change of +1 cell. (b) The z ′ and x ′ domains of 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.
Upon taking the derivatives ∂ 2 ∂x∂y , we find the PDE given in Eq. 2.1 where the total division rate is defined by β(x, y, t) := x 0β (x, y, z, t)dz. For the boundary condition at y = 0, we take the derivative ∂/∂x and set y → 0 + to find The other boundary condition defined by construction is n(x, x, t) = 0.
The above derivation provides an explicit boundary condition representing newly born cells that may be asymmetric in birth size. Quantities such as the total cell population N (t) and the mean total biomass M (t) (the total volume over all cells) can be easily constructed from the density n(x, y, t): Higher moments of the total volume can also be analogously defined. By applying these operations to Eq. 2.1 and using the boundary condition (Eq. 2.3), we find the dynamics of the total population and biomass Finally, we also define the distribution of division events over the size and added size variables, accumulated over a time T : 2.1. Division probability and connection to time-sizer model. 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 "sizertimer" 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 0 , the probability that the cell splits within time [a, a + da] is defined by γ(a;ā)da. In the absence of death, to ensure that any single cell will eventually split, ∞ 0 γ(a;ā)da = 1. Reasonable choices for γ(a;ā) are Gamma, lognormal, or normal distributions. Without loss of generality, we propose a simple gamma distribution for γ(a;ā): whereā is the mean division age and σ 2 a is the variance. This type of distribution can be derived from the sum of independent, exponentially distributed ages.
For determinisitic exponential growth g = λx, age a and the parameterā can be explicitly expressed in terms of x, y and possibly other fixed parameters: in which ∆ is the fixed added size parameter that represents the adder mechanism. With a(x, y) andā(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);ā(x, y)): Assuming this "hazard function" form of a growth law, cells born at small initial size x(0) = x 0 = x − y take longer time to divide, while cells born with large size split sooner. Using the gamma distribution, we find a division rate of the form (2.11) β(x, y) = Γ ā 2 (x,y) σ 2 a γ(a(x, y);ā(x, y)) 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 product of the overall division rate β(x, y, t) and a differential division probability. The simplest model is to assume Fig. 2: The size and added-size dependent rate β(x, y) constructed using a gamma distribution for the splitting probability γ (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 deviation leads to a lower overall cell division rate β. When x is large,ā defined in 2.9 is small, a nonzero division rate β(x, y → 0) > 0 arises indicating that large newborn cells divide quickly to control size across the population. This particular feature arises from our construction of β as a hazard function. Modifying birth rate at small values of y so that β(x, y = 0) → 0 will not qualitatively change the predicted densities as long as the birth rate peak persists at small y.
that the differential division probability h(r) is a function of only the ratio r between the size of the daughter cell 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, 1]. The boundary condition (Eq. 2.3) can thus be written in the form A reasonable model for h(r = x/x ′ ) is a lognormal form that is symmetric about r = 1/2: where the parameters δ and σ r determine the bias and spread of the daughter cell size distribution, and the normalization constant is Z(σ r , δ) = 1 0 (h 0 (r) + h 0 (1 − r))dr. 2.2. Numerical Implementation and Monte-Carlo Simulations. With the differential birth rate functionβ defined, we can now consider the implementation of numerical solutions to Eqs. 2.1 and 2.3 as well as event-based simulations of the underlying corresponding stochastic process. Since a typical initial condition may not be smooth, a classical solution to Eqs. 2.1 and 2.3 may not exist. Thus, we provide a proof of existence and uniqueness of the weak solution to Eqs. 2.1 and 2.3 in Appendix A. We show convergence of a discrete approximation 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 in which both x and y are discretized with step size h. We define locally averaged functions by (2.15) f i+ 1 2 ,j+ 1 where f (x, y, t) can represent n(x, y, t), g(x, y, t), or β(x, y, t). Similarly, in the domain i, j ≥ 0 and j, k < i. The discretization of the transport equation can be expressed as 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+ 1 2 ,i = 0 to prevent density flux out of the y < x domain. In Eq. 2.17, g i+1,j+ 1 2 (t) can be taken as g( The discretized version of the boundary condition (Eq. 2.3) can be expressed as 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 numerically computed deterministic distributions (see Appendix C). We construct a list of cells and their associated sizes and 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 process results in, at any given time, a collection of cells, each with a specific size and added size. This collection of cells represents a realization of the population that should be approximated by the distributions that are solutions to Eqs. 2.1 and 2.3.

Analysis and Extensions.
In this section, we numerically investigate the adder-sizer model and plot various cell population densities and birth event distributions under different parameter regimes. We also show the consistency of numerical solutions of the adder-sizer PDE with results from direct Monte-Carlo simulations of the corresponding stochastic process, which demonstrates that numerical solutions of the linear PDE model for cell population is in agreement with single-cell level stochastic models. After investigating birth rate parameters that can lead to blow-up of population-averaged cell sizes, we extend the basic adder model to include motherdaughter growth rate correlations and processes that measure added size from different points in the cell cycle, i.e., an initation-adder model.

Cell and division event densities.
We evaluated our adder-sizer PDE model by using the division rate given in Eq. 2.10 and first assuming the simple and well-accepted growth function g(x, y, t) = λx. Fig. 3 shows the numerical results for the densityn(x, y, t) = n(x, y, t)/N (t) at successive times t = 1, 4, 12, respectively. Stochastic simulations of the underlying process yield cells populations consistent with the deterministic densities derived from the PDE model. In Fig. 4, we compare the cell densitiesn(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 the broader the cell and event densities.

Cell Volume
Explosion. At the single-cell level, a stochastic map model by Kessler and Burov assumed a multiplicative noise and predicted that cell sizes can eventually grow without bound, in agreement with what was experimentally observed for filamentous bacteria [17]. However, stochastic maps of generational cell size do not capture population-level distributions in size or age. In this subsection, we will numerically explore how a possible "blowup" in the population-averaged cell volumes. Within PDE models that describe population distributions, timer and sizer mechanisms have been shown to exhibit blow-up depending on properties of the birth rate β(a, x) [1,7,16]. Analysis of the conditions on full differential division rateβ(x, y, z, t) that would Numerically computed densitiesn(x, y, t) = n(x, y, t)/N (t) using g(x, y, t) = λx andβ(x, y, z, t) defined by Eqs. 2.10, 2.8, and 2.14. For all plots, we use σ a = 0.1 in γ(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.1, δ = 0) and plotn(x, y, 1),n(x, y, 4), andn(x, y, 12), respectively. In (d-f), we plot the densities using a broad (in fact, double-peaked) differential division function h(r) with parameters σ r = 0.2, δ = 0.7. In all calculations, we assumed an initial condition corresponding to a single newly born (y = 0) cell with size x = 1. For more asymmetric cell division in (d-f), the density spreads faster. In these cases, the densities closely approach a steady-state distribution by about t = 12. Also shown in each plot are realizations of Monte-Carlo simulations of the discrete process. Individual cells are represented by blue dots which accurately sample the normalized continuous densitiesn(x, y, t).
result in blow-up in the "adder-sizer" PDE model is more involved. Here, we provide only a heuristic argument for sufficient conditions for blow-up. First, we characterize the shape of the densities in the adder-sizer model. In the analogous McKendrick equation [15] one can investigate the age profile defined by dividing the number density by the total population size. The long term age profile may be stable even when the total population size continuously increases. We take a similar approach here by analyzingn(x, y, t) = n(x, y, t)/N (t) where N (t) is given by Eq. 2.5. Writing the adder-sizer PDE in terms ofn, we find Integrating this equation over x, y leads toṄ /N = ∞ 0 dx x 0 dy βn, which can be substituted into the first term in Eq. 3.1 to yield the nonlinear PDE A number of standard approaches may be applied to analyze Eq. 3.2. For example, in [15], solutions are attempted by controlling the analogous non-linear integral term. In the adder-sizer problem, we can define β(t) = Ω βn in the above expression to find a self-consistent condition on β(t) . One can also assess the steady-staten ss by setting ∂nss ∂t = 0 and establishing convergence. One indication of blow-up is a diverging mean cell size x(t) = 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 in which q(t) := Ω gn. If β, g, andn =n ss are time-independent and a steady state mean cell size exists, we expect it to obey x(∞) = q(∞)/ β(∞) . For the special case of deterministic exponential growth g(x) = λx, we can write the time evolution of the mean size as 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 the precise conditions on β leading to cell volume explosion are difficult to find, we will explore this possible phenomena using numerical experiments. We numerically examine the density n(x, y, t → ∞) and the mean cell size x(t) using the β,β defined in 2.10, 2.8, and 2.14.
In  . The curve associated with the σ a = 0.2 saturates while the one corresponding to σ a = 1 exhibits blow-up. However, the blowup is suppressed if a death term (µ = ln 2) is included.

Mother-daughter growth rate correlation.
Recent experiments indicate that the growth rate of a mother cell is "remembered" by its daughter cells. For growth rates of the form g(x, y, t) = λx, the exponential growth 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-averaged growth rate while in [7], changes in growth rates across two consecutive generations are modeled as a Markov process in order to estimate a division rate function β. In this subsection, we first introduce a generalized adder-sizer PDE incorporating variability in λ and then explore the mother-daughter growth rate correlation affects the population dynamics.
A mother-daughter growth rate correlation between two consecutive generations can be described by where ξ is a random variable, 0 ≤ R < 1 is the successive-generation growth rate correlation, andλ is the mean long-term, or preferred growth rate. Given a growth rate λ i of a mother cell, Eq. 3.5 describes the predicted growth rate λ i+1 of its daughter cells. We assume that the random variable has mean zero and is distributed according to some probability density P (ξ), which vanishes for ξ ≤ (1 − R)λ to ensure that the growth rates remain 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 growth rate λ. We propose the extended PDE model A possible symmetric mean zero distribution that vanishes at −(1 − R)λ takes on a log-normal form: If we start with one newly born daughter cell at size x 0 and growth rate λ 0 , the initial condition in our PDE model would be n 0 (x, y, λ) = δ(x − x 0 )δ(y)δ(λ − λ 0 ). 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 λ(t) converge to different values. For larger correlation 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 cells tend to grow slower. Because it takes shorter time for faster cells to divide, they will produce more generations of faster-growing cells, leading to a larger average growth rate defined as .
On the other hand, for a fixed mother growth rate λ i , smaller correlations R lead to mean daughter cell growth rates λ i+1 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 λ(t) to values smaller thanλ. Fig. 6(c) explicitly shows that when R = 0, the mean growth rate approaches a value smaller thanλ = ln 2.

Initiation-Adder Model.
Recent experiments suggest a new type of adder mechanism for bacterial cell size control [25]. Rather than a fixed volume added between birth and division as the primary control parameter, new experimental evidence suggests that the control parameter in E. coli is the added volume between successive initiations of DNA replication. Initiation occurs when the ori sites in a cell's genome are separated, 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 assumes that a cell's volume per initiation site (the ori site in the genome) tends to add a fixed volume between two consecutive initiations.
If the number of ori sites in a cell is q, initiation increases the number to 2q. Immediately after division and DNA separation, the number of oris 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 have initiated DNA replication but that have not yet divided. We define n 1 (x, y, t)dxdy as the expected number of pre-initiation cells in with volume in [x, x + dx] and with added volume y < x in [y, y + dy]. Mean post-initiation cell numbers with volume in [x, x + dx] and added volume in [y, y + dy] are described by n 2 (x, y, t)dxdy. In the general initiation-adder process, when a pre-initiation cell commences DNA replication (initiates) can depend on the volume or added volume. Thus, we describe transitions from a pre-initiation cell transitions into a post-initiation cell by the rate k i (x, y, t). After initiation, the number of ori sites doubles and the added volume is reset to zero in the newly formed post-initiation cell. In analogy with the differential division rate in Eq. 2.1,  Fig. 7: Schematic for the initiation adder process. DNA replication is initiated (indicated by the red dot) before copied DNA is segregated and cell division. In this example, q = 1 and y 2 is and added volume per origination site for two origination sites. The density of cells with q = 1 copy of DNA (before DNA replication initiation) is denoted n 1 (x, y, t) while the density of cells post-initiation is denoted n 2 (x, y, t), where y denotes the volume added after initiation. The factor that controls y 1 + y 2 in the initiation-adder model is the volume ∆ added between 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, we assume that the added volume is divided proportionally to the volume of the daughter cells, i.e., if the mother cell's volume is x with added volume y since initiation, and if one daughter cell's volume is z < x and the other 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 two coupled densities n 1 and n 2 : x z x n 2 (z, yz/x, t)β(z, x, yz/x, t)dz, ∂n 2 (x, y, t) ∂t + ∂[g 2 n 2 ] ∂x + ∂[g 2 n 2 ] ∂y = −β(x, y, t)n 2 , n 1 (x, 0, t) = 0, g 2 n 2 (x, 0, t) = x 0 k i (x, y, t)n 1 (x, y, t)dy, β(x, y, t) = x 0β (x, z, y, t)dz, in which we have allowed for different growth rates in the different cell phases. Both n 1 and n 2 are defined in the domain {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, but after initiation. In [30], a strong size control acting on initiation initiation was proposed where all cells will have inititated DNA replication before reaching some fixed volume x i . This hypothesis can be implemented in our initiationadder model by setting k i (x → x i , y, t) → ∞. The probability that a cell born at time t 0 has not yet initiated, ki(x(s),y(s),s)ds = 0 for all t 0 , then the largest volume that any cell can attain will be e λτ0 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(df). Note that due to cell size control affecting the pre-initiation stage, initial daughter cell sizes stay small at initiation and n 1 (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 after division. We have checked numerically that for constant k i = 10 3 , the densities n 1 (x, y, t) are negligible while n 2 (x, y, t) approaches the density of the division adder shown in Fig. 3 (for the same differential division functions  N (1, 0.1) and the sameβ(x, y, z, t) as that used in Fig. 3(d-f). (a-c) shows the normalized densitiesn 1 (x, y, t) ≡ n 1 (x, y, t)/N (t) where N (t) = dy dx(n 1 + n 2 ).
(d-f) shows the normalized post-initiation densityn 2 (x, y, t). For the k i used in this example, the pre-initiation densities span larger volume and added volumes. The densities are indistinguishable 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 can be seen from the first of Eqs. 3.9 where n 1 can be neglected and is dominated by the two terms on the righthand-side. Substituting k i (x, y, t)n 1 ≈ 2 ∞ x dz x n 2 (z, yz/x, t) into the integral terms in the second equation, we find Eq. 2.1 for n 2 (x, y, t).

Summary and Conclusions.
In this paper, we used PDE models to describe population dynamical behavior under the adder division mechanism. Under certain conditions, this PDE for the adder mechanism can also be converted to the well-known size-and age-structured PDE. In the absence of death, we motivated models for the differential birth rate functionβ(x, y, z, t) that are consistent with normalized division probabilities In Appendix 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 and space are unbounded as this problem is related to other first order PDE models that have been studied in more detail.
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 simulations of the discrete process generate sample configurations. The observed configurations are are consistent with 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, under the same initial conditions, a broader cell density n(x, y, t) and a broader event density ρ D (x, y, T ). We also demonstrate numerically, the divergence of the mean cell size x(t) = M (t)/N (t). We showed that division probabilities that are broader in the age or added size (and smaller in magnitude) more likely lead to mean cell sizes that explode with time.
We then incorporated growth rate correlation between cells of successive generations [17] into our "addersizer" PDE model. By extending the dimension of the density function to include growth rates and allowing for variability in growth rate as new cells are born, we developed a PDE model that incorporated the stochastic nature of growth rate inheritance and that describes evolution of the growth rate distribution of cells. We found that the steady-state value of the mean growth rate depends on the correlation of growth rates between mother and daughter cells. This dependence arises from a subtle interaction between the shape of the growth rate distribution 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 population dynamics under a new initiation-adder mechanism suggested by recent experimental results. In the limit that the initiation rate k i of DNA replication is significantly faster than all other time scales in the problem, the numerical solutions of the initiation adder model (Eq. 3.9) converge to those of division adder model (Eqs. 2.1  and 2.3). Under proper assumptions that come from experimental findings, we found that the initiation adder would also lead to effective cell size control [30].
There are new cellular processes and size control mechanisms that have been recently discovered and that can be mathematically modeled. Thus, there is likely general mathematical topics that remain to be explored within PDE and PIDE models of structured populations. For example, a recent experimental study indicates that an adder mechanism may be the result of several consecutive processes in the cell division cycle, suggesting that a much more complicated coupled system of PDEs/PIDEs would be required. tional Institutes of Health (R01HL146552) and the Army Research Office (W911NF-18-1-0345).
where the independent variables (x, y, t) ∈ R 2 ∩ {y < x} × R + . First, we assume that and nondimensionalize the size and added size by ∆, the added size parameter defined in Eq. 2.9. We also impose an additional assumption on g: We also assume the initial distribution n 0 (x, y) is compactly supported in (0, Ω) × [0, Ω), Ω < ∞. From this assumption and A.3, the closure of n(x, y, T )'s support is compact for any finite time T since n = 0 only when y < x and x(s) ≤ Ce Ks − (1 + T ) from Grönwall's Inequality, where C < 1 + T + Ω is given by the initial 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 support of n, we can assume the closure of g, β,β's support to be compact. One can generalize the definition of the weak solution n to [(R + ) 2 ∩ {y < x}] × [0, ∞) as in [22]. ∂Ψ ∂t + g(x, y, t) ∂Ψ ∂x + g(x, y, t) ∂Ψ ∂y − β(x, y, t)Ψ(x, y, t) holds for all test functions Ψ ∈ C 1 (([0, Ω(T )]) 2 ∩ {y ≤ x}) × [0, T ]) satisfying Ψ(x, y, T ) ≡ 0, Ψ(Ω(T ), y, t) = 0 and Ψ(x, x, t) = 0, where we set g,β, β = 0 for x ≥ Ω(T ), x ≤ y or x ≤ z. Upon using the boundary condition in A.1, the right-hand-side becomes Note that if n ∈ C 1 (((R + ) 2 ∩ {y < x}) × R + ) is a classical solution to the PDE (Eq. A.1), then it must also satisfy Eq. A.4 in any time interval [0, T ]. We refer to [22] for a proof of the existence and uniqueness of a weak solution of a related, simpler renewal equation. However, our adder-sizer PDE is more complicated. The proof of uniqueness requires very different techniques from the sizer PDE; yet the proof of existence is similar to 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). Taking the difference between using these purported solutions, we obtain Proof. We can transform the above equation into an ODE along the characteristic line and use contraction mapping, which is a standard practice in functional analysis to prove existence and uniqueness of the solution to 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 ), 0 ≤ y t < x t , x t − y t = c. SinceΨ(c, T ) = 0 andΨ(c, t) = Ψ(x t , x t − c, t), We can see that if x ≤ y or , s ≤ t ≤ T , and choose x s small enough such that x s < min{1, It is easy to verify that T is a contraction mapping for Ψ(x t , 0, t) and thus there exists a unique solution Ψ 0 satisfying A.6 in D 0 defined as D 0 = {(x, t)|s ≤ t ≤ T, 0 ≤ x ≤ x(x s , t)}. We then let x 1 s > x s and define s , t)} such that the difference of the area between regions D 1 and D 0 is less thanB −1 . Next, define a second mapping T 1 by T 1 is also a contraction mapping and we can obtain a Ψ 1 on D 1 such that T (Ψ 1 ) = Ψ 1 . Denote and it is easy to verify that Ψ is C 1 continuous on D 0 ∩ D 1 by first proving it is continuous and then taking the partial derivatives, and Ψ satisfy A.6 in the region D 0 ∪ D 1 . Following the same procedure, we can extend Ψ to satisfy A.6 in the region t ∈ [s, T ]. Then, for [0, s], we choose as close enough to s and use the same strategy by defining T 2 as (A.11) z, r) .
We finally obtain a unique function Ψ satisfying A.6 in [0, T ] × [0, ∞). From A.8, the value of Ψ is determined byS, Ψ(x, 0, t), U and we conclude that there exists a unique C 1 solution for A.6.
A.1.2. Uniqueness of weak solution for the adder-sizer model. From Section A.1.1 we obtain the existence and uniqueness of Ψ of the adjoint problem. Given any time T and S(x, y, t) ∈ C 1 (R + × (R + ) 2 ) satisfying the condition in Theorem A.1, since we can set g, β,β's support to be compact in [0, T ], we can find a unique C 1 continuous Ψ satisfying A.6. By substituting A.6 into A.5, we obtain (A.12) for any S(x, y, t) ∈ C 1 (R + × R + 2 ) satisfying S(x ≤ y, t) = S(x ≥ Ω(T ), y, t) = 0, which implies n ≡ 0 a.e. in 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 + ) 2 ∩{y < x}. In [22] such work is done for the renewal equation. We do not discuss this generalization in detail here.

A.2. Existence of the weak solution.
We construct a series of functions {n i } with a limit n for this series satisfying A.6 for all test functions Ψ. We use a semi-discrete approximation to discretize the PDE and obtain piecewise solutions. As the mesh size becomes smaller, we expect the piecewise solution to converge to a function n satisfying A.4. The idea of constructing a series of piecewise constant solutions and proving their convergence 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 > 0 fixed in both x and y axis and let time t be continuous. We denote Here, β i+ 1 2 ,j+ 1 . Given a fixed time T , we wish to find a solution of pointwise function n h (t), which takes values on the grid points (x i+ 1 2 , y j+ 1 2 ). Then n h can be seen as a vector function. According to our assumption there exists Ω such that the initial value n 0 is nonzero within the region {(x, y)|y < x, x < Ω}, and from our previous calculation there exists Ω(T ) < ∞ such that n is nonzero within the region {(x, y)|y < x, x < Ω(T )}. Eventually, we will set h(k) = Ω(T )/k and let the mesh size h → 0 by letting 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) By multiplying the second equation in A.14 by sign(n i+ 1 2 ,j+ 1 2 ) for each i, j ≤ i pair and summing over index We can simplify the above expression to We then have d dt Now, we estimate |n i+ 1 2 ,j+ 1 2 (t)| by first defining P (t) ≡ max 0≤i≤L−1,0≤j≤i−1 At a fixed time 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 maximum occurs at i − 1 > j > 0, we have while if the maximum occurs at j = i − 1 > 0, we have 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 minimumt < t such that P (v) > S(v)e −g ′ v for v ∈ (t, t]. Ift = 0, and since P (t) is nonincreasing from Eq. A.21, giving the second conclusion in Theorem A.2 that the L ∞ bound is uniform and independent of h. (t)1(ih ≤ x < (i + 1)h, jh ≤ y < (j + 1)h), where above, h = h(k) and 1 is the indicator function. Since there is an upper bound for both β andβ, and both β,β are continuous, we have the following result 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 1 ∩L ∞ ([0, Ω(T )] 2 ∩ {y < x} × [0, T )), so n h(k) are all L 2 functions. There exists a function n ∈ L 2 ([0, Ω(T )] 2 ∩ {y < x} × [0, T )) and a subsequence k i → ∞ that satisfies n h(ki) ⇀ n. Since L 2 ([0, Ω(T )] 2 ∩ {y < x} × [0, T )) implies L 1 integrability, we can deduce that n is an L 1 function as desired.
To prove n h(ki) ⇀ n, we need only to verify that there exists a subsequence n h(ki) such that for all test functions f ∈ L 2 , which gives the result n h(ki) ⇀ n. We can show that n is a weak solution by multiplying the first two of Eqs. A.14 by a test function Ψ ∈ C 1 ([0, Ω(T )] 2 × [0, T )), Ψ(x, y, T ) = 0, Ψ(x, y, t) = 0, y ≥ x for which For a given L ∈ N + and h = Ω(T ) L , ds ∂Ψ ∂s (s, y, t), |n h | is uniformly bounded while g is C 1 continuous. From above we can pick a subsequence in {n h(k) }, denoted by n h(ki) ⇀ n. We use n h = n h(ki) in the above formula. Since Ψ ∈ C 1 [0, T ] × [0, Ω(T )] 2 , given any Ψ we have a positive upper bound R(Ψ) < ∞ for Ψ and any of its first derivatives. Thus, ∂Ψ(x, y, t) ∂t . As ∂t | → 0 since ∂Ψ ∂t and n h(ki) are all bounded. Moreover, so that the first term in Eq. A.26 tends to the limit in Eq. A.27. By the same procedure and using the condition that g is uniformly continuous in [0, T ] × [0, Ω(t)] 2 (g is C 1 ), it is easy to verify that the second and third terms on the LHS of Eq. A.26 tend to It is also easy to verify that the first and second terms on the RHS of A.26 tend to − x 0 dy β(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 ODE system A.14 it is an approximation to the solution of A.14. Provided h, ∆t → 0 satisfies the CFL condition 2 g ∞ ∆t < h, conclude that at least a subsequence of the numerical solutions converge to the unique weak solution of A.1. Furthermore, recently, the existence to an eigenpair of the adder-sizer PDE A.1 under specific smooth conditions satisfied by the coefficients g, β,β has been proved in [10], allowing for studying asymptotic behavior of the solution.
Appendix B. Numerical Scheme. We denote u(t) = {n 1 (t), n 2 (t), . . . , n L−1 (t)} T where n j (t) = in which b i is a lower bidiagonal matrix with diagonal and d i is a diagonal matrix diag(d i ) = g j+ 1 2 ,i−1 (t) dt h , j = i, i + 1, . . . , L − 1. 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 1 (x i , y i , t, b 1 ), ..., 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, 1).
Suppose we have a β 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 dx dt = λx. We require that G i = ai 0 γ(a ′ )da ′ < b i at the beginning of each step for every i. • Step 2: For each cell, calculate G i = ai+∆t 0 γ(a ′ )da ′ . If G i ≥ b i , then we numerical calculate a ∆t i such that ai+∆ti 0 γ(a ′ )da ′ ≈ 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 and let 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, go to step 5. • Step 4: Remove cell i from S(t ′ ), record its volume x at t ′ , and generate the random numbers r from the distribution h(r) and b m , b m+1 from U(0, 1). Then, add two new cells in S(t ′ ) labeled by c m (rx, 0, t, b m ) and c m+1 (x − rx, 0, 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 automatically satisfied at t = 0. For our runs, we used 10 cells of initial volume 0.5 and t max = T is the same as the maximum time for the numerical PDE experiments. We can also generalize the model to incorporate the mother-daughter growth coefficient correlation by including a new label λ i to each cell.