PDE models of adder mechanisms in cellular proliferation
Mingtao Xia Chris D. Greenman Tom Chou
q-bio.CB math.AP 35Q92, 92C37









PDE models of adder mechanisms in cellular proliferation
    Mingtao XiaDepartment of Mathematics, UCLA, Los Angeles, CA, USA 
Chris D. GreenmanSchool of Computing Sciences, University of East Anglia, Norwich, UK
Tom ChouDepts. of Computational Medicine and Mathematics, UCLA, Los Angeles, CA, USA
    version of December 30, 2023
================================================================================================================================================================================================================================================

equationsection


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.
  




PDE, structured populations, cell size control




 35Q80, 92B05, 92C37




§ 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
<cit.>.  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 <cit.>, or added volume since birth y
<cit.>. 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 <cit.>.

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 <cit.>. The McKendrick
“transport” equation for the cell density takes the form
∂_tn(a,t) + ∂_an(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∫_0^aβ(s)n(s,t) s describes the birth of zero-age
cells. Fully demographically stochastic versions of the timer model
have also been recently developed
<cit.>.

The timer (or age-dependent) model does not explicitly track cell
sizes, but PDE models incorporating sizer mechanisms have been
developed <cit.>.  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 <cit.>.  Existence and uniqueness of
weak solutions to timer and sizer models have been proved for certain
boundary and initial conditions.  These types of structured population
equations can be partially solved using the method of characteristics
but the boundary conditions can only be reduced to a Volterra-type
integral equation <cit.>.

Much like a general growth law g(a,x,y,t) that can depend on age,
size, added size, and time, the three distinct mechanisms of cell
division need not be mutually exclusive. In this paper, we mainly
focus, at the cell 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, has been 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 <cit.> and can be
motivated by an initiator accumulation mechanism distinct from those
used to justify sizers or timers <cit.>.

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 <cit.>. 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
<cit.>. 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 the latest
“initiation adder” mechanism, which states that the added mechanism
takes effect on the cell's size at initiation instead of division in
<cit.>. 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
<cit.>. These models describe how cell sizes
evolve with generation and can interpolate among timer, sizer, and
adder mechanisms.  Kessler and Burov <cit.> 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. <cit.> 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.




§ 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) x  y as the mean number of cells with size
in [x,x+ x] and added volume in [y, y+ y].  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 <cit.> for the sizer model yields a transport equation of
the form


    [ ∂ 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. <ref>.

To explicitly outline our general derivation, consider the total
population flux into and out of the size and added size domain
Ω shown in Fig. <ref>(a) and define
β̃(x',y',z',t)  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'+ z'] and the other with size within [x'-z',
  x'-(z'+ 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 of differing
sizes (asymmetric division), a process that has been observed in
numerous contexts <cit.>. We also assume
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'>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. <ref>(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. <ref>(b)).  Upon simplifying the above birth terms by using
∫_0^x' z' = ∫_0^x z'+ ∫_x^x' z'
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


    ∫_0^x x' ∫_0^y y'
     ∂ n(x',y',t)∂ t
    +    ∫_0^x x'  g(x',y,t)n(x',y,t)  + ∫_0^y y' 
    g(x,y',t) n(x,y',t) 
         = ∫_0^∞ y' 
    ∫_0^x x' ∫_0^x' z'
     β̃(x',y'z',t)n(x',y't) 
         + ∫_y^∞ y' ∫_0^x x'
    ∫_0^x' z' β̃(x',y,z',t)n(x',y',t) 
         + 2∫_0^∞ y' ∫_x^∞ x'
    ∫_0^x z' β̃(x',y',z',t)n(x',y',t).


Upon taking the derivatives ∂^2∂ x ∂
  y, we find the PDE given in Eq. <ref> where the total
division rate is defined by β(x,y,t) ∫_0^xβ̃(x,y,z,t) z. For the boundary condition
at y=0, we take the derivative ∂/∂ x and set y→
0^+ to find


    g(x,y=0,t)n(x,y=0,t) = 2∫_x^∞ x'
    ∫_0^x' 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/2),
and boundary condition of the adder-sizer model reduces to 


    g(x,y=0,t)n(x, y=0,t)= 4∫_0^2xβ(2x, y', t)n(2x, y', t)  y'.


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):


    N(t)  =∫_0^∞ x∫_0^x y  n(x, y,t),   
    M(t)  =∫_0^∞ x ∫_0^x y  x n(x, y, t).


Higher moments of the total volume can also be analogously defined.  By
applying these operations to Eq. <ref> and using the boundary
condition (Eq. <ref>), we find the dynamics of the total
population and biomass


    N(t)/ t = ∫_0^∞ x ∫_0^x y β(x, y, t)n(x, y, t),   M(t)/ t = 
    ∫_0^∞ x ∫_0^x 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, accumulated over a time T:


    ρ_ d(x,y,T)  = ∫_0^Tβ(x,y,t)n(x,y,t) t/∫_0^T t ∫_0^∞ x'∫_0^x y' β(x',y',t)n(x',y',t).





 §.§ 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
“sizer-timer” structured population model <cit.>.  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+ a] is defined by γ(a; a̅) a. In
the absence of death, to ensure that any single cell will eventually
split, ∫_0^∞γ(a; a̅) a = 1.  Reasonable
choices for γ(a; a̅) are Gamma, lognormal, or normal
distributions. Without loss of generality, we propose a simple gamma
distribution for γ(a; a̅):


    γ(a;a̅) =
    1/aΓ((a̅/σ_a)^2)exp[-aa̅σ_a^2+(a̅σ_a)^2ln(aa̅σ_a^2)],


where a̅ is the mean division age and σ_a^2 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 a̅ can be explicitly expressed in terms of x,y and
possibly other fixed parameters:

    a(x,y) = 1/λln(x x-y),   a̅(x, y) = 1λln(x-y+Δ/x-y),


in which Δ is the fixed added size parameter that represents
the adder mechanism.

With a(x,y) and a̅(x,y) defined in Eqs. <ref>, 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)):


    β(x,y,t) = γ(a(x,y); a̅(x, y))/1-∫_0^a(x,y) a' γ(a'; a̅(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


    β(x,y) = Γ(a̅^2(x,y)σ_a^2)γ(a(x,y); a̅(x, y))Γ(a̅^2(x,y)σ_a^2,
               a(x,y)a̅(x,y)σ_a^2),


where Γ(·, ·) is the upper incomplete gamma function.
We plot two examples of the time-independent rate β(x,y) in
Fig. <ref>.



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


    β̃(x,y,z,t) = β(x,y,t) h(z/x)/x,


where r ≡ z/x ∈[0,1].  The boundary condition
(Eq. <ref>) can thus be written in the form


    g(x,0,t)n(x,0,t) =2 ∫_x^∞ x' ∫_0^1
    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/2:


    h(r)  = h_0(r) + h_0(1-r) Z(σ_r,δ),    
    h_0(r)  = e^-(-δ + ln
        r)^2/2σ_r^2e^-ln^2(1-r)/2σ_r^2,


where the parameters δ and σ_r determine the bias and spread of
the daughter cell size distribution,  and the normalization constant is
Z(σ_r,δ)= ∫_0^1 (h_0(r)+h_0(1-r))  r.




 §.§ Numerical Implementation and Monte-Carlo Simulations

With the differential birth rate function β̃ defined, we
can now consider the implementation of numerical solutions to
Eqs. <ref> and <ref> 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. <ref> and <ref> may not exist.  Thus, we provide a
proof of existence and uniqueness of the weak solution to
Eqs. <ref> and <ref> in Appendix <ref>.  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


    f_i+1/2,j+1/21/h^2∫_ih^(i+1)h x ∫_jh^(j+1)h y  f(x,y,t),


where f(x,y,t) can represent n(x,y,t), g(x,y,t), or
β(x,y,t).  Similarly,


    β̃_i+1/2,j+1/2((s+1/2)h, t) 
    = h^-3∫_ih^(i+1)h x
    ∫_jh^(j+1)h y ∫_kh^(k+1)h z β̃(x,y,z,t)


in the domain i, j ≥ 0 and j,k < i.  The discretization of the
transport equation can be expressed as


    [ n_i+1/2,j+1/2(t+) - 
    n_i+1/2,j+1/2(t)/ + 
    g_i+1,j+1/2ñ_i+1,j+1/2-g_i,jñ_i,j+1/2/h
    + g_i+1/2,j+1ñ_i+1/2,j+1
    -g_i+1/2,jñ_i+1/2,j/h;                                                                                               = -β_i+1/2,j+1/2 n_i+1/2,j+1/2(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+1/2,i=0 to prevent density
  flux out of the y<x domain.  In Eq. <ref>, g_i+1,
  j+1/2(t) can be taken as g((i+1)h, (j+1/2)h, t)
while ñ_i+1, j+1/2(t) = ∫_jh^(j+1)hy  n((i+1/2)h, y, t) is a finite-volume numerical
approximation to ∫_jh^(j+1)hy  n((i+1)h, y,
t). The discretized version of the boundary condition
(Eq. <ref>) can be expressed as


    g_i+1/2,0n_i+1/2,0(t)
    =2h^2∑_k=i+1^L∑_j=0^k-1β̃_k+1/2,j+1/2((i+1/2)h, t)
    n_k+1/2,j+1/2(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 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 . 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. <ref> and <ref>.



§ 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 mother-daughter 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. <ref> and first assuming the simple and
well-accepted growth function g(x,y,t) = λ x. Fig. <ref>
shows the numerical results for the density n̅(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. <ref>, 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 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 <cit.>.

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)
<cit.>. Analysis of the
conditions on full differential division rate β̃(x,y,z,t)
that would 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 <cit.> 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 analyzing n̅(x,y,t) = n(x,y,t)/N(t)
where N(t) is given by Eq. <ref>. Writing the adder-sizer PDE in
terms of n̅, we find







    ∂n̅/∂ t + n̅/N N/ t + 
    ∂ (gn̅)/∂ x+∂(gn̅)/∂ y 
    = -βn̅.


Integrating this equation over x,y leads to Ṅ/N =
∫_0^∞ x ∫_0^x y  βn̅,
which can be substituted into the first term in Eq. <ref> to
yield the nonlinear PDE


    ∂n̅/∂ t + ∂ (gn̅)∂ x
    + ∂ (gn̅)∂ y = -(β + 
    ∫_Ωβn̅)n̅.


A number of standard approaches may be applied to analyze
Eq. <ref>. For example, in <cit.>,
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-state n̅_
  ss by setting ∂n̅_ ss/∂ 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. <ref> by x
and integrating (using the boundary condition and symmetry of the
β̃ distribution) we find

    ⟨ x(t) ⟩ t +
    ⟨β(t)⟩⟨ x(t)⟩ = q(t),


in which q(t)∫_Ω gn̅.  If β, g, and
n̅ = 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


    ⟨ x(t)⟩/ t=
    [λ - ⟨β(t)⟩] ⟨ x(t)⟩,   ⟨β(t)⟩≡∫_0^∞ x∫_0^x 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. <ref>), 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
<ref>, <ref>, and <ref>.



In Fig. <ref>(a) and (b) we plot the marginal distribution
n̅(x,t)∫_x^∞ y  n(x,y,t)/
∫_0^∞ x ∫_x^∞ y  n(x,y,t) for
different values of the division rate variability σ_a at
different times. The associated division rates correspond to those
plotted in Fig. <ref>(a) and (b).  In Fig. <ref>(c)
we plot the mean cell sizes ⟨ x(t)⟩ = M(t)/N(t)
corresponding to the distributions in (a) and (b).  For sufficiently
broad division probabilities γ(a) (large σ_a), the
division rates β are small, and ⟨ x(t)⟩ fails to
saturate and diverges.



 §.§ 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
<cit.>.  In <cit.>,
fluctuations in λ have been discussed at the single-cell level
to explore their effects on the population-averaged growth rate while
in <cit.>, 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


    λ_i+1 = (λ_i-λ̅)R + λ̅ + ξ,


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


    {∂ n(x,y,t,λ)/∂ t +
    ∂ (gn)/∂ x   
    + ∂ (gn)/∂ y
    =-β(x, y,t)n(x, y, t, λ), 
    
    g(x,0,t,λ)n(x, 0, t, λ)=     2∫_0^∞λ' ∫_x^∞ x'
    ∫_0^x' 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, 0, λ) =      n_0(x, y, λ),
    .





A possible symmetric mean zero distribution that vanishes at
-(1-R)λ̅ takes on a log-normal form:


    P(ξ) ∝exp[-ln^2(ξ+(1-R)λ̅)/2σ_ξ^2
      -ln^2((1-R)λ̅ -ξ)/2σ_ξ^2].



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. <ref> shown in
Fig. <ref> 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


    ⟨λ(t) ⟩ = ∫_0^∞x∫_0^xy∫_0^∞λ  λ n(t, x, y, λ)/∫_0^∞x∫_0^xy∫_0^∞λ n(t, x, y, λ).


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. <ref>(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 <cit.>.  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)xy as the expected number of
pre-initiation cells in with volume in [x,x+x] and with added
volume y<x in [y, y+y].  Mean post-initiation cell numbers
with volume in [x,x+x] and added volume in [y, y+y] are
described by n_2(x, y, t)xy. 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. <ref>, 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:


    ∂ n_1(x,y,t)/∂ t 
    + ∂ [g_1 n_1]/∂ x + 
    ∂ [g_1n_1]/∂ y     
    = -k_ i(x,y,t)n_1 + 2∫_x^∞z/xn_2(z,yz/x,t)
    β̃(z, x, yz/x,t)z, 
    ∂ n_2(x, y, t)/∂ t + 
    ∂ [g_2n_2]/∂ x+ ∂ [g_2n_2]/∂ y     
    = -β(x,y,t)n_2, 
    
    n_1(x,0,t) = 0,    g_2n_2(x, 0, t)     
    = ∫_0^x k_ i(x,y,t)n_1(x, y, t)y, 
    β(x, y, t)     = ∫_0^xβ̃(x, z, y, t) z,


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
{ℝ^+^2∩{y<x}}×ℝ^+.  These coupled
PDEs are different from the PDE associated with the standard
“division-adder” described in Eqs. <ref> and <ref>.
Here, the added volume is reset to zero not after division, but after
initiation.

In <cit.>, 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 initiation-adder model by setting
k_ i(x→ x_ i, y, t)→∞. The probability that a
cell born at time t_0 has not yet initiated, e^-∫_t_0^t
  k_ i(x(s), y(s), s)s, always vanishes for all (t_0,
x_t_0, y_t_0) before some finite time t and x(t)< x_ i.
Thus, n_2(x, 0, t) is nonzero only in [0, x_ i] for all t.
If there exists a constant τ_0 such that
lim_τ→τ_0 e^-∫_t_0^t_0+τk_
    i(x(s), y(s), s)s=0 for all t_0, then the largest volume
that any cell can attain will be e^λτ_0x_ i,
leading to strict size control and no blowup.



Fig. <ref> shows numerical solutions to Eq. <ref> using the
same birth rate function as that used in Fig. <ref>(d-f). 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. <ref> (for the
same differential division functions β̃). Thus, the
initiation adder model converges to the standard division adder model
when k_ i→∞.  This can be seen from the first of
Eqs. <ref> where n_1 can be neglected and is dominated by
the two terms on the right-hand-side. Substituting k_
  i(x,y,t)n_1≈ 2∫_x^∞ z/xn_2(z,yz/x,t)
into the integral terms in the second equation, we find
Eq. <ref> 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 <cit.> into our “adder-sizer” 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. <ref>) converge to those of division adder model
(Eqs. <ref> and <ref>). Under proper assumptions that
come from experimental findings, we found that the initiation adder
would also lead to effective cell size control <cit.>.

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.




PDE-type models can be used to model cell densities that
  evolve according to timer, sizer, or adder mechanisms, as well as
  combinations of mechanisms such as the sizer-timer model and the
  “adder-sizer” model studied here. Under a deterministic cell
  growth assumption, one might propose a growth rate function
  g(a,x,y) and birth rate β(a,x,y) that depend on all three
  variables, age a, size x, and added size y.  Thus, one might
  propose a full sizer-timer-adder model of the form
  (∂_t+∂_a)n + ∂_x(gn)+∂_y(gn) =
  -β(a,x,y,t)n(a,x,y,t) with purported boundary condition
  g(0,x,0,t)n(0,x,0,t) = 2∫_0^∞
  a'∫_x^∞ x' ∫_0^x' y' β̃(a',x',y',x,t)n(a',x',y',t).  However, the three
  variables ,a,x,y are not all independent.  For example, if the
  deterministic added size y(t) is monotonic in time t, the age
  after birth a and the added size y are functions of each
  other. More generally, if we can determine the evolution of all
  three variables (a,x,y) given two of them, we cannot construct a
  meaningful 3+1-dimensional PDE model.  One can understand the loss
  of independence by noticing that when a cell divides, both its
  daughter cells' ages and their added volumes reset to a=y=0. For
  deterministic growth, the age and added size are bijective and are
  thus not independent.  Subsequent deterministic growth of the
  daughter cells is described by their sizes and either their
  age a or added size y.

Thus, the age and added size variables are not independent
  and given timer, sizer, adder mechanisms, there are only two
  structurally different PDE models, the sizer-timer PDE, and the
 “adder-sizer" PDE studied in detail here.  The sizer-timer model can
  be reduced to just a sizer model or a timer model (McKendrick
  equation). For example, assuming age-independent division and growth
  rates β(x,t) and g(x,t), integrating Eq. <ref> over
  ages a, and using Eq. <ref>, we find ∂_tn(x,t) +
  ∂_x(gn) = -β(x,t)n(x,t) +
  2∫_0^∞β(x',x,t)n(x',t) x' and the boundary
  condition n(x=0,t) = 0, while the "adder-sizer" model can be
  reduced to a simple adder model. Thus, assuming deterministic cell
  growth, one can consider only the timer, sizer, adder, sizer-timer,
  or "adder-sizer" PDE models.  An adder-only model can be defined
  only when both the division growth rates depend only on added size
  y and not on size x.  However, if the growth is itself
  stochastic, one might propose higher order models that can include
  all types of cell division mechanisms.


Within the PDE model, the boundary condition is nonzero in only one
dimension at certain time t_0, we can denote the nonzero line
n(0, x, 0, t_0) as l(t_0). As t increases l(t_0) will result
in a two dimensional nonzero space in the four dimensional space,
denoted by S(t_0). And take different t_0 will eventually result
in a three dimensional nonzero hyperplane based on the boundary
condition. We can see that eventually the nonzero space in the four
dimensional space (a,x,y,t) will collapse to a three
dimensional hyperplane, indicating one of the variables is redundant.

Also, some other extensions in this topic include derive similar PDE
models for other cell division mechanisms like the timer and sizer
models, and then prove for the natures of the solutions for those
models. And, given specific division rate, it will interesting to use
different models, namely for cells of different mechanisms to simulate
and observe the difference between numerical solutions' behavior of
those PDE models and compare them with experimental data. And
conversely, given initial condition and all the parameters, we may be
able to tell which division strategy the cell takes by implementing
different PDE models and observe which one's solution is most similar
to the actual data. From our comparison of those three models, we can
already distinguish the adder mechanism by numerically solving the
adder model PDE and capture the long time behavior of the solution
since it will result in the volume explosion of cells, yet for the
other two cells' size appears to be bounded, which is in
correspondence to the probability result in <cit.>. So we could
expect these PDE models can successfully predict the total population
change, as well as the volume distribution and other features given
initial conditions.





§ ACKNOWLEDGMENTS

This work was supported in part by grants from the NSF (DMS-1814364),
the National Institutes of Health (R01HL146552) and the Army Research
Office (W911NF-18-1-0345).  
 











§ 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. The full problem is defined as 


    {∂n/∂t + ∂(ng)/∂x
    + ∂(ng)/∂y =-β(x, y, t) n(x, y, t),
    
    g(x,0,t)n(x, 0, t)=2∫_x^∞ x'∫_0^x'y β̃(x', y, x, t)n(x', y, t),
    β(x, y, t) ∫_0^xβ̃(x, y, z, t) z,
    β̃(x, y, z', t) = β̃(x, y, z-z', t),  β̃(x, y, 0, t) = 0,   n(x, x, t) =0,
    
    n(x, y, t=0)  n_0(x, y).
    .


where the independent variables  (x, y,t)∈ℝ^2∩{y<x}×ℝ^+.

First, we assume that


    0 < g_min≤g∈   C^1({(ℝ^+)^2∩{y≤x}}×ℝ^+), 
    
    n_0(x, y) ∈   L^1∩L^∞∩C^1(ℝ^+∩{y<x}), 
    
    0 ≤β̃∈   L^∞∩L^1∩C^1({(ℝ^+)^3∩{y<x, z<x}}×ℝ^+), 
    β(x, y, t) ∈   L^∞∩L^1∩C^1({(ℝ^+)^2∩{y≤x}}×ℝ^+),

and nondimensionalize the size and added size by Δ,
the added size parameter defined in Eq. <ref>.
We also impose an additional assumption on g:


    |g(x, y, t)| < K(t + x + 1),    K < ∞.


We also assume the initial distribution n_0(x,y) is compactly
supported in (0, Ω)×[0, Ω), Ω<∞.  From this
  assumption and <ref>, 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 [(ℝ^+)^2∩{y<x}]× [0,
      ∞) as in <cit.>.

Definition A.1 Given time T<∞ and assuming <ref>, for
a function n∈L^1((([0,
  Ω(T)])^2∩{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-sizer PDE in the weak sense in time [0, T],
    if 

    -∫_0^T t∫_0^∞ x∫_0^x y  n(x,
    y, t)[∂Ψ/∂t+g(x,y,t)
    ∂Ψ/∂x+g(x,y,t)∂Ψ/∂y
    -β(x,y,t)Ψ(x,y,t)] 
    
    =∫_0^∞ x∫_0^x y  n_0(x, y)Ψ_0(x, y) +
    ∫_0^T t∫_0^∞ x Ψ(x, 0, t)n(x, 0, t)g(x, 0, 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
<ref>, the right-hand-side becomes


    ∫_0^∞ x∫_0^x y  n_0(x, y)Ψ_0(x, y) +
    2∫_0^T t ∫_0^∞ x∫_0^x y∫_0^x z Ψ(z,0,t)β̃(x,y,z,t)n(x,y,t).


Note that if
n∈C^1(((ℝ^+)^2∩{y<x})×ℝ^+)
is a classical solution to the PDE (Eq. <ref>), then it must
also satisfy Eq. <ref> in any time interval [0, T].  We refer
to <cit.> 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 <cit.>.



 §.§ Uniqueness

First, we prove uniqueness of the solution to <ref>. Assume
there are two weak solutions n^(0) and n^(1) for the adder-sizer PDE
satisfying <ref> with the same initial condition
n_0^(0)(x, y)=n_0^(1)(x, y). Taking the difference between using these
purported solutions, we obtain 


    -∫_0^T t∫_0^∞ x∫_0^x y Δ n(x,y,t)
    [∂Ψ/∂t+g(x,y,t)
    ∂Ψ/∂x+g(x,y,t)
    ∂Ψ/∂y-β(x,y,t)Ψ(x,y,t)]
    
    =2∫_0^T t ∫_0^∞ x ∫_0^x y ∫_0^x z Ψ(z, 0, t)β̃(x,y,z,t)Δ n(x,y,t),

where Δ n = n^(1) - n^(0). 




  §.§.§ Adjoint Problem

We consider the adjoint problem for Ψ in the given time interval
[0, T] and with a source term S(x,y,t):


    ∂Ψ/∂t+g(x,y,t)
    ∂Ψ/∂x+g(x,y,t)
    ∂Ψ/∂x-β(x,y,t)Ψ(x,y,t) 
       = -2∫_0^xΨ(z, 0, t)β̃(x, y, z, t) z - S(x,y,t),   
    0≤y < x 
    Ψ(x,y,T)    =0,   Ψ(Ω(T), y, t)=0,   Ψ(x,x,t)=0.


Theorem A.1 Assume <ref>, and S∈C^1([0, Ω(T)]^2×[0,
  T]), S(Ω(T), y, t)=0, and S=0 when
x≤y. Then there exists a unique C^1 solution to the
adjoint problem.

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. <ref>, 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 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

where 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 above
equation becomes an ordinary differential equation. Given x_t, we
define


    {Ψ̃(c, s)Ψ(X(c, s), s)e^-∫_t^sβ(X(c, v), v) dv,                                        
     
    U(c, z, s) 2β̃(X(c, s), z, s)
    e^-∫_t^sβ(X(c, v), v) dv,   S̃(c, s) S(X(c, s), s)e^-∫_t^sβ(X(c, v), v) dv.
    .

Thus, along the characteristic line we can write <ref> as

    ∂/∂sΨ̃(c,s)
    =-∫_0^x(c,s)Ψ(z, 0, s)U(c, z, s) z-S̃(c,s).


Since Ψ̃(c, T)=0 and Ψ̃(c, t)=Ψ(x_t,x_t-c, t),


    Ψ(x_t, x_t-c, t)=∫_t^TS̃(c, s) s +∫_t^T s
    ∫_0^x(c, s) z Ψ(z, 0, s)U(c, z, s),    0 < 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≤0 or
x_t>Ω(T). Using c=x_t, Eq. <ref> becomes


    Ψ(x_t, 0, t)=∫_t^TS̃(x_t, s) s + ∫_t^T s
    ∫_0^x(x_t, s) z Ψ(z, 0, s)U(x_t,z,s).


From condition <ref> we obtain
x(s)≤(x_t+1+T)e^K(s-t)-(1+T). From condition A.3, we define
B̃ = 2β̃_∞ < ∞.
Next, we choose s=max{T-1/Kln(1+1/2B̃(1+T)),
T-1/Kln2, T-1} such that
e^K(T-t)≤1+1/2B̃(1+T), s≤t≤T, and choose x_s
small enough such that x_s<min{1, 1/8B̃(T-s)}. We
denote a mapping T defined on the functional space as


    T(Ψ)(x_t, 0, t) = ∫_t^TS̃(x_t, s)s+∫_t^T s ∫_0^x(s,
      x_t) z Ψ(z, 0, 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,
0, t) and thus there exists a unique solution Ψ_0 satisfying
<ref> in D_0 defined as D_0={(x, t)|s≤t≤T,
0≤x≤x( x_s, t)}. We then let x_s^1>x_s and define D_1=
{(x, t)|s≤t≤T, 0≤x≤x(x_s^1, t)} such that the
difference of the area between regions D_1 and D_0 is less than
B̃^-1. Next, define a second mapping T_1 by


    {
    T_1(Ψ)(x_t, 0, t) = ∫_t^T s ∫_x(x_s, s)^x(x_t, s) z Ψ(
    z, 0, s) U(x_t, z, s) + I(x_s, t),     t∈[s, T],    x_t∈[x(t,
      x_s), x_s^1],
    
    I(x_s, t)=∫_t^T s S̃(x_t, s) +∫_t^T s 
    ∫_0^x(x_s, s) z  Ψ_0(z,0,s) U(x_t, z, s).                          .


T_1 is also a contraction mapping and we can obtain a Ψ_1 on
D_1 such that T(Ψ_1)=Ψ_1. Denote



    Ψ(x, 0, t)={Ψ_0(x, 0, t),  (x, t)∈D_0,
    Ψ_1(x, 0, t),  (x, t)∈D_1,
    .

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 <ref> in the region
D_0∪D_1.

Following the same procedure, we can extend Ψ to satisfy
<ref> in the region t∈[s, T]. Then, for [0, s], we
choose a s̃ close enough to s and use the same strategy by
defining T_2 as


    {
    T_2(Ψ)(x_t, 0, t)= ∫_t^s r S̃(x_t, r)+
    ∫_t^s r ∫_0^x(x_t, r) z Ψ(z, 0, r) U(x_t, z, r)
    +Ĩ(t, x_s),    t∈[s̃, s], 
    Ĩ(x_s, t) = ∫_s^T r S̃(x_t, r)
    +∫_s^T r ∫_0^x(x_t, r) z Ψ(z, 0, r) U(x_t, z, r)                      .
    .

We finally obtain a unique function Ψ satisfying <ref> in [0,
  T]×[0, ∞).

From <ref>, the value of Ψ is determined
by S̃, Ψ(x, 0, t), U and we conclude that there exists a 
unique C^1 solution for <ref>.




  §.§.§ 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(ℝ^+×(ℝ^+)^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 <ref>.
By substituting <ref> into <ref>, we 
obtain


    ∫_0^T t ∫_0^Ω(T) x ∫_0^x y Δ n(x,y,t)S(x,y,t) = 0

for any S(x,y,t)∈C^1(ℝ^+×ℝ^+^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,
∞)×(ℝ^+)^2∩{y<x}. In <cit.> such
work is done for the renewal equation. We do not discuss this
generalization in detail here.



 §.§ Existence of the weak solution

We construct a series of functions {n_i} with a limit n for this
series satisfying <ref> 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
<ref>.  The idea of constructing a series of piecewise constant
solutions and proving their convergence to a weak solution is similar
to that in <cit.>.



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


    (x_i, y_j)=(ih, jh), (x_i+1/2,  y_j+1/2)=
    ((i+1/2)h, (j+1/2)h),     j< i ∈N, 
    β_i+1/2, j+1/2(t) = 1/h^2∫_ih^(i+1)h y 
    ∫_jh^(j+1)h x β(x, y, t),    j<i ∈N,
    β̃_i+1/2, j+1/2((s+1/2)h, t) =
    1/h^3∫_ih^(i+1)h z ∫_jh^(j+1)h y 
    ∫_sh^(s+1)h x β̃(x,y,z,t),    s≤i,
    
    g_i, j(t)=g(ih, jh, t),    j<i ∈1/2N.


Here, β_i+1/2,
  j+1/2(t)=h∑_s=0^iβ̃_i+1/2,
  j+1/2((s+1/2)h,t). 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. <ref>, 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):


    h^2 n_i+1/2, j+1/2(t)/ t+h(g_i+1, j+1/2(t)
    n_i+1/2, j+1/2(t)-g_i, j+1/2(t)n_i-1/2, j+1/2(t))  
         + h(g_i+1/2, j+1(t)n_i+1/2, j+1/2(t)
    -g_i+1/2, j(t)n_i+1/2, j-1/2(t)) + h^2β_i+1/2, j+1/2(t)
    n_i+1/2, j+1/2(t)=0,    0≤j<i-1 
         h^2 n_i+1/2, j+1/2(t)/ t 
    + hg_i+1, j+1/2(t)n_i+1, j+1/2(t)  
         -hg_i+1/2, j(t)n_i+1/2, j-1/2(t) 
    + h^2β_i+1/2, j+1/2(t)
    n_i+1/2, j+1/2(t)=0,    0≤j=i-1 
         g_i+1/2, 0(t) n_i+1/2,
      -1/2(t) = 2h^2∑_ℓ=i^L-1∑_j=0^ℓ-1β̃_ℓ+1/2, j+1/2
    ((i+1/2)h,t) n_ℓ+1/2,
      j+1/2(t),  
         n_i+1/2,
      j+1/2(0)=1/h^2∫_x_i^x_i+1 y
    ∫_y_j^y_j+1 x  n_0(x, y),    n_i+1/2, i+1/2(t)=0,


where we henceforth omit the h superscript in the proof.
In the two-dimensional upwind scheme, derivatives in one direction are
neglected on neighboring sites in the other direction:
n_i,j±1/2=n_i-1/2, j±1/2,
n_i±1/2, j=n_i±1/2, j-1/2.  The
boundary condition n(x, x, t)=0 is implemented by n_i+1/2,
  i+1/2(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^1 continuous,
which means that there exists a unique solution in time [0, T],
T<∞.

Theorem A.2 For t∈[0, T] and assuming <ref> holds, 
we find the bound

    ∑_i=1^L-1∑_j=0^i|n_i+1/2, j+
    1/2(t)|≤e^M t∑_i=1^L-1∑_j=0^i|n_i+1/2, j+1/2(0)|,


where B̃=2β̃_∞, M=2B-b,
B=β_∞, and b=min_tmin_i,
  jβ_i+1/2, j+1/2(t). The
L^∞ bound is given by n^h(t)_∞≤max{1/g_ minB̃e^MTn(0)_1,
n^h(0)_∞}e^2g̃'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 <ref> by
sign(n_i+1/2, j+1/2) for each i, j≤ i,


    h^2/ t|n_i+1/2, j+1/2(t)|+
    hg_i+1, j+1/2(t)|n_i+1/2, j+1/2(t)|+
    hg_i+1/2, j+1(t)|n_i+1/2, j+1/2(t)|+
    h^2β_i+1/2, j+1/2(t)|n_i+1/2, j+1/2(t)|≤
    
    hg_i, j+1/2(t)|n_i-1/2, j+1/2(t)|+
    hg_i+1/2, j(t)|n_i+1/2, j-1/2(t)|


By multiplying the second equation in <ref> by
sign(n_i+1/2, j+1/2) for each i, j≤ i pair
and summing over index ∑_i=1^L-1∑_j=0^i-1,


    h^2∑_i=1^L-1∑_j=0^i-1|n_i+1/2,
      j+1/2(t)|
    
    
    
    
    +h∑_j=0^i-1g_L,
      j+1/2(t)|n_L-1+1/2,
      j+1/2(t)|
     + h^2∑_i=1^L-1∑_j=0^i-1β_i+1/2,
      j+1/2(t)|n_i+1/2, j+1/2(t)|
    ≤h∑_i=0^L-1g_i+1/2, 0(t)|n_i+1/2,
      -1/2(t)|.



We can simplify the above expression to


    h^2d/dt∑_i=1^L-1∑_j=0^i-1| n_i+1/2, j+1/2(t)|    +     h^2∑_i=1^L-1∑_j=0^i-1β_i+1/2, j+1/2| n_i+1/2, j+1/2(t)|
        ≤    2h^3∑_i=0^L-1| ∑_ℓ=i^L-1∑_j=0^ℓ-1β̃_ℓ+1/2, j+1/2((i+1/2)h,t)
    n_ℓ+1/2, j+1/2(t)| 
        ≤     2h^2∑_ℓ=1^L-1∑_j=0^ℓ-1
    |β_ℓ+1/2, j+1/2(t)||n_ℓ+1/2, j+1/2(t)|.


We then have


    / t∑_i=1^L-1∑_j=0^i-1
    |n_i+1/2, j+1/2(t)|≤(2B-b)
    ∑_i=1^L-1∑_j=0^i-1|n_i+1/2, j+1/2(t)|,

which yields

    ∑_i=1^L-1∑_j=0^i-1|n_i+1/2, j+1/2(t)|
    ≤e^Mt∑_i=1^L-1∑_j=0^i-1
    |n_i+1/2, j+1/2(0)|.


<ref> states that the l^1 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_1≤i≤L-1
|n_i+1/2,1/2(t)|e^-g̃'t for t∈[0,
  T]. For the maximum value of S at some index i, we find


    [ h^2|n_i+1/2, 1/2(t)|/ t+h(g_i+1,
      1/2(t)|n_i+1/2, 1/2(t)|-g_i,
      1/2(t)|n_i-1/2, 1/2(t)|) + h(g_i+1/2, 1(t)|n_i+1/2,
      1/2(t)|-g_i+1/2, 0(t)|n_i+1/2,
      -1/2(t)|)≤ 0,;                                                           h^2|n_i+1/2,
        1/2(t)|/ t+hg_i+1, 1/2(t)|n_i+1/2,
      1/2(t)|-g_i+1/2, 0(t)|n_i+1/2,
      -1/2(t)| ≤ 0,    i=1. ]


and


    (|n_i+1/2,
        1/2(t)|e^-g̃'t)/ t+h^-1g_i+1/2,
        1(t)|n_i+1/2,
      1/2(t)|e^-g̃'t≤ h^-1g_i+1/2,
        0(t)|n_i+1/2, -1/2(t)|e^-g̃'t,


By the assumption that g(x, y,t)≥g_ min(t)≥g_ min>0 and
g< K (T+1+Ω(T)), we have


    (|n_i+1/2,
        1/2(t)|e^-g̃'t)/ t+ h^-1g_ min(t) |n_i+1/2,
      1/2(t)|e^-g̃'t≤ h^-1(g_ min(t)/g_ min)max_1≤i≤L-1|g_i+1/2,
      0(t)n_i+1/2, -1/2(t)|.









Finally, defining G(t)=h^-1∫_0^t g_ min(s) s yields


    (|n_i+1/2,
        1/2(t)|e^-g̃'te^G(t))/ t≤1 h(g_ min(t)/g_ min)
    max_1≤i≤L-1
    |g_i+1/2,0(t)n_i+1/2, -1/2(t)|e^G(t).

From the L^1 bound, we can deduce


    max_tmax_1≤i≤L-1|g_i+1/2(t)n_i+1/2,
      -1/2(t)|≤h^2B̃e^MTn^h(0)_1
    ≤B̃e^MTn(0)_1,   
    t>0


and conclude that for the function S(t)e^G(t)

    S(t)e^G(t)≤ S(0) + 1/g_ minB̃e^MTn(0)_1(e^G(t)-1),

and S(t)≤max_1≤i≤L-1{n_i+1/2,1/2(0),1/g_
    minB̃e^MTn(0)_1}, which then gives the
L^∞ bound for the pointwise solution n^h when
j=0.


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{|n_i+1/2,
    j+1/2(t)|e^-2g̃'t}. 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











    h/ t(|n_i+1/2,
      j+1/2(t)|e^-2g̃'t) = - [
    g_i,j+1/2(t)|n_i+1/2, j+1/2(t)| - g_i+1,
        j+1/2(t)|n_i+1/2, j+1/2(t)|. 
    . +g_i+1/2, j(t)|n_i+1/2, j+1/2(t)| -
      g_i+1/2, j+1(t)|n_i+1/2,
        j+1/2(t)| + 2hg̃'|n_i+1/2,j+1/2(t)|] 
    e^-2g̃'t≤ 0,


while if the maximum occurs at j=i-1>0, we have


    / t(|n_i+1/2,
      j+1/2(t)|e^-2g̃'t)≤[h^-1(g_i+1/2,
        j(t) - g_i+1, j+1/2(t))|n_i+1/2,
      1/2(t)|-2g̃'|n_i+1/2,
      j+1/2(t)|] e^-2g̃'t≤ 0.


In Eqs. <ref> and <ref>, 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. <ref>,
P(t)≤ P(0)= n^h(0)_∞. If t>t̃ >0,
P(t)≤P(t̃)≤S(t̃)≤max_0≤t≤TS(t),
while if t̃=t,
P(t)=S(t)≤max_0≤t≤TS(t). Thus, 
P(t) =  n^h(t)_∞
e^-2g̃'t≤max{max_0≤t≤T{S(t)},
n^h(0)_∞} and


    n^h(t)_∞≤max{max_0≤t≤T{S(t)},
    n^h(0)_∞}e^2g̃'t,


giving the second conclusion in Theorem A.2 that the
  L^∞ bound is uniform and independent of h.



  §.§.§ Existence of the weak solution



For a given time T<∞, we can take the grid size
h(k)=Ω(T)/k → 0 by letting the integer k →∞. Spatially piecewise constant functions can then be defined
based on the sequence of vector functions {n^h(k)}.  By setting
n^h_i+1/2, i+1/2(t)=0, we define n^h(x,y,t),
β^h, and β̃^h as


    n^h(x, y, t) =    ∑_i=0^k-1∑_j=0^i-1n^h_i+1/2, j+1/2(t)
    1(ih≤x<(i+1)h, jh≤y<(j+1)h), 
    β^h(x, y, t)=    ∑_i=0^k-1∑_j=0^iβ_i+1/2, j+1/2(t)1(ih≤x<(i+1)h, jh≤y<(j+1)h), 
    β̃^h(x, y, z, t) =    ∑_i=0^k-1∑_j=0^i-1∑_ℓ=0^i-1β̃_i+1/2, j+1/2((ℓ +1/2)h, t)
    1(ih≤x<(i+1)h, jh≤y<(j+1)h, ℓ h≤z<(ℓ +1)h), 
    
    n^h(x, 0, t)=      n^h_i+1/2, -1/2(t),    ih≤ x<(i+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


    lim_k→∞β^h(k)(x, y, t)
    →   β(x, y, t)   a.e.    0≤β^h(k)≤β_∞<∞,
    lim_k→∞β̃^h(k)(x, y, z, t)
    →   β(x, y, z, t)   a.e.    
    0≤β̃^h(k)≤β̃_∞<∞,
    lim_k→∞n^h(k)(x, y, 0)→   n(x, y, 0)  a.e..


Then, we can apply Theorem A.2 to the 
piecewise constant solutions n^h(k) of Eqs. <ref>.

Corollary A.3 Under the conditions of Theorem A.2,
for any t∈[0, T] and any h,


    ∫_0^Ω(T) y ∫_0^Ω(T) x  |n^h(x, y,
    t)| ≤e^Mt∫_0^Ω(0) y∫_0^Ω(0)
    x  |n^h(x, y, 0)|


and 


    n^h(t)_∞≤max{| n(0)|_∞,
       Be^MT| n(0)|_1}e^2g̃'̃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^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(k_i)⇀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(k_i)⇀n, we need only to verify that
there exists a subsequence n^h(k_i) such that for all test
functions f∈L^2, ∫_0^T
t∫_0^Ω(T) x ∫_0^x y   n^h(k_i)f
→∫_0^T t∫_0^Ω(T) x∫_0^x
y  nf. Since L^2 space is separable, we have a countable
set of basis function {b_i(x, y, t)} for the space
L^2([0,Ω(T)]^2∩{y<x}× [0,T)).  Thus,
  every n^h(k) can be decomposed as n^h(k) =
  ∑_i=1^∞α^k_ib_i.  The sequence {n^h(k)} is
  uniformly L^∞ bounded, so ∑α_k^2 are all
  uniformly bounded. We can then select a subsequence {n^h(k_i)}
  from {n^h(k)} satisfying
  lim_i→∞α^k_i_j=α_j so that
  ∑_i=1^∞α_j^2 <∞. If we set
  n=∑_i=1^∞α_i b_i, then, by decomposing any test
  function Ψ∈L^2([0,Ω(T)]^2∩{y<x}×
  [0,T)) by Ψ = ∑_i=1^∞γ_i b_i, we have


    lim_i→∞| ∫_0^T t ∫_0^Ω(T) x 
    ∫_0^x y  (n^h(k_i)-n)
    Ψ|
    = | ∑_s=1^∞(α^k_i_s-α_s)γ_s|=0,

which gives the result n^h(k_i)⇀n.  

We can show that n is a weak solution by multiplying the
  first two of Eqs. <ref> by a test function
  Ψ∈C^1([0,Ω(T)]^2× [0,T)),
    Ψ(x,y,T)=0, Ψ(x,y,t)=0, y≥x for which


    Ψ_i+1/2,
      j+1/2(t)≡1/h^2∫_x_i^x_i+1
    x∫_y_j^y_j+1 y  Ψ(x,y,t),    j≤i.


For a given L∈N^+ and h=Ω(T)/L,


    ∫_0^T t ∑_i=1^L-1∑_j=0^i-1(
    h^2 n^h_i+1/2,
        j+1/2(t)/ tΨ_i+1/2, j+1/2(t) +
    h[g_i+1, j+1/2(t)n^h_i+1/2, j+1/2(t)-g_i,j+
    1/2(t)n^h_i-1/2,
        j+1/2]Ψ_i+1/2, j+1/2(t) 
    
    + h[g_i+1/2, j+1(t)n^h_i+1/2,
        j+1/2(t)-g_i+1/2, j(t)n^h_i+1/2,
        j-1/2]Ψ_i+1/2, j+1/2(t)+
    h^2β_i+1/2, j+1/2(t)n^h_i+1/2,
      j+1/2Ψ_i+1/2, j+1/2(t)) 
    
    
    = ∫_0^T t ∑_i=1^L-1h
    g_i+1/2, i(t)n^h_i+1/2, i-1/2(t)Ψ_i+12, i-12(t), 
        n^h_i+1/2, i+1/2 = 0.

Integrating the above equation by parts with respect to time, we find


    ∫_0^T t [
    ∑_i=1^L-1∑_j=0^i-1 h^2n^h_i+1/2,
      j+1/2(t)Ψ_i+1/2,
        j+1/2(t)/t + h∑_i=1^L-2∑_j=0^i-1 
    g_i+1,j+1/2(t)n^h_i+1/2,j+1/2(t)
    (Ψ_i+3/2, j+1/2(t)-Ψ_i+1/2, j+1/2(t)). 
    . + h∑_i=1^L-1∑_j=0^i-2 g_i+1/2,
      j+1(t)n^h_i+1/2, j+1/2(t)(Ψ_i+1/2,
      j+3/2(t)-Ψ_i+1/2, j+1/2(t))] 
     =
    
    -h^2∑_i=0^L-1∑_j=0^i-1n^h_i+1/2,
      j+1/2(0)Ψ_i+1/2, j+1/2(0) 
    - h∫_0^T t ∑_i=1^L-1g_i+1/2,
      0(t)n^h_i+1/2, -1/2(t)Ψ_i+1/2,
      1/2(t)
    
    -∫_0^T t ∑_i=1^L-1h
    g_i+1/2, i(t)n^h_i+1/2, i-1/2(t)Ψ_i+12, i-12(t)
    
    +h∫_0^T t [∑_j=0^L-2g_L, j+1/2(t)n^h_L-1/2,
      j+1/2(t)Ψ_L-1/2,
      j+1/2(t) + ∑_i=1^L-1∑_j=0^i-1h^2β_i+1/2,
      j+1/2(t)n^h_i+1/2,
      j+1/2(t)Ψ_i+1/2, j+1/2(t)] .


Since Ψ_i+3/2, j+1/2(t)-Ψ_i+1/2,
  j+1/2(t) = ∫_ih^(i+1)h
x∫_jh^(j+1)h y ∫_x^x+h s ∂Ψ/∂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(k_i)⇀n. We use
n^h=n^h(k_i) 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,


    |∫_0^T t ∑_i=1^L-1∑_j=0^i-1(h^2n^h(k_i)_i+1/2,
    j+1/2(t)Ψ_i+1/2,
    j+1/2(t)/ t) -
    ∫_0^T t ∫_0^Ω(T) x ∫_0^x y  
    n^h(k_i)(x, y, t)∂Ψ(x,y,t)/∂t| 
    ≤
    ∫_0^T t ∑_i=0^L-1∫_ih^(i+1)h x ∫_ih^x y |n^h(k_i)(
    x, y, t)∂Ψ(x, y, t)/∂t|.

As h→ 0, |∫_0^T
t∑_i=0^L-1∫_ih^(i+1)h x ∫_ih^x y 
n^h(k_i)(x,y, t)∂Ψ(x, y, t)/∂t| → 0
since ∂Ψ/∂t and n^h(k_i) are all
bounded. Moreover,


    ∫_0^T t ∫_0^Ω(T) x∫_0^x y  
    h^2n^h(k_i)(x,y,t)∂Ψ(x,y,t)/∂t→∫_0^T t∫_0^Ω(T) x∫_0^x y  
    h^2n(x, y, t)∂Ψ(x, y, t)/∂t


so that the first term in Eq. <ref> tends to the limit in
Eq. <ref>. 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. <ref> tend to ∫_0^T t
∫_0^Ω(T) x∫_0^x y 
(gn)(x,y,t)∂Ψ/∂x and ∫_0^T
t∫_0^Ω(T) x∫_0^x y  (gn)(x,
y,t)∂Ψ/∂y, respectively.

It is also easy to verify that the first and second terms 
on the RHS of <ref>





tend to 
 
-∫_0^T t∫_0^Ω(T)
x∫_0^x y  n(x,y,0)Ψ(x,y,0) and -2∫_0^T t
∫_0^∞ x∫_0^x y∫_0^x z Ψ(z,0,t)β̃(x,y,z,t)n(x,y,t), respectively.  The third
term on the RHS of <ref> h∫_0^T t
∑_i=1^L-1 g_i+1/2,i(t)
n^h_i+1/2,i-1/2(t)Ψ_i+1/2,i-1/2(t)
tends to 0 since Ψ is C^1 continuous and is 0 on the
boundary x=y. Since Ψ is continuous and is 0 at x=Ω(T),

 
h∫_0^T∑_j=0^L-2g_L,
  j+1/2(t)n^h_L-1/2,j+1/2(t)
Ψ_L-1/2,j+1/2(t) t → 0 as
h→ 0. Finally, the last term on the RHS of of <ref>
h∫_0^T∑_i=1^L-1∑_j=0^i-1h^2β_i+1/2,
  j+1/2(t)n^h_i+1/2,
  j+1/2(t)Ψ_i+1/2, j+1/2(t)t
→∫_0^T t ∫_0^Ω(T)
x∫_0^x 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
<ref>. Since the numerical solution obtained by the scheme in
Appendix B is a discretization in time for the ODE system
<ref> it is an approximation to the solution of
<ref>. Provided h,Δ t → 0 
satisfies the CFL condition 2g_∞Δ t < h,
conclude that at least a subsequence of the numerical
solutions converge to the unique weak solution of
<ref>.










Furthermore, recently, the existence to an eigenpair of the
adder-sizer PDE <ref> under specific smooth conditions
satisfied by the coefficients g, β, β̃ has been
proved in <cit.>, allowing for studying asymptotic
behavior of the solution.



§ NUMERICAL SCHEME


We denote (̆t) = {_1(t), _2(t), …,
_L-1(t)}^ where _j(t) =
{n_1/2,j-1/2, n_1+1/2,j-1/2,
…, n_L-1/2,j-1/2} and n_i≤ j = 0.
Equations <ref> and <ref> can then be written in the
form (̆t+Δ t) =  A(t)(̆t), where


    A(t) = [[ B_1+ C_1      C_2      C_3      C_4        ⋯    C_L-2    C_L-1;      D_2      B_2        0        0        ⋯        0        0;        0      D_3      B_3        0        ⋯        0        0;        ⋮        ⋮        ⋮        ⋮        ⋮        ⋮        ⋮;        0        0        0        0        ⋯    B_L-2        0;        0        0        0        0        ⋯    D_L-1    B_L-1;          ]],


is made up of the following L-1 L× L matrices



    B_i=
    [
        [                ;                ;                ;                ;         b_i    ;                 ]],   C_i=
    [
        [                ;                ;                ;                ;         c_i    ;                 ]],





    D_i=
    [
        [                ;                ;                ;                ;         d_i    ;                 ]],


in which b_i is a lower bidiagonal matrix with diagonal 


    diag( b_i)  = 1-1 hg_j+1,i-1 2(t) t
      -1 hg_j+1/2,i(t) t - 
    β_j+1 2,i-1 2(t) t,    j=i,i+1,…, L-1,


and lower off diagonal ( b_i)_-1 = g_j,i-1
    2(t) t  h,   j=i+1,…, L-1,







    ( c_i)_sj = {[ β̃_i-12+j,i-12
    ((s+12)h, t)t,           i+j-s-1>0,  i+j ≤ L;                             0                    otherwise. ].






    c_i =[[      β̃_i+1 2, i-1 2(3 2h, t)t      β̃_i+3 2, i-1 2(3 2h, t)t                              ⋯      β̃_L-1 2, i-1 2(3 2h, t)t;                              0      β̃_i+3 2, i-1 2(5 2h, t)t                              ⋯      β̃_L-1 2, i-1 2(5 2h, t)t;                              0                              0                              ⋯      β̃_L-1 2, i-1 2(7 2h, t)t;                              ⋮                              ⋮                              ⋮                              ⋮;                              0                              0                              ⋯ β̃_L-1 2, i-1 2((L- 3 2)h, t)t;                              0                              0                              ⋯                              0 ]],




and d_i is a diagonal matrix
diag( d_i) = g_j+1 2, i-1(t)  t  h,   
  j=i, i+1, …, L-1.



§ 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 <ref> and
β̃ of the form <ref>. 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=∫_0^a_iγ(a')a' < b_i at the beginning of each
  step for every i.

  * Step 2: For each cell, calculate
  G_i=∫_0^a_i+Δtγ(a')a'. If G_i≥ b_i,
  then we numerical calculate a Δt_i such that
  ∫_0^a_i+Δt_iγ(a')a'≈b_i.

  * Step 3: Choose the smallest Δt_i among all possible
  Δt_is 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.


siamplain



