DEMOGRAPHIC MODELING OF TRANSIENT AMPLIFYING CELL POPULATION GROWTH

Quantitative measurement for the timings of cell division and death with the application of mathematical models is a standard way to estimate kinetic parameters of cellular proliferation. On the basis of label-based measurement data, several quantitative mathematical models describing shortterm dynamics of transient cellular proliferation have been proposed and extensively studied. In the present paper, we show that existing mathematical models for cell population growth can be reformulated as a specific case of generation progression models, a variant of parity progression models developed in mathematical demography. Generation progression ratio (GPR) is defined for a generation progression model as an expected ratio of population increase or decrease via cell division. We also apply a stochastic simulation algorithm which is capable of representing the population growth dynamics of transient amplifying cells for various inter-event time distributions of cell division and death. Demographic modeling and the application of stochastic simulation algorithm presented here can be used as a unified platform to systematically investigate the short term dynamics of cell population growth.

1. Introduction.The kinetic balance between inflow and outflow of cells contributes to maintain tissue homeostasis.A stem cell possesses self-renewal ability to produce a copy (copies) of itself by self-renewal cell division for some cases [29].By differentiating cell division, a stem cell population generates transient amplifying progenitor cell(s) [30].Although progenitor cells are generally capable of undergoing at most a finite number of cell divisions, for epidermal cells as an example, the mass population size is maintained by active and rapid proliferation of progenitor cells [3].Cellular proliferation can occur if a cell receives stimulus such as growth factors and stress.For example, naive T cells are activated and then proliferate after receiving appropriate signals via T cell receptors from antigen presenting cells [32].
At a homeostatic state, the number of T lymphocytes is maintained at a constant level in adult individuals.This fact indicates that the balance of cell division and death should be tightly regulated.Note that a T lymphocyte population is composed of several heterogeneous sub-populations (effector, naive and memory) in terms of biological function and cellular kinetics.One of most important biological questions to address is how the balance of the population number is constantly maintained despite the population is heterogeneous.To address this question, tracking experiments of T cells using quantitative labeling technique have been conducted to reveal cellular kinetics [7].Due to the development of label-based measurement methods such as using Carboxy Fluorescein Succinimidyl Ester (in short, CFSE), it is possible to track the time series process of cell division and death in a quantitative manner.The number of cell divisions that a cell has experienced can be measured by fluorescent intensity of CFSE.
Time-series quantitative measurement of cellular proliferation in combination with the use of mathematical models enables us to estimate kinetic parameters of cell population growth dynamics.By constructing mathematical models, quantitative data representing cell population growth can be used to estimate growth kinetic parameters of lymphocytes (the rate for cell division and death).Several different types of mathematical models have been proposed and extensively investigated; formulation by way of stochastic birth death process [43,28], partial differential equations [25,24,1], a simple model [22], statistical model [12].See [6] for a comprehensive review on the study for estimating growth kinetic parameters of lymphocytes by using CFSE.Most of previous mathematical studies regarding CFSE labeling experiments provide the estimation of kinetic parameters based on curve fitting with quantitative time series data representing the change of population number.These estimated parameter values provide useful information since it succeeded in quantitatively describing an effector phase of lymphocyte proliferation: the cell number increases due to clonal expansion, and then decreases due to elevated apoptosis rate as the number of replications increases.Because mathematical models used for parameter fitting differ among previous studies, comparison of results is made by careful biological interpretations.There exist possibilities that mathematical study can further contribute to provide quantitative description of T lymphocyte population dynamics especially in the following two issues.
One is to provide a flexible framework which allows to consider several different types of cell division and death in a unified manner.In order to describe the growth dynamics of a cell population, it is necessary to specify underlying stochastic rules for inter-event time distributions of cell division and cell death.Most of mathematical models considered in previous study are formulated by ordinary differential equations.The formulation by way of ordinary differential equations postulates that inter-event time distributions of cell division and cell death follow exponential distributions.However, this implicit assumption is not always appropriate to represent the process of cellular proliferation.In fact, several quantitative measurement experiments based on the tracking of cell division and cell death at the single cell level suggest that inter-event time distributions of cell division can be well characterized by long-tailed distributions such as log-normal distributions rather than exponential distributions.Moreover, inter-event time distributions of cell death are well approximated by Weibull distributions [15].This implies that death rate of cells is not constant as postulated in the formulation of ordinary differential equations or Poisson counting process.These observations suggest that an internal state change may affect the dynamics of population growth.
Second is a proposal of useful quantity which quantitatively captures growth dynamics of transient amplifying cell populations.Estimated parameters by themselves do not explicitly provide information of whether a cell population is growing or not: population dynamics is determined by the balance between the increase by cell division and loss by cell death.A quantity which represents the ratio between successive two divisions can be useful to quantitatively capture the dynamics of a growing or declining cell population.
To address the two issues above, in the present paper we would like to stress that each existing quantitative mathematical model describing transient amplifying cell population growth is always reformulated as a variant of the parity progression model.The parity progression model was originally proposed in human demography to describe the birth interval dynamics in human fertility [9,18,17].The parity progression model has been extensively used to estimate various demographic indicators in practice.In the context of cell demography, we newly define the generation progression ratio (GPR) as an alternative concept to the well known basic reproduction number R 0 , which provides a quantitative measure for the expected number of daughter cells generated from their mother cell 1 .
The organization of this paper is as follows.In Section 2, we show that existing quantitative mathematical models can be all seen as a specific case of the generation progression model.We define the generation progression ratio and calculate it for each specific cell population growth model.As an application, in Section 2.6, GPR is calculated for some experiment on the single cell level cell tracking measurement [15].Finally, we make use of an algorithm developed in [33] to implement stochastic simulations to represent the short term population growth dynamics of transient amplifying cells.In Appendix, we show the relation between R 0 and GPR by formulating the generation progression model as a multi-group McKendrick equation system.

2.
Transient amplifying cell population growth models.Throughout this paper, we consider population dynamics of transient amplifying cells.A schematic representation of cell division and death is depicted in Figure 1.2.1.Generation progression model.We describe the growth dynamics of a transient amplifying cell population by a system of linear renewal equations.In this paper, we assume that the number of divisions that a cell has experienced defines the generation.The generation is a group of cells which have experienced the same number of cell divisions.Transient amplifying cells are generally capable of undergoing a finite number of cell divisions.Hence it is reasonable to assume that there exist totally N + 1 distinct generations (j = 0, 1, 2, ..., N ).Let us denote the state of source-type cells such as a stem cell population as the 0-th generation.We distinguish the 0-th generation from the other generations because a stem cell population or unstimulated resting cell population such as naive T cells can have different proliferative kinetics from that of transient amplifying cells.Transition from j-th generation to (j + 1)-th generation occurs if a cell undergoes cell division (j = 0, 1, ..., N − 1).It follows from the definition of the generation that there is no transition from the j-th generation to the other generations other than j + 1 (j = 0, 1, 2, ..., N − 1).Although undivided transient amplifying progenitor cells can be generated from a source type cells, in the present paper, we always assume that there is no external recruitment of cells.
Let τ denotes the duration of a cell in a particular generation elapsed from the time that the cell enters the current generation, so simply τ is called "age" of a cell.Let λ j (τ ) denote the force of cell division which represents the rate of cell division.Incidence probability of cell division Λ j (τ ) is defined as a probability of undergoing cell division from j-th generation to (j + 1)-th generation with the rate λ j (τ ) during infinitesimal period (τ, τ + ∆τ ).Λ j (τ ) is mathematically formulated as Assume that cells at j-th generation are removed due to cell death with the rate µ j (τ ).Let F j (τ ) denote the survival probability of cells at j-th generation (j = 1, 2, ..., N ): Let b n (t) denote the cell population production rate (n = 0, 1, ..., N ) at time t which is an equivalent notation to the population birth rate in mathematical demography and epidemiology.For n = 1, 2, ..., N , the explicit form of generation progression model is given by (2.3) where b 0 is a given initial data.System (2.3) has already been proposed in [4] and investigated in [44] in details 2 .
Hereafter we focus on different aspects of system (2.3) that have not been fully investigated in the previous papers yet.Note that system (2.3) takes a quite similar form of the parity progression model [18,17,9].The difference between the parity progression models and the generation progression model (2.3) for j = 1, 2, ..., N is that factor "2" representing the number of daughter cells produced from a single mother cell is multiplied in each equation of (2.3).
We define the generation progression ratio (GPR) at n-th generation which describes the expected number of daughter cells produced from a single mother cell at n-th generation under given incidence probability of cell division Λ n (τ ) and survival probability F n (τ ).The generation progression ratio is mathematically formulated as (2.4) By the definitions of Λ n (τ ) and F n (τ ), we find that which shows that the generation progression ratio is twice the probability that a cell survive the n-th generation.Moreover, generation progression ratio (2.4) represents the ratio between one generation to the next generation.To see this aspect, define the size of generation P n by (2.6) By integrating the both side of (2.3) from 0 to ∞ with respect to t and changing the order of integration, we obtain that which shows that GPR gives the ratio of size of successive generations.It follows from (2.7) that the population size would increase on average if GPR(n) > 1, whereas it would decrease on average if GPR(n) < 1. Hence the generation progression ratio determines a threshold of expansion or shrink of population size between generations.
However, it should be noted that the basic reproduction number R 0 for total cell population composed of N + 1 distinct generations is zero if N is finite (see Appendix), so GPR is an index of reproductivity of each generation in the transient phase.On the other hand, if N = ∞, lim n→∞ GPR(n) is, if it exists, equivalent to the basic reproduction number R 0 for total cell population.This fact is a particular characteristics to transient amplifying cell populations because generation in the present context is prescribed by the number of cell divisions that a cell has experienced, and each "generation" is, in general, different from each other, so the generation progression process is not necessarily a self-renewal process.On the other hand, if we consider population dynamics of a stem cell population, reproduction of a stem cell occurs via self-renewal.The basic reproduction number R 0 of stem cell populations is therefore defined as an expected number of daughter cells still having self-renewal capability which are reproduced from a single stem cell during its survival period and can be identified with mother cells.In order to distinguish the basic reproduction number and the generation progression ratio, we remain to use the notation of generation progression ratio to represent the expected fraction of population expansion or shrink between two generations.
If both Λ n (τ ) and F n (τ ) are independent of generation number n, then generation progression ratio (2.4) becomes independent of n.In the following subsections, we do not explicitly write the dependence of n if it is not necessary.

Synchronous progression model.
Assume that all cells divide synchronously every time period T .This implies that the incidence probability of cell division is given by the Dirac's delta function ( Then the model for cell population growth is given by (2.9) Note that system (2.9) defines a system of difference equations for each fixed time t.
We derive an explicit solution of (2.9) for a special case.Assume that the survival probability of cells is independent of n and has the form of an exponential function with rate parameter µ: (2.10) Consider the following special case for 0-th generation: where C represents the initial population size of 0-th generation.Then we obtain the following explicit solution of (2.9): (2.12) The GPR for (2.9) is calculated as (2.13) The explicit solution (2.12) indicates that the population number exactly increases or decreases with the generation progression ratio 2e −µT at every cell division (see also the right panel of Figure 2).However, it is not reasonable to assume that all cells synchronously undergo cell division.In the following subsections, we consider different types of probability distribution for the inter-event time of cell division and death.

2.3.
Single step constant progression model.Assume that inter-event time distribution of cell division and death follow an exponential distribution with rate parameter λ and µ, respectively.Then incidence probability of cell division Λ n (τ ) and survival probability F n (τ ) at n-th generation are explicitly given by and For 0-th generation, we assume that survival probability at 0-th generation is given by an exponential function with rate parameter µ + λ.Let C denote the initial population size.Then system (2.3) is rewritten as By changing variable τ to t − u and differentiating the both sides of (2.16) with respect to t, the following ordinary differential equations are derived from (2.16): (2.17) System (2.17) corresponds to the mathematical model proposed in Revy et al. [35,4].Since cell division or cell death as an event occurs at a constant rate and the event is represented by a single process, we call system (2.16) as single step constant progression model.The generation progression ratio for single step constant progression model (2.16) is calculated as (2.18) The term λ λ+µ in (2.18) represents the conditional probability that a cell survives in a current generation and then transits to the next generation by cell division (see also [41, page 11]).Hence the generation progression ratio (2.18) represents the expected number of daughter cells reproduced from a mother cell.
2.4.Multi-step constant progression model.Since cell division process is composed of sequential multiple sub-processes, it is reasonable to assume that cell division process is further distinguished into several sequential steps.Assume that cell division process can be described by m distinct sequential sub-processes.A natural candidate to represent such inter-event time distribution of cell division is a gamma distribution.The Gamma distribution with rate parameter λ and shape parameter k (k = 1, 2, ..., m) is given by An exponential distribution with rate parameter λ is derived from (2.19) as a special case k = 1.We assume that the incidence probability of cell division and the survival probability of cells at n-th generation are respectively given by and (2.21)For n ≥ 1 and k = 1, 2, ..., m − 1, let x n|k (t) denote the density of cells at n-th generation in the k-th intermediate sub-process.The explicit form of x n|k (t) is given by (2.22) For m = 1, we set x n|0 (t) = b n (t).In the similar manner to the linear chain trick method [26], system (2.3) is reduced to the following system of ordinary differential equations (k = 1, 2, ..., m − 1): For m = 1, we obtain single step constant progression model (2.17).
The GPR for multi-step constant progression model (2.23) is given by Finally we show that synchronous progression model (2.9) defined in subsection 2.2 can be derived as a limiting case of multi-step constant progression model (2.23).More precisely, we consider a situation that there exist many sequential subprocesses in cell division process, and the rate of transition from one sub-process to another occurs rapidly.This situation can be represented as sufficiently large m and λ in the Gamma-distribution (see also the left panel of Figure 4).Assume that λ/m is constant.Let T denote the ratio of m to λ: λ = m/T .By applying the method developed in [42,Appendix], we can show that as m → ∞.Correspondingly, the GPR for (2.23) converges to the GPR for (2.9) as m → ∞.In fact, by substituting λ = m/T into (2.24) and taking the limit m → ∞, we find that Although synchronous progression model (2.9) seems to be less realistic, it can be derived from multi-step constant progression model (2.23) as a limiting case provided that cell division process can be represented by fast and many sequential sub-processes.The similar idea can be applied to another process in which a fixed time interval as time delay is incorporated to represent an inter-event time (see the next subsection).
The survival probability of cells at n-th generation is given by an exponential function with rate parameter µ: (2.28) Let A n (t) and B n−1 (t) denote the cell density at n-th generation in the A-phase and B-phase, respectively.Note that A n (t) and B n−1 (t) forms a pair (cf., [44]).For 0-th generation, we assume that b 0 (t) is given by an exponential function with rate parameter λ + µ.Then system (2.3) is deduced to We have replaced b n (t) with A n (t) to keep consistency of notation.The subordinate renewal equation for B n (t) is defined as (2.31) If T = 0, then Smith-Martin type model (2.31) is reduced to ordinary differential equations.In particular, variable A n (t) corresponds to single step constant progression model (2.17).
The GPR for Smith-Martin type model (2.29) is given by (2.32) Note that the GPR for the Smith-Martin type model corresponds to the GPR for single step constant progression model (2.18) multiplied with the survival probability of cells in B-phase e −µT .As shown in subsection 2.4, discrete-type time delay T can be derived from multistep constant progression model (2.23) as a limiting case if cell division process can be represented as fast and many sequential sub-processes.Hence the Smith-Martin type model can be approximated by multi-step constant progression model (2.23).A major advantage of using multi-step constant progression model (2.23) is to flexibly bridge two typical formulations for incidence probability of cell division: exponential distribution and fixed time interval represented by a discrete-type time delay.2.6.Hawkins' quantitative measurements.In this subsection, we reproduce the growth dynamics of a cell population studied in [15].In subsections 2.2-2.5, it is assumed that survival probabilities are described by exponential functions.In single cell measurement experiments performed in Hawkins et al., the timings of cell division and death are traced at the single cell level.Tracking data tor the timings of cell division and death suggest that incidence probability of cell division is well fitted by a log-normal distribution or Gamma distribution, whereas the survival probability of the cell population is better fitted by a Weibull distribution or lognormal distribution rather than exponential distributions.These observations imply that the time required for accomplishing cell division considerably varies.Moreover, age-dependent cell death would be more feasible than age-independent (constant rate) cell death to describe death process of the population studied in [15].
Since incidence probability of cell division and survival probability have already been determined in [15] by means of a model selection method, we directly substitute them into generation progression model (2.3) with a slight modification.More precisely, we assume that incidence probability of cell division at n-th generation Λ n (τ ) is given by a log-normal distribution with average µ n and variance σ 2 n of the form:

.33)
According to the single cell level measurement experiments [15], the mean and variance for the first-time cell division are given by µ 0 = 37.21 and σ 0 = 3.76, respectively.On the other hand, the mean and variance for more than second cell divisions are given by µ n = 9.3 and σ n = 2.54 (n ≥ 1), respectively.We assume that survival probability of cells at n-th generation F n (τ ) is given by the Weibull distribution of the form: The explicit values of shape and rate parameters m and 1/η are given by m = 1.86 and η = 33.02,respectively.Note that in Figure 1I of the original paper [15], survival probability curves for each generation are drawn.In Figure 1I Since we do not know whether the right hand side of integration in (2.35) has an analytic solution, we numerically compute the value of the GPR for totally eight generations n = 1, 2, ..., 8 (see Table 2.6).The integration is solved by the double exponential integration method [40].The numerical computation results indicate that the number of the population would increase on average until the fourth generation.While the population number would decrease after the fifth generation due to elevated death rate.We perform linear regression analysis to estimate linear relationship between GPRs and generations.Let x and y denote the continuous variable representing generations and GPRs, respectively.The data in Table 2. Since the experimental data is provided as a sum of cells in all generations, we do not know when a cell population starts shrinking.According to the calculated values of GPRs in Table 2.6, the cell population begins to decline from the fifth division, since the calculated value of GPR for the fifth generation becomes less than 1.In the experimental condition of [15], calculated values of GPRs suggest that death of cell known as activation induced cell death may occur within five times of replications.Since activation induced cell death in B cells is determined by the complex interplay of cell survival and apoptosis signals [8], outcome can change depending on experimental conditions.Calculation of GPRs can be integrated to a mathematic-assisted practical in vitro assay which enables us to evaluate the strength of BCR signaling, or the transcription level of BCL-2, a well known transcription factor that acts as an anti-apoptotic effects to B cells.

Stochastic simulation.
If time-series quantitative data such CFSE fluorescence intensity are available, then it becomes possible to estimate kinetic parameter values of cellular proliferation by iteratively solve a system of differential equations under an appropriate optimization method.In subsections 2.3-2.5, system of renewal equations is reduced to (delay) differential equations.For ordinary and delay differential equations, several useful packages or softwares are available to approximately solve differential equations [39] 1.Generation progression ratios for the Hawkins' quantitative measurement data.integration methods to numerically solve renewal equations as a special class of the second type Volterra integral equations [5,31], to our knowledge, there are no practical packages which are directly applicable to generation progression model (2.3).Hence a systematic approach to estimate kinetic parameter values of cellular proliferation for a given system of (renewal) equations is difficult to undertake in practice.Moreover, it is important to take into account for demographic stochasticity in cell population growth if population size is small.
Recently, one of the authors developed an algorithm for implementing stochastic simulation, called the individual-based Gillespie algorithm [33].The individualbased Gillespie algorithm can be applied to a fairly general class of structured population models including age-structured population models.Because the Gillespie algorithm is a numerical representation of Poisson counting process [13], it is not directly applicable to general counting process such as age-dependent death process.The developed algorithm by the authors is capable of reproducing survival curve of the Weibull type distribution [33].In this subsection, we demonstrate that for each generation progression model studied in subsections 2.2-2.6, the developed algorithm is applied to perform stochastic simulations by specifying incidence probability of cell division Λ n (τ ) and survival probability F n (τ ).
To implement stochastic simulations for (2.3), the timings of cell division or death should be calculated based on a given probability distribution for the inter-event time distribution of cell division or cell death.The choice of an event (cell division or death) that will occur next and the time step ∆t to proceed are determined by the rule prescribed in the algorithm [33].Basically, pseudo random variable to determine the timing of cell division or death is calculated from a given probability density function for cell division or death.If a cell at n-th generation is chosen to undergo cell division, then two daughter cells appear at (n + 1)-th generation, and the number of cells at (n + 1)-th generation increases.While if a cell is chosen to undergo apoptosis, the cell disappears and the number of cell in the corresponding generation decreases.In each time-step procedure ∆t, every cell gets older.
On each left panel of each Figures 2-6, incidence probability of cell division Λ n (τ ) and survival probability F n (τ ) are drawn (see also Table 2.7 as a summary of model ingredients and GPRs).While on each right panel of Figures 2-6, we draw a simulation result for the corresponding specific generation progression model.In all numerical simulations, we set N = 8 to specify the maximum number of cell divisions.Figure 2 represents a simulation result for synchronous progression model (2.9) with the parameter values: C = 100 (initial population size), µ = 0.05 and T = 10.Since all cells divide synchronously every time period T , the solution of system (2.9) exhibits an artificial trajectory.Figure 3 represents a simulation result of single step constant progression model (2.16) with the parameter values: C = 100 , µ = 0.05 and λ = 0.1.In contrast to the result for synchronous progression model (2.9), no synchronous patterns are observed.Figure 4 represents a simulation result of multi-step constant progression model (2.3) with incidence probability of cell division (2.20) and survival probability (2.21).Parameter values are given as follows: C = 100, λ = 50 and m = 500.The numerical values of λ and m are set sufficiently large enough to find the tendency of convergence to synchronous progression model (2.9).In fact, Figure 4 exhibits a similar time-course pattern as drawn in Figure 2. Figure 5 represents a simulation result of Smith-Martin type model (2.29) with the parameter values: C = 100, µ = 0.025, λ = 0.1 and T = 10.Note that the timing of transition from A-phase to B-phase is determined by a random variable generated from an exponential distribution, while the transition from B-phase to A-phase occurs after fixed duration T has past.Finally, a simulation result of the generation progression model (2.3) with incidence probability of cell division (2.33) and survival probability (2.34) as adopted for Hawkins' quantitative measurement data are drawn on Figure 6.The parameter values written in subsection 2.6 are used to carry out numerical simulations.Although the assumption postulated in survival probability of cells (2.34) is necessary to reconsider, the simulation result shows qualitatively good agreement with experimental data (see Figure 1A of the original paper [15] in which cell population growth (blue solid line) and accumulated cell death count (red solid line) for 150 hours starting from 100 initial cells are drawn: Open access, available at the following website: http://www.pnas.org/content/106/32/13457.figures-only).
The difference among stochastic simulations for specific models in subsections 2.2-2.6 shows the importance of inherent stochasticity in a cell-cycle phase.If cellcycle is precisely regulated such as for the egg of xenopus laevis (African clawed frog) [34], the expected dynamics can be represented by synchronous model (Figure 2).On the other hand, if there is stochasticity in a cell cycle period, we observe qualitatively different (exponential) growth curves depicted in Figures 3, 5 and 6.Statistical representation of the incidence probability of cell division and death as proposed in this paper provides a useful and versatile framework to compare the effect of stochasticity that naturally exists in cell division and death.2. Comparison of incidence probability of cell division Λ n (τ ) and survival probability F n (τ ) for different generation progression models with the generation progression ratio (GPR).
3. Discussion.We showed that several existing quantitative mathematical models can be reformulated as a generation progression model, a variant of parity progression models developed in the field of mathematical demography.In Section 2, we focused on several typical transient cell population growth models which were proposed to estimate growth kinetics parameters for cellular proliferation of lymphocytes.By specifying incidence probability of cell division Λ n (τ ) and survival probability F n (τ ), we showed that existing models described by differential   equations are equivalently reformulated as a specific case of generation progression models.Hence it becomes possible to compare differences in terms of underlying implicit assumptions among existing mathematical models as the differences in terms of incidence probability of cell division and survival probability.Although these considerations have already been discussed in [44] in part, we further showed that single step constant progression model (2.17) and Smith-Martin type model (2.29) can be derived from multi-step constant progression model (2.23).Moreover, generation progression ratios (GPRs) were calculated for each generation.The generation progression ratio provides a quantitative estimate for the expected fraction of population increase or decrease between two generations.We applied the stochastic simulation method developed in [33] to numerically simulate the dynamics of cell population growth for each generation progression model.In implementing numerical computations to differential equations, a specific numerical integration method should be chosen depending on the type of differential equations (ODEs, DDEs or PDEs).While the stochastic simulation method presented in [33] is applicable to all generation progression models independent of the choice of a specific numerical integration method.For the quantitative data investigated in subsection 2.6, it could be more reasonable to use stochastic simulations since demographic stochasticity due to the low number of population size might be expected.Although we focused on formulating a systematic approach for construction of specific generation progression models and implementation of stochastic simulations, more detailed analysis on the quantitative time-series data at the single cell level [15] would expect to provide quantitative description for the population growth dynamics of transient amplifying cells.Moreover, it is important to derive other existing quantitative mathematical models from generation progression model (2.3) such as the model proposed by Deenick et al [7], the cyton model [16], and timedependent Smith-Martin type model (IL-2 consumption model) [11].
Age-structured cell population dynamics have been extensively studied during 1960s to 1980s [37,10,2,21,36,14].Age-structured population model was firstly proposed by McKendrick in 1926 to investigate the spread of infectious agents, known as the McKendrick equation [27] (see also [23], a comprehensive review paper for the brief summary of past existing study on McKendrick equations).
McKendrick equations and the Lotka renewal equation as the equivalent alternative formulation have been extensively studied in mathematical demography and epidemiology.Several important demographic indicators such as the basic reproduction number R 0 can be defined for the McKendrick equation.In Appendix, we discuss the relation between a multi-group McKendrick equation system and generation progression models.Moreover, the relation between the basic reproduction number and generation progression ratio is discussed.If the maximum number of cell division (generation) N is infinite, then it is shown that R 0 is obtained as a limit of GPR(n) as n → ∞.
All of quantitative mathematical models focused in the present paper do not consider recruitment of cells.For relatively long-term dynamics of cell population growth, cells would be recruited from a somatic stem cell population via differentiating division.Since collection of quantitative and time-series data for the population growth of stem cells is difficult, it is practically useful if one could estimate kinetic parameters for the population growth model of stem cells from quantitative time-series data of their progenitor cells.For example, estimation of growth kinetics parameters of a mathematical model for hematopoietic stem cell population growth dynamics in bone marrow from some type of white blood cells in peripheral blood or a specific tissue could be practically important.As an another example, we have formulated a mathematical model which will be used to estimate growth kinetic parameters of epidermal cell population under the steady state condition of three-dimensional skin culture systems [33].Generation progression model (2.3) can be extended to include a class of stem cell population as a source of progenitor population.Moreover, it is important to incorporate the effect of growth factors in cell division and death.In this case, extension of generation progression model (2.3) to time-dependent generation progression model is an appropriate choice to incorporate the effect of growth factors (cf., [44]).It is important in practice to define generation progression ratios under the influence of growth factors.
Appendix A. Relation between R 0 and GPR.In order to see the relation between R 0 and GPR, we show that generation progression model (2.3) can be formulated as a multi-group McKendrick equation system ( [18], [20]).As mentioned in Section 2, each generation is defined as a cell group which is produced through the same number of cell divisions.
Let n(t, τ ) = (n 0 , n 1 , • • •n N ) T3 denote an N + 1-dimensional vector whose jth element represents an age density function of transient amplifying cells of j-th generation at time t.Let matrix Q be a (N + 1) × (N + 1) diagonal transition matrix as: , where λ N = 0 if the N -th generation does not reproduce daughter cells anymore.Moreover let M be a reproduction matrix as .
Then corresponding to our basic model (2.where Ψ(τ ) = M (τ )U (τ, 0) is explicitly given as As is well known in the theory of Volterra integral equation, the solution of (A.4) is given by the generation expansion: where each "generation" G j is given by the iteration process It is easy to see that G j is a vector whose elements are zero except for (j + 1)-th element and the (j + 1)-th element is given by b j (t).That is, the iteration (A.5) is no other than our generation progression model (2.3).
Then the next generation matrix (next generation operator if N is infinite) of the multi-group model is given by which is the generation expansion of the total cell population density.If N is the maximum of parity, G j = 0 and n j = 0 for j ≥ N + 1.Even in case that N = ∞, the generation expansion is a finite sum as long as reproductive period of a cell is finite.
Finally we consider the case that N = ∞.Then (A.4) becomes a 1 -valued abstract integral equation, so we need a careful treatment based on functional analysis.On the other hand, if λ j = λ and µ j = µ, that is, there is no difference among generations with respect to their reproduction and survival ability, we can reduce the infinite-dimensional problem into a scalar problem.In fact, the total density of cells p(t, τ ) satisfies the Lotka-McKendrick system ( [36]): Then in the linear phase, we can calculate the initial growth rate of the cell population based on the knowledge of the reproduction kernel ψ.
For more general case, if N = ∞, the basic reproduction number is the asymptotic per generation growth factor, so it would be obtained by lim n→∞ GPR(n) if it exists ( [19]).

Figure 1 .
Figure 1.A schematic representation for population dynamics of transient amplifying cells.There exists a maximum number of cell division.A cell undergoes cell division to enter a next generation.A cell will disappear due to cell death.See the main text for detailed definition and explanation.
6 are used to estimate the coefficients a and b of regression line y = ax + b.Estimated values for a and b are given by a = −0.2881and b = 2.2538.The slope a approximately represents the stepwise decline of proliferative ability of the cell population across generations.
[38,4,11]-Martin type model.Single step constant progression model(2.16)doesnotaccountfor the minimum fixed time required for accomplishing cell division.Therefore quantitative estimation for growth kinetic parameters of cell division by means of the single step constant progression model would become imprecise in particular for slowly dividing cells.The mathematical model proposed by Smith and Martin incorporates a fixed time interval required for accomplishing cell division[38,4,11].In the Smith-Martin type model, a cell division period is distinguished into two phases, A-phase and B-phase.A-phase approximately corresponds to G 1phase, whereas B-phase corresponds to the concatenation of S-phase, G 2 -phase and M-phase of cell division.It is assumed that the transition from A-phase to B-phase is a stochastic event.More precisely, it is assumed that the expected time to exit Aphase and enter to B-phase follows an exponential distribution with rate parameter λ.On the other hand, accomplishment of B-phase requires a fixed time T .Let us derive the Smith-Martin type model from(2.3).Incidence probability of cell division at n-th generation Λ n (τ ) is given by