Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

The impact of maturation time distributions on the structure and growth of cellular populations

  • Here we study how the structure and growth of a cellular population vary with the distribution of maturation times from each stage. We consider two cell cycle stages. The first represents early G1. The second includes late G1, S, G2, and mitosis. Passage between the two reflects passage of an important cell cycle checkpoint known as the restriction point. We model the population as a system of partial differential equations. After establishing the existence of solutions, we characterize the maturation rates and derive the steady-state age and stage distributions as well as the asymptotic growth rates for models with exponential and inverse Gaussian maturation time distributions. We find that the stable age and stage distributions, transient dynamics, and asymptotic growth rates are substantially different for these two maturation models. We conclude that researchers modeling cellular populations should take care when choosing a maturation time distribution, as the population growth rate and stage structure can be heavily impacted by this choice. Furthermore, differences in the models' transient dynamics constitute testable predictions that can help further our understanding of the fundamental process of cellular proliferation. We hope that our numerical methods and programs will provide a scaffold for future research on cellular proliferation.

    Citation: Asma Alshehri, John Ford, Rachel Leander. The impact of maturation time distributions on the structure and growth of cellular populations[J]. Mathematical Biosciences and Engineering, 2020, 17(2): 1855-1888. doi: 10.3934/mbe.2020098

    Related Papers:

    [1] Yijun Lou, Bei Sun . Stage duration distributions and intraspecific competition: a review of continuous stage-structured models. Mathematical Biosciences and Engineering, 2022, 19(8): 7543-7569. doi: 10.3934/mbe.2022355
    [2] Jordi Ripoll, Jordi Font . Numerical approach to an age-structured Lotka-Volterra model. Mathematical Biosciences and Engineering, 2023, 20(9): 15603-15622. doi: 10.3934/mbe.2023696
    [3] Tyler Cassidy, Morgan Craig, Antony R. Humphries . Equivalences between age structured models and state dependent distributed delay differential equations. Mathematical Biosciences and Engineering, 2019, 16(5): 5419-5450. doi: 10.3934/mbe.2019270
    [4] Ming Mei, Yau Shu Wong . Novel stability results for traveling wavefronts in an age-structured reaction-diffusion equation. Mathematical Biosciences and Engineering, 2009, 6(4): 743-752. doi: 10.3934/mbe.2009.6.743
    [5] Frédérique Billy, Jean Clairambault, Franck Delaunay, Céline Feillet, Natalia Robert . Age-structured cell population model to study the influence of growth factors on cell cycle dynamics. Mathematical Biosciences and Engineering, 2013, 10(1): 1-17. doi: 10.3934/mbe.2013.10.1
    [6] Eminugroho Ratna Sari, Fajar Adi-Kusumo, Lina Aryati . Mathematical analysis of a SIPC age-structured model of cervical cancer. Mathematical Biosciences and Engineering, 2022, 19(6): 6013-6039. doi: 10.3934/mbe.2022281
    [7] Qiuyi Su, Jianhong Wu . Impact of variability of reproductive ageing and rate on childhood infectious disease prevention and control: insights from stage-structured population models. Mathematical Biosciences and Engineering, 2020, 17(6): 7671-7691. doi: 10.3934/mbe.2020390
    [8] Edoardo Beretta, Dimitri Breda . Discrete or distributed delay? Effects on stability of population growth. Mathematical Biosciences and Engineering, 2016, 13(1): 19-41. doi: 10.3934/mbe.2016.13.19
    [9] Anthony Tongen, María Zubillaga, Jorge E. Rabinovich . A two-sex matrix population model to represent harem structure. Mathematical Biosciences and Engineering, 2016, 13(5): 1077-1092. doi: 10.3934/mbe.2016031
    [10] Hongying Shu, Wanxiao Xu, Zenghui Hao . Global dynamics of an immunosuppressive infection model with stage structure. Mathematical Biosciences and Engineering, 2020, 17(3): 2082-2102. doi: 10.3934/mbe.2020111
  • Here we study how the structure and growth of a cellular population vary with the distribution of maturation times from each stage. We consider two cell cycle stages. The first represents early G1. The second includes late G1, S, G2, and mitosis. Passage between the two reflects passage of an important cell cycle checkpoint known as the restriction point. We model the population as a system of partial differential equations. After establishing the existence of solutions, we characterize the maturation rates and derive the steady-state age and stage distributions as well as the asymptotic growth rates for models with exponential and inverse Gaussian maturation time distributions. We find that the stable age and stage distributions, transient dynamics, and asymptotic growth rates are substantially different for these two maturation models. We conclude that researchers modeling cellular populations should take care when choosing a maturation time distribution, as the population growth rate and stage structure can be heavily impacted by this choice. Furthermore, differences in the models' transient dynamics constitute testable predictions that can help further our understanding of the fundamental process of cellular proliferation. We hope that our numerical methods and programs will provide a scaffold for future research on cellular proliferation.


    In this paper we investigate how the distribution of time spent in the various cell cycle stages impacts the growth and stage structure of a cellular population. For brevity, we will refer to the distribution of time spent in a stage as the maturation time distribution of the stage. In addition to providing fundamental insights into the process of cellular proliferation, this research can inform diverse models of cellular population growth. Indeed, in the context of medical research, a cellular population's stage structure can impact the efficacy of drug therapy [1,2,3].

    Within the mammalian cell cycle there exists a checkpoint, known as the restriction point or G1/S checkpoint, which controls entry into S phase [4,5,6,7]. This checkpoint is regulated by growth factor signaling. As such, mammalian cells can coarsely be divided into those that have and those that have not received sufficient growth signals (mitogenic signals) to begin the process of cellular division [4,5,6,7]. Hence, we consider a model with two cell cycle stages, representing the stage prior to restriction point passage (early G1) and the stage delineated by restriction point passage and mitosis (late G1 through M) [6]. As mentioned above, we are especially interested in how the maturation time distributions impact the stage structure of the population, i.e., the fractions of cells that have and have not passed the restriction point.

    This work builds on previous research [1,8] which considered how division time distributions impact the age or generation structure of a cellular population. For example, [8] investigates the sensitivity of the generation structure to the distribution of division times (i.e., the intermitotic time distribution). While similar in many respects, an important difference between the model system considered by [8] and our own model system is that while generation number increases indefinitely, cell cycle stages occur in a cycle. As a result, the two models have different boundary conditions, which may considerably impact their dynamics. Another important difference between this work and [8] is that we consider general maturation time distributions and more general initial age distributions, which yield models that are not analytically solvable. A second paper of interest is [1] which develops methods for incorporating age dependency into models of cellular populations and demonstrates the utility of this approach for the study of drug therapy. Our work extends [1] by modeling, in addition to age, cell cycle stage. As a result of this extension, our model could be used to investigate both age- and stage-dependent effects in relation to drug therapy. Indeed, cell cycle dynamics and the stage structure of a cellular population are thought to be important for drug therapy [2]. Finally, this work also relates to [2,3] where stage-structured models are used to study the impact of drug therapy on pancreatic cancer cells. These models assume that, in the absence of treatment or crowding, maturation between cell cycle phases is governed by an exponential distribution, i.e., that cells experience a constant per capita maturation rate. In contrast, we investigate the impact of non-constant maturation rates on the stage structure of a cellular population. In summary, to the best of our knowledge, this is the first paper to consider the impact of maturation time distributions on the stage structure of a cellular population.

    In addition to providing theoretical results relating stable age and stage structures to maturation time distributions, we perform extensive numerical analyses. Our analyses compare and contrast the dynamics and stable stage- and age-distributions of two distinct maturation time distributions: the exponential and the inverse Gaussian distribution. Over several decades, numerous model maturation time distributions have been proposed [9,10,11,12,13,14,15,16,17,18,19,20]. Still, the most appropriate maturation time distribution remains unknown and is likely context dependent. We focus on the exponential maturation time distribution because it is almost certainly the most frequently used to model progression between cell cycle stages. This is a practical choice. Because the exponential distribution assumes that a cell has a constant probability per unit time of maturing into the next cell cycle stage, the maturation rate for this model is constant. As a result, cell age is an extraneous variable, in that the population stage structure can be modeled as a system of ordinary differential equations in time alone. For this reason, the exponential model may be considered an ageless model. A drawback of the exponential model is that the assumption of a constant transition rate is generally considered to be biologically unreasonable. Moreover, the exponential distribution provides a poor fit to data (the distribution mode is zero, whereas few, if any, cells are observed to mature from a cell cycle stage at age zero). As a reasonable alternative to the exponential distribution, we consider the inverse Gaussian distribution. This distribution describes the maturation time as a first exit time [21], so that the time to exit a cell cycle stage corresponds to the time it takes a random variable to achieve a threshold value. The corresponding maturation rate is thus age-dependent, and the distribution mode can be shifted away from zero. Moreover, this model has been demonstrated to provide an accurate approximation of the intermitotic time distributions of several cell lines [20]. In considering these two very different distributions, we hope to clarify how the qualitative behavior of the population varies with the choice of maturation time distribution.

    Section 2 of this paper introduces the system of partial differential equations (PDEs) used to model the stage-structured cellular population and presents results on the existence of solutions via the method of characteristics. Section 3 derives the steady-state growth rate of the population and presents the stable stage and age structure of the population. Section 4 presents the numerical algorithm for simulating the model's solution and addresses the challenges of simulating the model with age-dependent maturation rates. Section 5 presents the results of numerical simulations and compares and contrasts the growth and structure of populations governed by exponential and inverse Gaussian maturation rates.

    The PDE population model is presented below.

    gt(a,t)+ga(a,t)=βg(a)g(a,t);for a0, t0 (2.1)
    ft(a,t)+fa(a,t)=βf(a)f(a,t);for a0, t0 (2.2)

    with boundary conditions:

    g(a,0)=g0(a)fora0, (2.3)
    g(0,t)=20βf(a)f(a,t)dafort0, (2.4)
    f(a,0)=f0(a)fora0, (2.5)
    f(0,t)=0βg(a)g(a,t)da.fort0, (2.6)

    Here g gives the density of cells in the first stage, f gives the density of cells in the second stage, a denotes the ''age" of a cell relative to the time it entered its current stage, and t denotes time. According to (2.1) the rate of change of the density of cells in the first stage (gt) is determined by continuous aging (ga) and the age-dependent, per capita rate at which cells enter the second stage (βg(a)). Similarly, (2.2) describes how the age density of cells in the second stage evolves through time. Boundary conditions (2.3) and (2.5) give the initial age density at the beginning of the experiment (it is assumed that this density is known). Finally, (2.3) and (2.5) describe how new cells (with a=0) enter the first and second stages. Cells in the first stage are ''born" when cells in the second stage divide, thereby giving rise to two cells in the first stage, with the age-dependent, per capita rate of βf(a). Cells in the second stage are created when cells in the first stage mature with an age-dependent, per capita rate of βg(a).

    In this section we outline the method of proof of the existence of solutions to model equations (2.1)–(2.6) and present the related theorems. Detailed proofs are presented in the appendix.

    We employ an iterative method, in which an approximating sequence is shown to converge to a solution. This method of proof is similar to that from [25], where global existence was shown for a size-structured model with a single stage and bounded size. We denote the terms of the approximating sequence by gn and fn. These functions are defined as solutions of the following system of partial differential equations:

    gn+1t(a,t)+gn+1a(a,t)=βg(a)gn+1(a,t);for a0, t0 (2.7)
    fn+1t(a,t)+fn+1a(a,t)=βf(a)fn+1(a,t),;for a0, t0 (2.8)

    subject to the boundary conditions:

    gn+1(a,0)=g0(a)fora0, (2.9)
    gn+1(0,t)=20βf(a)fn(a,t)dafort0, (2.10)
    fn+1(a,0)=f0(a)fora0, (2.11)
    fn+1(0,t)=0βg(a)gn(a,t)dafort0, (2.12)

    Notice that gn is approximating g, which is the distribution of the cells in the first stage of the cell cycle, while fn is approximating f, which is the distribution of the cells in the second stage of the cell cycle. Note that the solution value on the boundary where a0 is determined by the previous iterate. Solutions of (2.7)–(2.12) can be found via the method of characteristic equations. Using this method, we arrive at the following solution formulas. Here Bg(s1,s2):=s2s1βg(α)dα and Bf(s1,s2):=s2s1βf(α)dα.

    Case 1: For 0a0<t0

    gn+1(a0,t0)=eBg(0,a0)(20βf(α)fn(α,t0a0)dα), (2.13)
    fn+1(a0,t0)=eBf(0,a0)(0βg(α)gn(α,t0a0)dα). (2.14)

    Case 2: For 0t0<a0

    gn+1(a0,t0)=g0(a0t0)eBg(a0t0,a0), (2.15)
    fn+1(a0,t0)=f0(a0t0)eBf(a0t0,a0). (2.16)

    Note that in the region t<a, gn and fn are independent of n and satisfy (2.7)–(2.8) together with the boundary conditions (2.9) and (2.11) provided βf and βg are continuous and g0 and f0 are differentiable. Under additional assumptions, it can be shown that gn and fn satisfy (2.7)–(2.8) together with the boundary conditions (2.10) and (2.12) also for a<t. Indeed we have the following theorem, the assumptions of which may be strengthened to achieve global solutions (please see the appendix).

    Theorem 2.1. Suppose

    (i) f0 and g0 are nonnegative and continuously differentiable for a>0,

    (ii) f0L1[0,), g0L1[0,), f0L1[0,), and g0L1[0,) are finite,

    (iii) f0 and g0 are finite,

    (iv) βf(α) and βg(α) are nonnegative, bounded and continuous, and

    (v) there exists A>0, such that for every α>A, f0(α) and g0(α) are negative and increasing,

    then for T sufficiently small there exists solutions of (2.1)–(2.6) on Ω=[0,)×[0,T), continuously differentiable, except possibly on the line a=t.

    Since we have already found a solution formula for a>t, our proof is focused on the set Ω1={(a,t)|0at<T}, where the solution formula is given by

    gn(a0,t0)=2eBg(0,a0)0βf(α)fn1(α,t0a0)dα.

    Establishing the continuity and differentiability of the integral,

    0βf(α)fn1(α,t0a0)dα, (2.17)

    is our primary task. Standard textbook theorems on this topic do not directly apply due to the requirement that there exist an L1 function, M, such that, for every t |fn1(α,t)||M(α)|. For this reason, we have adopted condition (v) of Theorem 2.1, and adapted previous proofs [26,27] to work under this alternate condition. Details are presented in the appendix.

    In this section we derive the stable age and stage distributions of the model. The stage distribution of the population at time t gives the fraction of the total cells in each stage. For example, the fraction of cells in the first stage is

    0g(a,t)da0f(a,t)+g(a,t)da.

    Similarly, the fraction of cells in the second stage at time t is

    0g(a,t)da0f(a,t)+g(a,t)da.

    Meanwhile, the age distribution of cells in the first stage at time t is

    g(a,t)0g(a,t)da,

    and that of cells in the second stage is

    f(a,t)0f(a,t)da.

    If there exists a solution of (2.1)–(2.6) such that the age distributions are constant through time, those distributions are said to be stable age distributions for the model. Hence, the stable age distribution of cells in the first stage is characterized as the solution of the partial differential equation:

    0=tg(a,t)0g(a,t)da. (3.1)

    Similarly, for cells in the second stage. It can be seen that solutions of (3.1) are necessarily separable. Indeed,

    tg(a,t)0g(a,t)da=0g(a,t)dagtg(a,t)0gtda(0g(a,t)da)2 (3.2)
    =gt(a,t)0g(a,t)da+g(a,t)0ga(a,t)+βg(a)g(a,t)da(0g(a,t)da)2, (3.3)

    where we have assumed the solution is regular enough to allow differentiation through the integral. In regard to this assumption, note we have shown that such regular solutions exist. Indeed by Theorem 2.1, there exist solutions which are continuously differentiable for at, thus we may differentiate through the integrals t0g(a,t)da and t0f(a,t)da. Moreover, in the proof of Theorem 2.1 (see Appendix 1), we show that one can differentiate through the integrals tg(a,t)da and tf(a,t)da. However, the solution may not be continuous along the line a=t. In case, 0gtda should be replaced 0gtda+limatg(a,t)limat+g(a,t). Similarly, 0ftda should be replaced by 0ftda+limatf(a,t)limat+f(a,t). Because the later expressions, like the former expressions, are independent of a, this substitutions does not alter the proof that follows.

    Now, setting the right-hand-side of (3.3) to zero and supposing, in addition, that limag(a,t)=0 (which follows from conditions ii and v of Theorem 2.1) we find

    gt(a,t)g(a,t)=g(0,t)0βg(a)g(a,t)da0g(a,t)da (3.4)

    Integrating through the equality (with respect to t) and noting that the right-hand-side is independent of a, we find

    lng(a,t)g(a,0)+c(a)=t0g(0,t)0βg(a)g(a,t)da0g(a,t)dadt:=H(t), (3.5)

    where c(a) is the 'constant' of integration with respect to t (note it may depend on a). This yields

    g(a,t)=g(a,0)ec(a)eH(t). (3.6)

    Thus, a solution g for which the age distribution is stable is necessarily separable. Similarly, if f exhibits a stable age distribution, then f must be separable. Hence, in characterizing the stable age distribution we express g and f as g(a,t)=g1(t)g2(a) and f(a,t)=f1(t)f2(a). Under this assumption, the partial differential equation for g simplifies as

    g1(t)g2(a)+g1(t)g2(a)=βg(a)g1(t)g2(a), (3.7)

    which separates as

    g1(t)g1(t)=g2(a)g2(a)βg(a). (3.8)

    The previous equation yields

    g1(t)=g1(0)ecgt, (3.9)
    g2(a)=g2(0)ea0cg+βg(s)ds, (3.10)

    where cg is a constant representing the stable growth rate of the population in stage g. Similarly, we find

    f1(t)=f1(0)ecft, (3.11)
    f2(a)=f2(0)ea0cf+βf(s)ds, (3.12)

    where cf is a constant representing the stable growth rate of the population in stage f.

    Note that the value of the product g1(0)g2(0) uniquely determines the solution (g(a,t)=g1(t)g2(a)) given by (3.9) and (3.10). This being the case, we set g1(0)=1. Similarly, we take f1(0)=1. Then, by (2.4)

    g(0,t)=ecgtg2(0)=20βf(a)ecftf2(a)da. (3.13)

    From which it follows that cg=cf:=c. Dividing (3.13) by ect we have

    g2(0)=20βf(a)f2(a)da=20βf(a)f2(0)e(a0c+βf(s)ds)da. (3.14)

    Similarly,

    f2(0)=0βg(a)g2(0)e(a0c+βg(s)ds)da. (3.15)

    Substituting the later into the former yields the consistency condition

    12=(0ea0c+βg(s)dsβg(a)da)(0ea0c+βf(s)dsβf(a)da). (3.16)

    This equation can be written more simply in terms of the probability densities of maturation times, Ig(a) and If(a). Indeed, using

    βg(a)=Ig(a)aIg(s)ds=dda(lnaIg(s)ds) (3.17)

    (see Methods), it follows that

    Ig(a)=βg(a)aIg(s)ds=βg(a)ea0βg(s)ds. (3.18)

    Similarly, we find

    If(a)=βf(a)ea0βf(s)ds. (3.19)

    Substituting the previous two expressions into 3.16 we find

    12=(0eacIg(a)da)(0eacIf(a)da)=L{Ig}(c)L{If}(c), (3.20)

    which we wish to solve for c. Here we consider two cases, in which Ig and If are exponential or inverse Gaussian probability density functions. In the case that Ig and If are exponential probability density functions, it is easy to solve (3.20) for c. Doing so, we obtain

    c=(βg+βf)+(βg+βf)2+4βgβf2, (3.21)

    where βg and βf are the constant per capita maturation rates for stages g and f. In the case that Ig and If are inverse Gaussian probability density functions, we prefer to obtain the solution of (3.20) numerically.

    In summary, to compute the stable age distribution one first solves (3.20) to obtain the asymptotic growth rate, c. Substituting c for cg and 1 for g2(0) in (3.10) we find g2(a). Next, f2(0) is computed using (3.15), and this value is used to compute f2(a) according to (3.12) with c in place of cf. Note that in this procedure, one can vary the population size by altering the value of g2(0). Figures 8 and 9 below display the models' stable age distributions and verify their stability.

    In the following section, we will examine how stable age distributions, stable stage distributions, and growth rates vary between these two models. In so doing it is helpful to note that the stable solution has a single degree of freedom g2(0), which determines the initial population size.

    In simulating the model, we need to compute the per capita maturation rates, βf and βg, which can be calculated in terms of the maturation time probability densities If and Ig. For this, let R(a,y0) denote the probability that a cell matures (transitions) to the next stage after age a, given that the cell's internal state had value y0 at a=0. We may sometimes fix y0 and just write R(a). Now let β(a)δa+o(δa) (where limδa0o(δa)δa=0) be the probability that a cell matures over the interval [a,a+δa], given that it has not maturation at age a. That is, β(a) is the maturation rate. Then on the one hand

    R(a+δa)=R(a)(1β(a)δao(δa)).

    That is, the probability that a cell matures after age, a+δa, is the probability that the cell does not mature over [a,a+δa], given it did not mature up until age a, times the probability that the cell did not mature up until age a. On the other hand,

    R(a+δa)=R(a)+R(a)δa+o(δa).

    Equating these two expressions for R(a+δa) and canceling like terms, we have

    β(a)R(a)δa=R(a)δa+o(δa).

    Dividing by δa and taking the limit as δa goes to zero, we find

    β(a)=R(a)R(a).

    Thus, we can determine the maturation probability in terms of R(a). Note that, R(a)=aI(s)ds. Therefore

    β(a)=I(a)aI(s)ds. (4.1)

    In general, there may not be a closed form for the maturation rate, β(a), so that it must be approximated numerically. In this case, as a grows, the numerator and denominator in the expression for β(a) approach zero, and the computation is challenging due to limits on floating-point precision. As a result, it is useful to have a better qualitative characterization of β(a) for the purpose of validating numerical simulations. In the following paragraphs, we characterize several important features of the β(a) which corresponds to the inverse Gaussian probability density.

    First note we may use L'Hôpital's rule to compute the asymptotic value of β(a). For the inverse Gaussian distribution

    I(a):=12σ2πa3e(μa1)22aσ2,

    so that

    limaβ(a)=limaI(a)aI(s)ds=limaI(a)I(a)=lima321a12σ21a2+μ22σ2=μ22σ2. (4.2)

    Our characterization of β also involves the ratio

    q(a)=I(a)I(a)=321a12σ21a2+μ22σ2.

    Theorem 4.1. The per capita maturation rate, β(a), has exactly one critical point, at which it takes a maximum value. Moreover, if a is the age at which β(a) takes its maximum value, β(a)<q(a) for a>a.

    Proof. With β(a)=I(a)aI(s)ds and I(a)=I(a)q(a), where q(a)=321a12σ21a2+μ22σ2, as above we find that

    β(a)=aI(s)dsI(a)q(a)+I2(a)(aI(s)ds)2<0 (4.3)
    β(a)<q(a). (4.4)

    Similarly,

    β(a)>0β(a)>q(a). (4.5)

    and

    β(a)=0β(a)=q(a). (4.6)

    Also note that as a0, β(a)0, and q(a). Therefore, β(a)>q(a) for a small, i.e., β(a)>0 for a small, i.e., β(a) is initially increasing.

    To show that β has a single critical point, we must also consider the behavior of q(a). Note that

    q(a)=321a2+1σ21a3=032a+1σ2=0a=231σ2.

    Set ˆa=231σ2. Considering the limits of q(a) as a0 and a we see that q(a)>0 for a<ˆa and q(a)<0 for a>ˆa. Therefore q(a) takes its maximum value at ˆa=23σ2.

    Suppose toward a contradiction that β(a)0 for all a>0, then by continuity, β(a)>0 for a>0. Thus, β(a)>q(a) for a>0 by (4.5). Fixing a0>ˆa>0, we have β(a)>q(a) and β(a)>0>q(a) for a>a0>ˆa. Thus,

    0=limsβ(s)+q(s)=limssa0β(a)+q(a)da+β(a0)+q(a0)>β(a0)+q(a0)>0, (4.7)

    which is a contradiction. Therefore, we may define a as follows,

    a:=inf{a>0β(a)=0}.

    Note by continuity β(a)=0 and β(a)=q(a). Note also that a>0, since β(a)>0 for a small. Therefore, by continuity, β(a)>0 for a<a, i.e., β is increasing for a<a.

    Case 1:

    Suppose that 0=β(a)>q(a). Then by continuity there exists an interval over which the inequality holds. Moreover, since β(a)=q(a), there exists δ>0 such that β(a)>q(a) for a<a<a+δ. Suppose toward a contradiction that there exists a>a such that β(a)q(a). Let s=inf{a>aβ(a)q(a)}, and note sa. Thus for a<a<s, β(a)>q(a), and for a<a<s, β(a)>0. However, since q has a single critical point at which is takes a maximum, we know q(a)<0 for a>a. Therefore β(a)>q(a) for a<a<s, and

    β(s)β(a)=saβ(a)da>saq(a)da=q(s)β(a). (4.8)

    That is, β(s)>q(s). However, by continuity and the definition of s, it must be that β(s)q(s). Thus we have reached a contradiction.

    It follows that {aa>a,β(a)q(a)} is empty. That is, β(a)>q(a) for a>a, i.e., β(a)>0 for a>a. Since we have already noted q(a)<0 for a>a, we find that

    0=limsβ(s)+q(s)=limssaβ(a)+q(a)da>0. (4.9)

    Hence, this case does not occur.

    Case 2:

    Suppose that 0=β(a)<q(a). Then, as in Case 1, there exists δ such that β(a)<q(a) for a<a<a+δ. Suppose toward a contradiction that there exists a>a, so that β(a)q(a). Let s=inf{a>aβ(a)q(a)}. Note we have s>a, and by continuity β(s)=q(s). However for a<a<s, β(a)<q(a), and hence for a<a<s, β(a)<0. Since we have already shown Case I cannot happen, we also know that 0=β(s)q(s). Since q has a single maximum, we see that in fact β(a)<0<q(a) for a<a<s. Thus β(s)<q(s), and we have reached a contradiction. Thus, there exists no a>a such that β(a)q(a). Therefore, there exists a unique age a at which β(a)=0, and β(a)<q(a) for a>a.

    Case 3:

    Assume that 0=β(a)=q(a). Then since q(a) has a single critical point at which it takes a maximum value, q(a)<0 for a>a. Therefore, it cannot happen that β(a)=q(a), for a>a, since we previously showed Case I cannot happen. Thus, in this case too, we see that there is a unique time a at which β(a)=0. Moreover, if there exists s>a so that β(s)>q(s), then by continuity and because β(a)q(a) for a>a, it must be that β(a)>q(a) for every a>a. Therefore, β(a)>0>q(a) for a>a. Contradicting that limaβ(a)=q(a). Thus, it must be that β(a)<q(a) for a>a, as desired

    Hence we have shown that, in any case, there is a unique age a at which β(a)=0. Moreover, β(a)<q(a) for a>a, so that β(a) is decreasing for a>a. Since we have already noted that β(a) is increasing for a>a we see that β(a) takes its maximum value at a as desired.

    Next we derive an estimate of the maximum value of β(a) and a lower bound on the age at which β(a) assumes its maximum value. For this note that since q(a)0 and β(a)<0 for a<a

    max{a>0}β(a)max{a>0}(q(a))=q(ˆa) (4.10)
    =98σ2+μ22σ2 (4.11)

    Lemma 4.2. For a as above, μ22σ2<β(a)98σ2+μ22σ2. Moreover, 131σ2<a.

    Proof. To obtain the estimate on β(a)=max{a>0}β(a) and the lower bound on a, we first consider the unique time, a, such that q(a)=μ22σ2. That is, we consider the unique finite time at which q(a) achieves the asymptotic value of β(a).

    We have

    q(a)=μ22σ2 (4.12)
    0=321a12σ21a2 (4.13)
    0=32a12σ2 (4.14)
    a=131σ2 (4.15)

    It follows that 131σ2=a<a. Indeed, since a<ˆa=23σ2, q is increasing for a<a. Were a<a, we would have β(a)=q(a)<q(a)=μ22σ2. However, this leads to a contradiction because β(a) is strictly decreasing for a>a and approaches μ22σ2 as a. Therefore, a<a as desired. Hence, β(a)>β(a)>q(a)=μ22σ2, where the final inequality follows from the fact that β(a) is increasing (i.e., β(a)>q(a) for a<a.)

    The practical computation of the maturation rate curve presents difficulties using fixed precision floating point techniques. This is due to the behavior of both the numerator and denominator of β(a) as a increases causing the result to become unstable and ultimately undefined.

    β(a)=I(a)aI(s)ds. (4.16)
    limaI(a)=0 (4.17)
    limaaI(s)ds=0 (4.18)

    While variable precision arithmetic, as provided by MATLAB, allows for the stable computation of the maturation rate, as the numerator and denominator become small, the required precision increases, which in turn increases the computational expense. Runtimes with an acceptably small mesh sizes were found to be excessive for our choice of parameters. This was resolved by using interpolation to reduce the number of points that need to be evaluated with variable precision arithmetic. This is justified as the maturation rate curve is smooth (both the numerator and denominator are smooth and nonzero). In the absence of a predetermined method for choosing knots, we use an adaptive method: At each iteration, we double the number of knots by evaluating points midway between the existing knots, remembering previously computed knots using a memoization technique. Having the knots, we then interpolate for each point on our desired mesh. The vector of approximate maturation rate values is compared with the previous vector of approximate maturation rate values using the 2-norm of the difference. When this error is reduced below a specified threshold, here chosen to be 0.1, the solution is accepted and returned. In this way, an accurate approximation of the maturation rate curve is obtained with acceptable performance (See Figures 1 and 2).

    Figure 1.  Convergence of maturation rates (β) as the number of knots increases.
    Figure 2.  Adaptive error reduction with increasing numbers of knots. Blue denotes rejected solutions and red denotes accepted solution.

    In order to solve the system of PDEs numerically we discretize age and time in order to numerically integrate along the model's characteristic curves. The numerical scheme is summarized below.

    Given initial data g0(a), f0(a), with g0(a)=f0(a)=0 for a>amax, we approximate the solution of our system for t<T as follows:

    Let a=(0,h,2h,Ih), and t=(0,h,2h,Jh), where Jh=T, and Ih=T+amax. Now define ˆg=(ˆgij)R(I+1)×(J+1) and ˆf=(ˆfij)R(I+1)×(J+1) as matrices such that

    ˆfi0=f0(ai),i=1,I, (4.19)
    ˆgi0=g0(ai)i=1,I, (4.20)
    ˆf0j=integral(0,aI,βgˆg:,j)j=1,J, (4.21)
    ˆg0j=2 integral(0,aI,βfˆf:,j)j=1,J, (4.22)
    ˆfi+1,j+1=ˆfijexp{integral(ai,ai+1,βf)} (4.23)
    ˆgi+1,j+1=ˆgijexp{integral(ai,ai+1,βg)} (4.24)

    where integral(ai,ai+1,βf) is the approximation of ai+1aiβf(α) dα using MATLAB's implementation of the trapezoid rule with a uniform grid over [ai,ai+1] with four points, and integral(0,aI,βfˆf:,j) is the approximation of aI0βf(α)ˆf(α,tj) dα using the MATLAB's implementation of the trapezoid rule with the grid values in a and the corresponding entries of ˆf:,j. The grid size h was initially set to .02 and was reduced by half until the relative point-wise error was less than 102. This algorithm is an extension of that presented in [22] to a PDE system. The consistency and stability of the algorithm can be shown similarly to [22], so we omit the details here. All simulations were performed in MATLAB. These codes can be found in the following repository: https://github.com/rnleander/Asma.

    Here we simulate the model using both inverse Gaussian and exponential maturation time distributions. Parameters for the inverse Gaussian distribution were chosen to fit data on the division times of MCF10A cells, as described in [20]. The parameters for the exponential distributions were chosen to match the mean of the inverse Gaussian distribution as parameterized by the MCF10A cell data. Thus the average time spent in early G1, late G1-M, and the full cell cycle are the same for both models.

    Figure 3.  Stage maturation times, Ig and If, and total intermitotic time, (IgIf)(a), distributions for both exponential and inverse Gaussian models.
    Table 1.  Means and variances of model stages.
    Model Meang Varianceg Meanf Variancef
    Exponential 4.00 16.0 15.6 244
    Inverse Gaussian 4.00 64.9 15.6 367

     | Show Table
    DownLoad: CSV

    To explore the models' transient dynamics, we initialize them with their stable stage structures and two unstable age densities (see Section 3); the uniform initial age density (Figures 4 and 6) and the Gaussian initial age density (Figure 5). For the inverse Gaussian model, the stage and age structures oscillate over multiple days when initialized with either unstable age density. However, the amplitude of the oscillations decreases through time, so that the population stage structure approaches the stable stage structure (Figures 4 and 5). In contrast, for the exponential model, the stage structure is insensitive to the initial age density: No matter the initial age density, the stage structure does not evolve through time provided it is initialized as the stable stage structure. The age distributions behave similarly. Figures 47 demonstrate that the inverse Gaussian age distribution stabilizes much more slowly than that of the exponential distribution. We also initialized the models with stable age densities and perturbed stage structures. For these simulations, we perturbed the ratio of cells in early G1 to late G1-M by a factor of two. Again, the exponential model stabilizes much more quickly that the inverse Gaussian model (Figures 11 and 12).

    Figure 4.  Inverse Gaussian model with uniform initial age density and stable initial age structure. Note that the stage structure oscillates as the age structure evolves.
    Figure 5.  Inverse Gaussian model with Gaussian initial age density and stable initial age structure. Note that the stage structure oscillates as the age structure evolves.
    Figure 6.  Exponential model with uniform initial age density and stable initial age structure. Note that the stage structure does not evolve, and the age density stabilizes quickly.
    Figure 7.  Exponential model with Gaussian initial age density and stable initial age structure. Note that the stage structure does not evolve, and the age density stabilizes quickly.
    Figure 8.  Inverse Gaussian model with stable initial age density and stage structure. Note that when initialized with the stable age density and stable stage structure, neither the stage structure nor the age density evolve through time.
    Figure 9.  Exponential model with stable initial age density and stage structure. Note that the stage structure and age density do not evolve through time.

    The two models also differ in their predicted stable stage distributions (Figures 8 and 9). Despite the fact that, as parameterized, the average time spent in early G1 and late G1-M is identical for both models, the stable exponential stage structure has about 30% of cells in early G1, while the stable inverse Gaussian stage structure has about 22% of cells in early G1.

    Finally Figure 10 and Table 2 present the predicted and observed asymptotic growth rates of the population for the two maturation time models. It is interesting to note that, despite exhibiting the same average intermitotic time, the models differ substantially in terms of their asymptotic growth rates. As a result, after 100 hours, a population with exponential maturation times is expected to be almost twice as large as that with inverse Gaussian maturation times.

    Figure 10.  Growth rate curves and exponential growth fits. Note that while oscillation in the inverse Gaussian model does affect the total population, exponential growth remains a good fit.
    Table 2.  Predicted and fitted exponential growth rates in hr1. Note that there is some variability in the fitted growth rates for the Inverse Gaussian model with uniform and Gaussian initial conditions due to oscillations that are not present with the stable initial conditions.
    Model Predicted Uniform Gaussian Stable
    Exponential 0.04462 0.04455 0.04455 0.04455
    Inverse Gaussian 0.03689 0.03673 0.03694 0.03685

     | Show Table
    DownLoad: CSV

    Our results demonstrate that maturation time distributions can significantly impact the transient dynamics, growth, and structure of a cellular population. Indeed, when employing the ageless exponential maturation time distribution, the population stage structure is insensitive to the initial age density, in that perturbations in the age density do not perturb the stage structure. Moreover, the stage structure and age density stabilize quickly on perturbation. In contrast, the inverse Gaussian model's stage structure exhibits sustained oscillations over several days when the age desnity is perturbed away from the stable age density. This difference represents a major qualitative distinction between the two models. Importantly, the models also differ in their mean-predicted growth rates. That is, when parameterized so that the average time spent in each stage, and thus the average time spent in the full cell cycle, is the same for both models, the asymptotic growth rate of the exponential model is considerably greater. This difference is likely due to the fact that few cells exhibit short maturation times under the inverse Gaussian model. In contrast, short maturation times are common in the exponential model. Thus, the difference in growth rates is likely the result of variability in the time to mature. As a result of these differences, we recommend that modelers take care when choosing a maturation time distribution. Indeed, we hope that our methods and programs can help researchers develop and simulate more accurate, predictive models of cellular population growth. We note that the cell cycle model presented here is well-suited to study the impact of CDK inhibitors, which impact restriction point passage (in our model, maturation into late G1) [23]. The model could also be adapted to study the impact of drugs which target S phase (e.g., gemcitabine) [3]. A second benefit of this work is that it provides testable predictions that can be used to determine if the exponential distribution provides a reasonable approximation to a population's maturation time distribution. In particular, the rapidity with which the exponential model achieves a stable stage distribution is predicted to be a robust feature of the model that would be readily observable via experimentation. In summary, we find evidence that a population's maturation time distribution can inform model development and refinement while deepening our understanding of the fundamental process of cellular proliferation.

    The authors would like to thank Dr. Glenn Webb for sharing his Mathematica Code.

    All authors declare no conflicts of interest in this paper.

    In proving Theorem 2.1 we consider C(Ω1), the Banach space of continuous, bounded functions on Ω1 with the norm

    ||h||:=supxΩ1h(x). (6.1)

    We begin by establishing the following lemma.

    Lemma 6.1. For f0, g0, βf(α), and βg(α) as in Theorem 2.1

    (i) gn and fnC(Ω1) and

    (ii) g:=limngn and f:=limnfn belong to C(Ω1).

    Proof. The proof is by induction. First note that for (a0,t0)Ω1,

    gn(a0,t0)=2eBg(0,a0)0βf(α)fn1(α,t0a0)dα.

    Since eBg(0,a0) is continuous, it suffices to show that 0βf(α)fn1(α,t0a0)dα is continuous in Ω1.

    Suppose fn1(a,t) is continuous in Ω1. Choose (a0,t0)Ω1, let ϵ>0, and let (ak,tk) be a sequence of points converging to (a0,t0) in Ω1.

    The integral of interest may be written as:

    |0βf(α)fn1(α,tkak)dα0βf(α)fn1(α,t0a0)dα|0βf(α)|fn1(α,tkak)fn1(α,t0a0)| dα=A+T0βf(α)|fn1(α,tkak)fn1(α,t0a0)| dα+A+Tβf(α)|fn1(α,tkak)fn1(α,t0a0)| dα, (6.2)

    where A is chosen so that

    A|f0(α)| dαϵ2βf.

    Note by our assumptions, tkak<T and t0a0<T. Hence, there exists τ<T so that tkakτ for k=0,1,. Also, since fn1(a,t) is continuous on the closed and bounded sets D1={(a,t)|0atτ} and D2={(a,t)|0taA+T,0tτ}, there exists a constant C so that fn1(a,t)<C on D1D2, and, for every α[0,t0a0)(t0a0,A+T], as (ak,tk)(a0,t0), fn1(α,tkak)fn1(α,t0a0). Hence, by Lebesgue's Dominated Convergence Theorem [24],

    A+T0βf(α)|fn1(α,tkak)fn1(α,t0a0)| dα0.

    Now note that the final term to the right of the equality in (6.2) is bounded by our choice of A. Indeed,

    βfA+T|f0(α(tkak))eBf(α(tkak),α)f0(α(t0a0))eBf(α(t0a0),α)| dαβfA+T|f0(α(tkak))| dα+βfA+T|f0(α(t0a0))| dαϵ (6.3)

    Since ϵ was arbitrary, we see that gn is continuous in Ω1. Similarly, it can be shown that fn is continuous in Ω1. This ends the proof that fn and gn are continuous in Ω1.

    Now we show that fn and gn have limits in C(Ω1). For a<t and n=0 we have

    |g1g0|(a,t)=NeBg(0,a)0f0(α)dαg0(a) (6.4)
    Nf0L1+g0, (6.5)
    |f1f0|(a,t)Ng0L1+f0. (6.6)

    Where N:=max{2βg,βf,2βf2,2βg2,2βfβg}. Hence, we can define M:=max{f1f0,g1g0}<.

    In general,

    gn+1(a,t)gn(a,t)=eBg(0,a)(02βf(α)(fn(α,ta)dαfn1(α,ta))dα). (6.7)

    Thus, for a<t<T

    |gn+1(a,t)gn(a,t)|=|eBg(0,a)02βf(α)(fn(α,ta)fn1(α,ta))dα|eBg(0,a)ta02βf(α)|fnfn1|dα=NeBg(0,a)(ta)fnfn1NTfnfn1 (6.8)

    That is,

    gn+1gnNTfnfn1. (6.9)

    Similarly,

    fn+1fnNTgngn1. (6.10)

    (Note that by the triangle inequality and since fn and gn are finite, fn+1 and gn+1 are finite as well.) We now have that

    gn+1gnNTfnfn1 (6.11)
    (NT)n1M, (6.12)

    and

    fn+1fnNTgngn1 (6.13)
    (NT)n1M. (6.14)

    Since

    gn=g0+ni=1(gigi1), (6.15)

    we see that

    g:=limngn=g0+i=1(gigi1) (6.16)

    exists, and the convergence is uniform on Ω1 by the Weierstrass M-test, provided

    T<1N. (6.17)

    Therefore, gC(Ω1). Similarly,

    f:=limnfn=f0+i=1(fifi1) (6.18)

    exists in C(Ω1). This concludes the proof of convergence, so we have established Lemma 2.2.

    Assuming that we may differentiate through the integral in (2.13) - (2.14), and accounting for the possible discontinuity at a=t we find:

    Case 1 : For 0a<t

    gn+1a(a,t)=2eBg(0,a)0βf(α)(βg(a)(fn(α,ta)+fnt(α,ta))dα2βf(ta)eBg(0,a)(limα(ta)fn(α,ta)limα(ta)+fn(α,ta)) (6.19)
    gn+1t(a,t)=2eBg(0,a)0βf(α)fnt(α,ta)dα+2βf(ta)eBg(0,a)(limα(ta)fn(α,ta)limα(ta)+fn(α,ta)) (6.20)
    fn+1a(a,t)=eBf(0,a)0βg(α)(βf(a)gn(α,ta)+gnt(α,ta))dαβg(ta)eBf(0,a)(limα(ta)gn(α,ta)limα(ta)+gn(α,ta)) (6.21)
    fn+1t(a,t)=eBf(0,a)0βg(α)gnt(α,ta)dα+βg(ta)eBf(0,a)(limα(ta)gn(α,ta)limα(ta)+gn(α,ta)) (6.22)

    Case 2 : For 0t<a

    gn+1a(a,t)=eBg(at,a)(g0(at)+βg(at)g0(at)βg(a)g0(at)) (6.23)
    gn+1t(a,t)=eBg(at,a)(g0(at)βg(at)g0(at)) (6.24)
    fn+1a(a,t)=eBf(at,a)(f0(at)+βf(at)f0(at)βf(a)f0(at)) (6.25)
    fn+1t(a,t)=eBf(at,a)(f0(at)βf(at)f0(at)). (6.26)

    From the above cases we see that gn and fn given by (2.13)-(2.16) will satisfy (2.7) - (2.8) together with the boundary conditions (2.9) - (2.12), provided we may differentiate through the integral.

    Shortly we will show that the first partial derivatives of gn and fn are given by (6.19)-(6.22), but first we will establish the continuity of the expressions to the right of each equality in (6.19)-(6.22). Moreover, we will show these expressions converge uniformly in Ω1. For convenience, we refer to the integral expressions above as the partial derivatives of fn and gn; however, in the following lemma and proof, we do not assume this to be the case. That is, in the following lemma gna,gnt,fna, and fnt stand for the expressions on the right-hand-side of (6.23)-(6.26), respectively.

    Lemma 6.2. Let f0, g0, βf(α) and βg(α) as in Theorem 2.1, and define gnt, fnt, gna, and fna by (6.20), (6.22), (6.19) and (6.21), respectively.

    (i) gnt, fnt, gna, and fna belong to C(Ω1)

    (ii) fnt, fna, gnt and gna converge uniformly on Ω1.

    Proof. The continuity of gnt, fnt, gna, and fna on Ω1 follows by induction as in the proof of Lemma 2.2.

    Now we show the sequences fnt, fna, gnt and gna converge uniformly on Ω1. (Note that these sequences are constant for a>t, and hence convergence is uniform in this region as well.)

    We see that

    |gn+1tgnt|NeBg(0,a)ta0|fntfn1t|dα (6.27)
    +NeBg(0,a)|limα(ta)fn(α,ta)limα(ta)+fn1(α,ta)|N(ta)fntfn1t+Nfnfn1 (6.28)
    NTfntfn1t+N(NT)n1M, (6.29)

    and

    |fn+1tfnt|NeBf(0,a)ta0|gntgn1tdα| (6.30)
    +NeBf(0,a)|limα(ta)gn(α,ta)limα(ta)+gn1(α,ta)|N(ta)gntgn1t+Ngngn1 (6.31)
    NTgntgn1t+N(NT)n1M. (6.32)

    Combining these two together :

    gn+1tgntNTfntfn1t+N(NT)n1M (6.33)
    fn+1tfntNTgntgn1t+N(NT)n1M (6.34)

    Let ˆM:=max{f2tf1t,g2tg1t}. Then,

    gn+1tgntNTfntfn1t+N(NT)n1M (6.35)
    ˆM(NT)n1+(n1)NM(NT)n1. (6.36)

    Similarly,

    fn+1tfntˆM(NT)n1+(n1)NM(NT)n1. (6.37)

    Thus, provided ˆM is finite, the sequences of partial derivatives converge uniformly for tT<1N.

    To show that ˆM is finite, we first consider the base case. For this, it is useful to recall

    f0(a,t)=f0(a). (6.38)
    g0(a,t)=g0(a), (6.39)

    Also, for n>0 and a<t,

    gn+1t(a,t)=2eBg(0,a)0βf(α)fnt(α,ta)dα+2βf(ta)eBg(0,a)(limα(ta)fn(α,ta)limα(ta)+fn(α,ta)) (6.40)
    fn+1t(a,t)=eBf(0,a)0βg(α)gnt(α,ta)dα+βg(ta)eBf(0,a)(limα(ta)gn(α,ta)limα(ta)+gn(α,ta)), (6.41)

    while for n>0 and t<a

    gn+1t(a,t)=eBg(at,a)(g0(at)βg(at)g0(at)) (6.42)
    fn+1t(a,t)=eBf(at,a)(f0(at)βf(at)f0(at)). (6.43)

    From (6.38) and (6.39)

    g0t(a,t)=0for(a,t)(0,)×(0,T) (6.44)
    f0t(a,t)=0for(a,t)(0,)×(0,T). (6.45)

    Also, f0(a,t)f0(a) is continuous. Hence by (6.40) and (6.41)

    g1t(a,t)=0fora<t (6.46)
    f1t(a,t)=0fora<t. (6.47)

    On the other hand, for t<a, g1t and f1t are given by (6.42) and (6.43), respectively. Having computed g1t and f1t we are ready to compute g2t and f2t. For a<t:

    g2t(a,t)=2eBg(0,a)(ta0βf(α)f1t(α,ta)dα+taβf(α)f1t(α,ta)dα)+2βf(ta)eBg(0,a)(limα(ta)f1(α,ta)limα(ta)+f1(α,ta))=2eBg(0,a)taβf(α)βf(α(ta))eBf(α(ta),α)f0(α(ta))dα2eBg(0,a)taβf(α)eBf(α(ta),α)f0(α(ta))dα+2βf(ta)eBg(0,a)eBf(0,ta)(0βg(α)g0(α)dαf0(0)) (6.48)

    For the first two terms to the right of the final equality above, let u=α(ta) and du=dα, so that we obtain:

    g2t(a,t)=2eBg(0,a)0eBf(u,u+(ta))βf(u+(ta))βf(u)f0(u)du2eBg(0,a)0eBf(u,u+(ta))βf(u+(ta))f0(u)du. (6.49)

    Thus for a<t

    |g2t(a,t)|2βf2f0L1+2βf f0L1+2βfβgg0L1+2βff0N( f0L1+ g0L1+ f0L1+ f0). (6.50)

    Therefore,

    g2tg1t<. (6.51)

    Similarly,

    f2tf1t<. (6.52)

    This shows that ˆM is finite, and the partial derivatives with respect to t converge uniformly to their limits in Ω1 for T<1N. Furthermore, from (6.19) and (6.21) we see that

    gn+1a(a0,t0)=gn+1t(a0,t0)βg(a0)gn+1(a0,t0), (6.53)
    fn+1a(a0,t0)=fn+1t(a0,t0)βf(a0)fn+1(a0,t0). (6.54)

    Therefore, the uniform convergence of gna and fna follows from that of gnt, fnt, fn and gn.

    Now we will show that for every nN, fn and gn are continuously differentiable with respect to t. Moreover, we can compute fnt and gnt by differentiating through the integral in (2.13) and (2.14). The proof is by induction.

    Proof. Suppose that fn is continuously differentiable with respect to t in Ω1. Let (a,t)Ω1 and choose δ>0 so that (a,t±δ)Ω1 (or, in case a=t, (a,t+δ)Ω1). Also, suppose δ>Δt>0. Since Ω1 is convex, we see that (a,t±Δt)Ω1 for any such Δt. Given ϵ>0, suppose that A>A is chosen such that

    βfA|f0(α)|dα+βf2A|f0(α)|dα<ϵ2.

    This is possible since f0 and f0 are L1. Now we establish convergence of the difference quotient:

    0βf(α)(fn(α,t+Δta)fn(α,ta))Δtdα=A+T0βf(α)(fn(α,t+Δta)fn(α,ta))Δtdα+A+Tβf(α)(fn(α,t+Δta)fn(α,ta))Δtdα (6.55)

    We will handle the first and second terms to the right if the inequality in (6.55) separately. The first term can be expressed as

    A+T0βf(α)(fn(α,t+Δta)fn(α,ta))Δtdα=ta0βf(α)(fn(α,t+Δta)fn(α,ta))Δtdα+ta+Δttaβf(α)(fn(α,t+Δta)fn(α,ta))Δtdα+A+Tta+Δtβf(α)(fn(α,t+Δta)fn(α,ta))Δtdα (6.56)

    Since fn(α,t) and fnt(α,t) are continuous on the sets D1={(α,t)|0αtta+δ} and D2={(α,t)|0taA+T,0tta+δ}, by the mean value theorem, the first and third terms to the right of the equality in (6.56) can be expressed as:

    ta0βf(α)fn(α,t+Δta)fn(α,ta)Δtdα=ta0βf(α)fnt(α,t(α))dα,
    A+Tta+Δtβf(α)fn(α,t+Δta)fn(α,ta)Δtdα=A+Tta+Δtβf(α)fnt(α,t(α))dα.

    Where t(a) is between ta and ta+Δt. Since fnt(α,t) is continuous on D1 and D2, we have that fnt(α,t(α))fnt(α,t) point-wise as Δt0. Moreover, since D1 and D2 are closed and bounded, there exists a constant C so that fnt(α,t)<C on D1D2. Therefore, by Lebesgue's dominated convergence theorem, as Δt0,

    ta0βf(α)fnt(α,t(α))dαta0βf(α)fnt(α,ta)dα

    and

    A+Tta+Δtβf(α)fnt(α,t(α))dαA+Tta+Δtβf(α)fnt(α,ta)dα.

    For the second term in (6.56) we have

    ta+Δttaβf(α)fn(α,t+Δta)fn(α,ta)Δtdα=1Δtta+Δttaβf(α)fn(α,t+Δta)dα1Δtta+Δttaβf(α)fn(α,ta)dα (6.57)

    Since fn(α,t) is uniformly continuous on D1, which contains the domain of integration for the first integral above, as Δt0,

    1Δtta+Δttaβf(α)fn(α,t+Δta)dαlimα(ta)βf(α)fn(α,ta).

    Since fn(α,t) is uniformly continuous on the closed and bounded region D2, which contains the domain of integration for the second integral above, as Δt0,

    1Δtta+Δttaβf(α)fn(α,ta)dαlimα(ta)+βf(α)fn(α,ta).

    Hence

    A+T0βf(α)(fn(α,t+Δta)fn(α,ta))ΔtdαA+T0βf(α)fnt(α,ta)dα+limα(ta)βf(α)fn(α,ta)limα(ta)+βf(α)fn(α,ta).

    Now we show the convergence of the second term in (6.55). Applying the mean value theorem to fn(α,t(α)) for αt,

    |A+Tβf(α)(fn(α,t+Δta)fn(α,ta))Δtdα|=|A+Tβf(α)fnt(α,t(α))dα|βfA+T|fnt(α,t(α))|dα (6.58)

    where t(a) is between ta and ta+Δt. The integral in the final term can be expanded as:

    βfA+T|fnt(α,t(α))|dα=βfA+T|eBf(αt(α),α)(f0(αt(α))+βf(αt(α))f0(αt(α)))|dαβfA+T|f0(αt(α))|dα+βf2A+T|f0(αt(α))|dα. (6.59)

    Since (αt(α))(α(t+Δta))A, by (v) of Theorem 2.1, we have that |f0(αt(α))||f0(α(t+Δta))| and |f0(αt(α))||f0(α(t+Δta))|. Thus,

    βfA+T|f0(αt(α))|dα+βf2A+T|f0(αt(α))|dαβfA+T|f0(α(t+Δta))|dα+βf2A+T|f0(α(t+Δta))|dαβfA|f0(α)|dα+βf2A|f0(α)|dαϵ2. (6.60)

    In summary,

    |A+Tβf(α)(fn(α,t+Δta)fn(α,ta))Δtdα|<ϵ2.

    Similarly

    |A+Tβf(α)fnt(α,ta)dα|βfA+T|eBf(α(ta),α)(f0(α(ta))+βf(α(ta))f0(α(ta)))|dαβfA+T|f0(α(ta))|dα+βf2A+T|f0(α(ta))|dαβfA|f0(α)|dα+βf2A|f0(α)|dαϵ2. (6.61)

    Thus,

    limΔt0|A+Tβf(α)(fn(α,t+Δta)fn(α,ta))ΔtdαA+Tβf(α)fnt(α,ta)dα|ϵ.

    Therefore the absolute value of

    0βf(α)(fn(α,t+Δta)fn(α,ta))Δtdα(0βf(α)fnt(α,ta)dα+limα(ta)βf(α)fn(α,ta)limα(ta)+βf(α)fn(α,ta))

    is less than ϵ. Since ϵ>0 was arbitrary, for any (a,t)Ω1,

    gn+1t(a,t)=2eBg(0,a)0βf(α)fnt(α,ta)dα+2βf(ta)eBg(0,a)(limα(ta)fn(α,ta)limα(ta)+fn(α,ta)). (6.62)

    From Lemma 2.3, the expression to the right of the equality above is continuous, hence we have shown that gn is continuously differentiable with respect to t in Ω1. Similarly, we find that fn and gn are continuously differentiable with respect to both t and a in Ω1, and their derivatives are given by (6.19)–(6.22). Moreover, we see that gn and fn satisfy (2.7)–(2.8). It then follows from the uniform convergence of the partial derivatives of fn and gn together with the convergence of the sequences fn and gn, that

    limnfnt=ft,

    and

    limngna=gt.

    Hence, taking the limit through (2.7)–(2.8), we see that f and g satisfy (2.1)–(2.2). Moreover, since the convergence is uniform, we have that f and g are continuously differentiable on Ω1.

    To complete the proof of Theorem 2.1, it remains to show that g and f satisfy the boundary conditions (2.4) and (2.6).

    Proof.

    Note that:

    gn+1(0,t)=20βf(a)fn(a,t)da=2t0βf(a)fn(a,t)da+2tβf(a)f0(at)eBf(at,a)da (6.63)

    Since the domain of integration for the first integral is contained in Ω1 where fn converges uniformly to f we have that:

    2t0βf(a)fn(a,t)da2t0βf(a)f(a,t)da.

    Thus,

    gn+1(0,t)2t0βf(a)f(a,t)da+2tβf(a)f0(at)eBf(at,a)da=20βf(a)f(a,t)da (6.64)

    as desired. Similarly, we find that f satisfies the boundary condition (2.6).

    Thus, the proof of Theorem 2.1 is complete.

    Theorem 6.3. Assume that in addition to conditions (i)(v) of Theorem 2.1, βg and βf are differentiable and

    (vi) For a>ˆA, βg and βf are non-postive and increasing,

    then there exist continuously differentiable solutions of (2.1)–(2.2) together with the boundary conditions (2.3)-(2.6) on Ω={(a,t)|0a,0t<T}, for all T>0.

    Proof. By Theorem 2.1 there exist solutions g and f of (2.1)–(2.2) together with the boundary conditions (2.4) and (2.6) on ˉΩ1={(a,t)|0atT}, for T=12N. We may set ˆf0(a)=f(a,T) and ˆg0(a)=g(a,T). Note that ˆf0 and ˆg0 are continuous and continuously differentiable for aT at which point they are continuous from the left and right, with a jump discontinuity. Also we have that ˆf0 and ˆg0 are L. This is because gng0 and fnf0 in ˉΩ1, and, in addition, f and g are L for 0ta by (2.15) and (2.16). Also, ˆf0 and ˆg0 are L1. This follows from the fact that ˆf0 and ˆg0 are L and given by (2.15) and (2.16) for a large, where f0 and g0 are L1.

    Also we have that ˆf0=fa(a,T) and ˆg0=ga(a,T) are L. This is because gnaga0 and fnafa0 in ¯Ω1, and, in addition, fa(a,T) and fa(a,T) are L for 0ta by (6.23) and (6.25). Also, ˆf0 and ˆg0 are L1. This follows from the fact that ˆf0 and ˆg0 are L and given by (6.23) and (6.25) for a large, where βf, βg are L and f0, f0, g0 and g0 are L1.

    Now we will verify that ˆf0 and ˆg0 satisfy condition v.

    ˆf0(a)=f0(aT)eBf(aT,a)f0(aT)[βf(aT)βf(a)]eBf(aT,a) (6.65)

    Since f0(a) satisfies v and βf is decreasing for a>ˆA we see that ˆf0(a) is non-positive for a>A:=max{A+T,ˆA+T}. Note also that the first term above,

    f0(aT)eBf(aT,a),

    is increasing for a>A. Indeed, eBf(aT,a) is decreasing for a>A, and f0(aT) is negative and increasing for a>A. Therefore for ˆa>a>A,

    f0(aT)eBf(aT,a)f0(ˆaT)eBf(aT,a)f0(ˆaT)eBf(ˆaT,ˆa)

    Now, taking the derivative of the second term,

    f0(aT)[βf(aT)βf(a)]eBf(aT,a),

    we get:

    f0(aT)(βg(aT)βg(a))eBg(aT,a)f0(aT)(βg(aT)βg(a)eBg(aT,a)+f0(aT)(βg(aT)βg(a))2eBg(aT,a),

    which is positive for a>A by conditions v and vi. Therefore ˆf0(a) is increasing for a large. Similarly ˆg0(a) is negative and increasing for a large. Therefore, ˆf0 and ˆg0 satisfy condition (v) of Theorem 2.1.

    Having verified these conditions, we can begin to solve (2.1) and (2.2) subject to the boundary conditions (2.4) and (2.6), with ˆf0 and ˆg0 in place of f0 and g0, and ˆfn and ˆgn in place of fn and gn. For simplicity, we may change to time variable to τ=tT, so that the initial data corresponds to τ=0. The characteristic equations are unchanged. The only change is the jump discontinuity in ˆf0 and ˆg0. Note that jump discontinuities do not impact the continuity of the expressions in (2.13) and (2.14). Indeed, the proof that these expressions are continuous assumed a jump discontinuity, at a=t. In our new variables, that discontinuity is at a=τ+T; ˆgn and ˆfn are, in fact, continuous at a=τ. Indeed, since

    f(0,T)=0βg(α)g(α,T)dα

    and

    g(0,T)=20βf(a)f(α,T)dα

    ˆf0 and ˆg0 satisfy the boundary conditions

    ˆg0(0)=20βf(α)ˆf0(α)dα

    and

    ˆf0(0)=0βg(α)ˆg0(α)dα.

    Therefore, when a=τ we have,

    limaτˆgn(a,τ)=2eBg(0,τ)0βf(α)ˆf0(α)dα=eBg(0,τ)ˆg(0)=limaτ+ˆgn(a,τ) (6.66)

    Similarly,

    limaτˆfn(a,τ)=limaτ+ˆfn(a,τ) (6.67)

    Thus, ˆfn and ˆgn converge to continuously differentiable solutions of (2.7) and (2.8) subject to (2.3)–(2.6), on {(a,τ)|0τT,0aT+τ}. Returning to our original variables, we extend our original solution to t2T. Continuing in this way, for all time, we can define solutions of (2.7) and (2.8) subject to the boundary conditions (2.3)–(2.6), continuously differentiable for at.

    Figure 11.  Inverse Gaussian model with stable initial age density and initial stage ratio perturbed by a factor of three. Note that the perturbed initial stage structure alone is sufficient to induce oscillations in both age density and stage structure.
    Figure 12.  Exponential model with stable initial age density and initial stage ratio perturbed by a factor of three. Note that the stage structure quickly stabilizes and does not overshoot or oscillate.


    [1] P. Gabriel, S. P. Garbett, V. Quaranta, D. R. Tyson, G. F. Webb, The contribution of age structure to cell population responses to targeted therapeutics, J. Theor. Biol., 311 (2012), 19-27.
    [2] S. S. Hamed, R. M. Straubinger, W. J. Jusko, Pharmacodynamic modeling of cell cycle and apoptotic effects of gemcitabine on pancreatic adenocarcinoma cells, Cancer Chemother Pharmacol, 72 (2013), 553-563.
    [3] X. Miao, G. Koch, S. Ait-Oudhia, R. M. Straubinger, W. J. Jusko, Pharmacodynamic modeling of cell cycle effects for gemcitabine and trabected in combinations in pancreatic cancer cells, Front. Pharmacol., 7 (2016), 421.
    [4] J. A. Alberts Bruce, Molecular Biology of the Cell: a Problems Approach, 4th edition, Garland Science, New York, NY, 2002.
    [5] E. S. Wenzel, A. T. Singh, Cell-cycle checkpoints and aneuploidy on the path to cancer, In Vivo, 32 (2018), 1-5.
    [6] A. Zetterberg, O. Larsson, K. G. Wiman, What is the restriction point?, Curr. Opin. Cell Biol., 7 (1995), 835-842.
    [7] C. Schwarz, A. Johnson, M. K. oivomägi, E. Zatulovskiy, C. J. Kravitz, A. Doncic, et al., A precise Cdk activity threshold determines passage through the restriction point, Mol. Cell, 69 (2018), 253-264.
    [8] A. Zilman, V. Ganusov, A. Perelson, Stochastic models of lymphocyte proliferation and death, PLoS One, 5 (2010), e12775.
    [9] A. V. Gett, P. D. Hodgkin, A cellular calculus for signal integration by T cells, Nature, 1 (2000), 239-244.
    [10] K. Léon, J. Faro, J. Caneiro, A general mathematical framework to model generation structure in a population of asynchronously dividing cells, J. Theor. Biol., 229 (2004), 455-476.
    [11] R. Callard, P. Hodgkin, Modeling T- and B-cell growth and differentiation, Immunol. Rev., 216 (2007), 119-129.
    [12] J. A. Smith, L. Martin, Do cells cycle?, Proc. Natl. Acad. Sci. U.S.A., 70 (1973), 1263-1267.
    [13] A. Golubev, Exponentially modified Gaussian (EMG) relevance to distributions related to cell proliferation and differentiation, J. Theor. Biol., 262 (2010), 257-266.
    [14] A. Golubev, Genes at work in random bouts, BioEssays, 34 (2012), 311-319.
    [15] A. Golubev, Applications and implications of the exponentially modified gamma distribution as a model for time variabilities related to cell proliferation and gene expression, J. Theor. Biol., 393 (2016), 203-217.
    [16] S. J. Cain, P. C. Chau, Transition probability cell cycle model part I-balanced growth, J. Theor. Biol., 185 (1997), 55-67.
    [17] S. Svetina, B. Žekš, Transition probability model of the cell cycle exhibiting the age-distribution for cells in the indeterministic state of the cell cycle, in Biomathematics and Cell Kinetics (eds. A. J. Valleron and P. D. M. MacDonald), Elsevier/North-Holland Biomedical Press, New York, 1978, 71-82.
    [18] S. Cooper, The continuum model: statistical implications, J. Theor. Biol., 94 (1982), 783-800.
    [19] S. Banerjee, K. Lo, M. K. Daddysman, A. Selewa, T. Kuntz, A. R. Dinner, et al., Biphasic growth dynamics during Caulobacter crescentus division, bioRxiv, (2017), 047589.
    [20] Z. W. Jones, R. Leander, V. Quaranta, L. A. Harris, D. R. Tyson, A drift-diffusion checkpoint model predicts a highly variable and growth-factor-sensitive portion of the cell cycle g1 phase, PLOS ONE, 13 (2018), 1-20.
    [21] J. Folks, R. S. Chikara, The inverse Gaussian distribution and its statistical application-a review, J. R. Statist. Soc. B, 40 (1978), 263-289.
    [22] O. Angulo, J. López-Marcos, Numerical integration of nonlinear size-structured population equations, Ecol. Model., 133 (2000), 3-14.
    [23] S. C. Tate, S. Cai, R. T. Ajamie, T. Burke, R. P. Beckmann, E. M. Chan, et al., Semi-mechanistic pharmacokinetic/pharmacodynamic modeling of the antitumor activity of LY2835219, a new cyclin-dependent kinase 4/6 inhibitor, in mice bearing human tumor xenografts, Clin. Cancer Res., 20 (2014), 3763-3774.
    [24] R. R. Goldberg, Methods of Real Analysis, John Wiley & Sons, Hoboken, NJ, 1976.
    [25] K. Ito, F. Kappel, G. Peichel, A fully discretized approximation scheme for size-structured population models, SIAM J. Numer. Anal., 28 (1991), 923-954.
    [26] C. C. Pugh, Undergraduate texts in mathematics, in Real Mathematical Analysis (eds. S. Axler, F. Gehring and K. Ribet), vol. 19, Springer International Publishing, New York, NY, 2003.
    [27] Wikipedia, Leibniz integral rule, 2019. Available from: https://en.wikipedia.org/wiki/Leibniz_integral_rule.
  • This article has been cited by:

    1. Olga Dolgusheva, Olena Pozharytska, Tetiana Yefymenko, Yana Prosiannikova, Anna Pogorila, Effectiveness of linguistic and extralinguistic tools across English language comics during the perception of scientific information, 2023, 26, 22232621, 7, 10.5782/.kjhss.2023.7.20
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(4172) PDF downloads(343) Cited by(1)

Figures and Tables

Figures(12)  /  Tables(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog