Counting generations in birth and death processes with competing Erlang and exponential waiting times.

Lymphocyte populations, stimulated in vitro or in vivo, grow as cells divide. Stochastic models are appropriate because some cells undergo multiple rounds of division, some die, and others of the same type in the same conditions do not divide at all. If individual cells behave independently, then each cell can be imagined as sampling from a probability density of times to division and death. The exponential density is the most mathematically and computationally convenient choice. It has the advantage of satisfying the memoryless property, consistent with a Markov process, but it overestimates the probability of short division times. With the aim of preserving the advantages of a Markovian framework while improving the representation of experimentally-observed division times, we consider a multi-stage model of cellular division and death. We use Erlang-distributed (or, more generally, phase-type distributed) times to division, and exponentially distributed times to death. We classify cells into generations, using the rule that the daughters of cells in generation n are in generation [Formula: see text]. In some circumstances, our representation is equivalent to established models of lymphocyte dynamics. We find the growth rate of the cell population by calculating the proportions of cells by stage and generation. The exponent describing the late-time cell population growth, and the criterion for extinction of the population, differs from what would be expected if N steps with rate [Formula: see text] were equivalent to a single step of rate [Formula: see text]. We link with a published experimental dataset, where cell counts were reported after T cells were transferred to lymphopenic mice, using Approximate Bayesian Computation. In the comparison, the death rate is assumed to be proportional to the generation and the Erlang time to division for generation 0 is allowed to differ from that of subsequent generations. The multi-stage representation is preferred to a simple exponential in posterior distributions, and the mean time to first division is estimated to be longer than the mean time to subsequent divisions.

Introduction

Cells of the immune system patrol our bodies for months or years[1,2]. During an adaptive immune response, a subset of specific cells, initially a small fraction of the total population, expands as cells undergo multiple rounds of division over a few days[3]. Although most of these cells die as the infection is overcome, lasting immunity is ensured by the transformation, or “differentiation” of individual cells to a memory phenotype. The most convenient mathematical and computational models of the dynamics of cell populations, which can consider heterogeneity at the single-cell level, are Markov models. In these models, the variables describe the number of cells of each type as a function of time, and cellular events such as division, death or differentiation are defined by their associated rates; each event corresponds to a possible fate of an individual cell and cells are independent of each other. In this formulation, inter-event times are exponentially-distributed random variables, with probability density maximised at zero.

The rapid expansion of cohorts of lymphocytes is recreated in laboratories, either by stimulation in vitro or by transferring cells to lymphopenic mice. By labelling cells with carboxyfluorescein succinimidyl ester (CFSE) or cell trace violet at the beginning of an experiment, and then using flow cytometry at a later time, a cohort of cells can be classified into generations because each round of division dilutes the intracellular dye[4-9]. At the beginning of an experiment, when cells are labelled, it is considered that all cells are in generation 0. In time-lapse microscopy experiments, individual cells are tracked and correlations within family trees identified[10-17]. Over the timescales of such experiments, hours to days, it is not appropriate to treat cell division as an instantaneous event. Rather, cells are “cycling” through gap, synthesis and mitosis phases (/, S and M), and daughter cells cannot immediately redivide[18]. To improve on the exponential distribution, Smith and Martin proposed a model in which the time between divisions is the sum of a fixed time spent in phase B, corresponding to S//M, and a variable time spent in phase A, corresponding to [19]. In the “single stochastic division” model of Hogan et al., the rate of transition from A to B phase depends on the T cell clonotype and on the number of cells competing for the same resources[20]. If there is a common molecular mechanism controlling the time spent in all phases of the cell cycle, then phase B may, instead, occupy a fixed proportion of the total time[10]. Takahashi[21,22] divided the cell cycle into four phases, with the duration of each drawn from a Pearson type III distribution. Weber et al. postulated a delayed exponential waiting time for each of three phases, corresponding to , S and /M[23].

The cyton model is a general framework for modelling proliferating lymphocytes, based on the idea that each cell has a set of competing clocks, determining its fate[13]. A number of features are incorporated in the model: generation-dependent parameters, heritable factors, and correlations between cells of the same generation[11,16,24]. When the fate of an individual cell is determined by competing internal clocks, the probability density of observed times between divisions is not the same as that of the division clock because division only happens if another fate does not. The probability density of division times is said to be “censored”[25]. When all clocks have exponential probability densities, the probability density of observed division times is also exponential due to the memoryless property of exponential random variables; this conservation of shape does not hold for non-exponential distributions, including the log-normal and Erlang distributions.

Kendall[26] introduced the idea of cell division occurring at the end of a sequence of k phases, with an exponentially-distributed time spent in each phase. Luzyanina et al.[27] made use of a similar model, where each phase identifies a cell generation, to study methodological and computational issues related to CFSE labelling data. Gamma and Erlang distributions have been considered to model a cell’s time to division[28,29]. Yates et al.[30] incorporated the idea of a sequence of exponentially-distributed phases, also referred as stages, before division in a multi-stage representation of a population of dividing cells. This yields an Erlang distribution of times to division, while retaining some of the mathematical and computational advantages of the exponential distribution[30]. The number of stages and their mean duration can be used as free parameters to compare with experimental data[31,32]. On the other hand, the internal stages are a mathematical construct that do not directly correspond to biological phases.

In this manuscript we adopt a multi-stage representation of the cell cycle and we include cell death as a competing fate. In the simplest case, obtained by assuming identical birth and death rates across stages, we derive an analytical expression for the expected number of cells in each stage as a function of time, and study the limiting behaviour of the system as . We also consider cell generations, in order to make theoretical predictions comparable to CFSE experimental data[20]. The expected number of cells in each generation is found in the case of constant birth and death rates, and number of stages across generations. Although arising from different motivations, our multi-stage model with cell generations can be accommodated within the cyton framework, with the progressor fraction equal to one, division and death clocks following Erlang and exponential distributions, respectively. We show the applicability of this stochastic approach by calibrating our multi-stage model with cell generations and its exponential version with CFSE data from two populations of murine T cells[20]. Model calibration is performed making use of Approximate Bayesian Computation Sequential Monte Carlo (ABC-SMC) approaches[33].

In the “Methods: multi-stage models of cell division and death with a Markovian framework” section, we introduce the multi-stage model including cell death and obtain analytical results for the multi-stage framework under some conditions on the model parameters. In the “Case study: lymphopenia-induced proliferation” section, we calibrate the multi-stage model with cell generations by making use of CFSE data from Hogan et al.[20]. Its performance is compared with a simple exponential model of cell division. A final discussion is provided in the “Discussion” section.

Methods: multi-stage models of cell division and death with a Markovian framework

We present a multi-stage (MS) model of the time between cell divisions. Cells pass through a sequence of N stages before dividing. The stages are not directly related to the biological phases of the cellular cycle. The time to progress from stage j to the next one, , is an exponentially-distributed random variable with mean . We will refer to these rates, , , as birth rates. Times to death are also distributed exponentially, with per cell death rate . Thus, at each stage, each cell may either proceed to the next one, with probability , or die, with probability . The time increment is a random variable, following the exponential distribution with mean .

Figure 1 illustrates the dynamics. Our multi-stage model is equivalent to considering two independent clocks for cell division and death, which compete to decide the cellular fate. The time-to-death clock follows an exponential distribution with rate , while the division time follows a continuous phase-type distribution with parameters and [34]. A particular choice of phase-type distribution is the , which is a concatenation of N identically distributed exponential steps, where all birth rates are equal: , . The case has been considered by Yates et al.[30].

Figure 1

Multi-stage model of cell division and death (MS model). The cell cycle is divided into N different stages. A cell has to visit N stages in order to divide. At each stage j, , the cell may proceed to the next stage, with birth rate , or die, with death rate .

The number of cells in stage j at time t, is the random variable , . Let , be the expected value of . The following set of differential equations may be obtained by considering the events that can happen in a short time interval:When we extend the MS model to assign a generation to each cell, we refer to the model as the MS-G model. In this way, mean quantities can be compared with CFSE experimental data[8]. Histograms of CFSE intensity display a series of peaks, each corresponding to a generation, or a number of divisions over the course of the experiment[35].

In the MS-G model, generation is split into different stages. The notation reflects the fact that the number of stages may depend on the generation g. A cell in generation g has to sequentially visit all compartments to divide. On the other hand, cells might also die at any stage of the cycle. As depicted in Fig. 2, if a cell belongs to generation g and is in compartment j, , it may proceed to the following stage, with birth rate , or die with death rate . Again the notation reflects the potential for these rates to depend on the generation. When a cell reaches the last stage, , of generation g and divides, its two daughters will join the first compartment of generation . In summary, given a cell in generation g, its time to division follows an Erlang distribution with parameters , whereas its time to death follows an exponential distribution with rate . These distributions correspond to two independent competing clocks to control cellular fate, similarly to those considered in Fig. 1.

Figure 2

Multi-stage model with cell generations (MS-G model). Each cell in the first stage of generation 0 has to visit all the compartments (or stages) in order to divide. When cells arrive at the last stage of generation 0, , they may divide with birth rate , or die with death rate . If a cell divides, its daughter cells join the first compartment of the next generation, and the process continues.

The number of cells in stage j of generation g at time t is the random variable , , . Let be the expected value of . The following set of differential equations may be obtained by considering the events that can happen in a short time interval:We are interested in computing the mean number of cells over time for the MS and MS-G models. Specifically the MS-G model will provide the mean number of cells in each generation, and thus, can be used together with CFSE data to obtain division and death rates. When division times are Erlang distributed (MS model), or if one considers that those Erlang distributions are identical across generations (MS-G model), it is possible to carry out a comprehensive analytical study. This is shown in the “Analytical results” section.

When convenient analytical solutions cannot be obtained, (2.1)–(2.2) can be solved numerically in different ways. For example, for the MS-G model, and keeping in mind our interest in modelling CFSE data, we assume there exists a maximum generation G that can be measured by the dye. Thus, one might be interested in following cells within generations . For these generations, Eq. (2.2) can be solved by making use of the matrix exponential. To this end, let be the column vector of the mean number of cells in each stage and generation as time evolves, i.e.,which has length , and where the column sub-vectors contain the mean number of cells across stages in generations . Let us also define the coefficient matrixwhere is a square matrix, whereas is a matrix. is then a real square matrix of dimension , and represents a null matrix with dimension . Given the vector of the initial conditions , which has length , the system of Eq. (2.2) can be rewritten as the following Cauchy problemThe solution of the system is given by , whererepresents the matrix exponential. For efficient ways of computing this matrix, see Refs.[36-39]. Finally, we note that since CFSE data describe the number of cells in each generation, one can then compute the mean number of cells in each generation over time asUsing (2.2) and (2.3), we compute the time evolution of as followsIn a related approach by Zilman et al.[29], a cell’s time to division is a gamma-distributed random variable, and time to death is exponentially distributed. Solutions are given in terms of integral equations. Here, with Erlang-distributed division times, we find a set of linear differential equations for the expected number of cells in each stage.

Analytical results

In this Section, we show how the Markovian framework of the proposed multi-stage models provides analytical tractability under some simplifying assumptions. Our aim is to compute the mean number of cells in each stage and generation over time, especially the limiting behaviour as .

MS model with Erlang division time

In this Section, we consider a simple case of the MS model, where identical birth rates are assumed across different stages; that is, . The phase-type distribution for the time to division in Fig. 1 is and the mean time to division is given by . Note that when the MS model becomes a Markov linear birth-and-death process, with birth rate, , and death rate, . Equation (2.1) becomesAs in Yates et al.[30], we introduce the new variables , , which satisfy the following ODEs:We find an Nth-order homogeneous differential equation for that does not depend on :together with a set of ODEs that relate to the derivatives of with respect to timeThe solution is[30] and , where is the first Nth root of unity, and (for ) are constants which depend on the initial conditions.

If, at time , there are cells in the first stage and zero cells in any other stage, we haveThen the solutions of the system (3.1) in terms of the original variables areTherefore, the expected total number of cells in the population at time t, M(t), is given by

Population extinction or growth

Let us consider the MS model with Erlang division time, and define , the total number of cells in the population at time , so that . One can definethe probability of population extinction at late times, for m cells initially in stage j at time . Since cells behave independently, it is clear that , so we omit the superindex from now on and denote . A first-step argument leads toThese equations can be solved recursively, leading towhere . When , ; when , .

The analytical solutions (3.5)–(3.6) provide another route to study the limiting behaviour as . The largest term in the summation of (3.5) is the one corresponding to . The combination , is positive if , negative, when , or zero if . Thus, extinction of the cell population is certain if . Figure 3 (centre) shows an example of extinction when and the initial number of cells . An example of behaviour when is shown in the left panel of Fig. 3.

Figure 3

Limiting behaviour when of a population with an initial number of cells, . Birth and death rates, and , have units of inverse time, . Left: Parameters: , , . The population of cells in stage j levels out to for sufficiently large times. Centre: Parameters: , , . The population of cells at any stage becomes extinct at late times. Right: Parameters: , , . The populations grow according to (3.11) and the relation between and given by Eq. (3.9) is satisfied. For example, at , .

If then, as ,andwhich is illustrated in Fig. 3 (right). The total population size is easily obtained usingAs ,The exponent , describing the late-time cell population growth, plotted in Fig. 4, is lower than would be expected if N steps with rate were equivalent to a single step of rate . As , we have . In terms of (3.10), since as . Because the cell population is unevenly distributed across stages, with a bias towards earlier stages in the long run, N steps with rate are not equivalent to a single step of rate .

Figure 4

The exponent that determines the asymptotic growth rate of the population is shown against the number of stages. The dotted line would be expected if N steps with rate were equivalent to a single step of rate .

Mean fraction of cells at each stage

We define the mean fraction of cells in each stage, , as the ratio between the mean number of cells in compartment j and the expected total number of cells in the population, i.e.,We make use of (3.1) and (3.10) to write,which have the following steady state solutionOne observes that , , which means (on average) the fraction of cells decreases stage by stage, independently of the initial distribution of cells. In fact, one can solve (3.14) to determine , as followswhich does not depend on or . Thus, at late times the fraction of cells in each stage only depends on the number of stages considered; the parameter sets the timescale of the dynamical system, and all cells are equally susceptible to death, regardless of the stage they are in.

MS-G model with identical Erlang division times across generations

The solutions of the system (2.2) can be written in a closed analytical form in particular cases. For example, one may consider a simplified scenario where the number of stages is equal to 1 for all the generations, i.e., for all . Then, if we consider that at time , there are cells in generation 0, so that , this leads to the following solutions:In this case the MS-G model becomes a birth-and-death process tracking cell generations, and becomes identical to that considered in Refs.[27,40-42], where the inter-event times of cell death and division are modelled as exponential random variables, rather than Erlang distributions.

In this Section we consider the case with identical number of stages, N, and rates, and , for each generation, so that division times are Erlang-distributed in each generation. Under these assumptions, it is possible to obtain an analytical expression for the mean number of cells in each generation. Then Eq. (2.2) becomesThese equations can be rewritten in terms of the new variables , for , . This is equivalent to multiplying (3.17) by the integrating factor . Thus, (3.17) becomesTo determine the solutions of (3.18), we focus here on the case and all the other compartments are empty at time . This implies that for , and by solving (3.18) recursively one getsThis expression allows one then to determine the mean number of cells in each stage of generation 1,By recursion of (3.18) the mean number of cells in each compartment j of generation g is given byGoing back to the original variables, , the solutions of (3.17) areFrom the previous equations, one can show thatsince cells in each generation and compartment either proceed to the next stage within their generation, divide (proceeding to the next generation), or die.

Once the mean number of cells in each compartment for a given generation is at hand, the expected number of cells in each generation can be determined according to (2.3). We can writeThis equation is consistent with the results of the exponential model ()[27]. On the other hand, if one is interested in the mean number of cells in each compartment, for , regardless of the generation they belong to, this can be computed as followsfor and . In practice, one could truncate the series above to get an approximation of the mean number of cells in each stage. However, we note that one can use instead the solution provided by (3.5), since the dynamics of the MS-G model is equivalent to the dynamics of the MS model, when the parameters N, and are generation-independent. It can be numerically checked that this indeed provides equivalent results. In fact, when or , one can analytically show the equivalence. In the former case (), it is enough to recall the power series of the exponential function. In the latter case (), we derive from (3.5)where we used the fact that . On the other hand, from (3.22) we obtainThis shows that the two models lead to the same expected number of cells in each stage.

Comparison between the MS-G model and the cyton model

The cyton model is a stochastic model proposed to describe the population dynamics of B and T lymphocytes[13]. Division and death times are regulated by two independent clocks, and the competition between both clocks determines the fate of the cell. When a cell divides, these clocks, which depend on the number of divisions the cell has undergone, are reset for each daughter cell. However, when analysing an in vitro experiment with this type of cells, there is evidence that not all cells either divide or die. For instance, a portion of them may not respond to the stimulation[43], or may respond without division[44]. This is the reason why a progressor fraction is defined in the cyton model. This progressor fraction represents for a given generation, the fraction of cells that are capable of undergoing further division. Each clock is described by a probability density function, and the parameters that define these probabilities are the free parameters in the model. Right skewed distributions, such as log-normal or gamma, are usually adopted to characterise the two independent clocks that regulate cell division and death. In summary, the cyton model is based on the following assumptions:In order to translate these assumptions into mathematical terms, let be the progressor fraction characterising cells having undergone g divisions, and let and represent the probability density functions for the time to division and death, respectively, for cells in generation g. The number of cells dividing for the first time, or dying, per unit time at time can be calculated, respectively, as[13]:where is the initial number of cells in the population. Consequently, the time evolution of the expected number of cells in generation 0, , obeys the differential equationThe number of cells in generation g dividing, or dying, per unit time at time t can be computed, respectively, asHence, the dynamics of the average number of cells in each generation, , is governed by the differential equationsIn the next sections we show how the cyton model is equivalent to our model for particular choices of the probability density functions of the division and death clocks, and , and the progressor faction .

death and division are random events, characterised by a probability density function for the time to divide or die, respectively,

these processes are independent, and compete to determine the fate of the cell,

the clocks responsible for these processes are reset when a cell divides,

only a fraction of the cells in each generation are capable to undergo further divisions, and

the machineries that regulate cellular fate depend on the cell’s generation.

Exponential time to division and death

We consider here the MS-G model with number of stages across generations equal to one, i.e., for all . This means that cells in generation g divide after an exponentially distributed time with rate , and die with rate . We note that this is different to a standard Markov birth-and-death process, since rates are generation-dependent. Equation (2.4) becomeIn this case, our model is equivalent to the cyton model with exponential times for division and death, and progressor fraction , . One can show this equivalence by proving that and , by induction on g. In the cyton model, the assumption of exponential time to division and death implies that and , . Therefore, according to (3.23) and (3.24), the number of cells at time t dividing for the first time or dying to exit generation 0 per unit time is given byWe know from (3.16) that . Therefore, we can write and , which proves the case . We assume and hold for generation g and we need to show they also hold for generation . We make use of (3.16) and (3.26) to writeFor the number of cells in generation dying, Eq. (3.27), together with (3.16) lead towhich concludes the proof. With the identities and in (3.25) and (3.28), one can show that and obey the same differential equations for all . Thus, the two models are equivalent.

Erlang time to division and exponential time to death

We now consider the more interesting case where the number of stages in each generation is greater than one, and the cell cycle can be described as a multi-stage process. We focus here on the case where identical number of stages N and birth and death rates, and , respectively, are considered across generations. Similarly to the previous case, we prove that and by induction on g. Since a cell’s time to division is Erlang distributed and a cell’s time to death is exponentially distributed, for all andwhere the progressor fraction is again set to 1 for each generation. Note that in this case the parameters in and are independent of the generation g, since the number of stages and the birth and death rates are identical for all generations. From (3.23) and (3.24), the number of cells dividing for the first time or dying to exit generation 0 per unit time at time t isThe dynamics of the expected number of cells in generation 0 is given by (3.25), as in the previous case. Using (3.20) and (3.21), we observe thatTherefore, and , which concludes the case . We make use of these identities in (3.25) to obtainwhich is the differential equation derived in (2.4) for . Now, let us suppose that the identities and hold for generation g and we prove them for generation . Using (3.26) and the induction hypothesis, we havewhere we have used (3.20) for the last step. The same arguments can be used to look at the number of cells in generation dying per unit of time,  (3.27). Together with the induction hypothesis, we can writewhere the last identity was obtained making use of (3.21). Hence, (3.28) becomeswhich is identical to (2.4) for , . This concludes the proof of the equivalence between the cyton model and the multi-stage model with generations when a cell’s time to divide is Erlang distributed with parameters and N, and a cell’s time to die is exponential with rate . In summary, the analysis presented in this section for the multi-stage model with Erlang division time and exponential death time leads to exact closed solutions for the cyton model with the previous choice of clocks.

Case study: lymphopenia-induced proliferation

In this Section we illustrate the applicability of the MS-G model to CFSE data, making use of an experimental study of lymphopenia-induced proliferation[20]. In particular, we compare the performance of the MS-G model to that of a simple exponential (or single stage) model with generations, which is equivalent to making for all g in the MS-G model.

Differences in T cell proliferation have been observed to vary between different T cell clonotypes (i.e., the set of T cells with the same T cell receptor). Hogan et al.[20] transferred CFSE-labelled OT-I or F5 T cells intravenously to lymphopenic mice. A certain number of days (3, 4, 5, 6, 7, 10, 12 and 18 days) after the transfer, spleens and lymph nodes were recovered from the mice and analysed by flow cytometry to quantify the expression levels of CD8, CD5, CD44, and CFSE dilution[20]. For each time point, the number of mice analysed was between 3 and 7. We note that two independent transfer experiments, carried out under identical conditions, were performed: one for OT-I cells and a second one for F5. In Fig. 5 both data sets are shown: for each time point the number of cells is plotted for each mouse and generation (identified via the CFSE dilution measurement). On the left (right), OT-I (F5) cells are represented by the green (blue) histograms. In order to infer model parameters, we will consider all cells which have divided five or more times as a single class, denoted . This is similar to the approach considered in Refs.[4,29,35]. The rationale behind this choice is to reduce errors in the quantification of labelled cells with low CFSE fluorescence, as is the case for five or more divisions.

Figure 5

Data set of murine T lymphocytes from Hogan et al.[20]. Left: OT-I T cells. Right: F5 T cells. For each time point, the number of cells is plotted for each mouse and generation.

Figure 5 clearly shows that OT-I T cells proliferate faster than F5 cells, so that by day 7 there are OT-I cells in generation 10, whereas for F5 cells the maximum generation observed at day 7 is 6. This greater proliferative capacity of OT-I cells eventually leads, after one week, to competition for resources (e.g., IL-7 cytokine) and the OT-I population approaching its carrying capacity[20]. Since our model does not account for competition, it can only appropriately describe the dynamics of OT-I cells during the first week of the experiment. Thus, for OT-I cells we will only make use of the data set up to that time (one week). Yet for the F5 population we will use the entire data set. In Hogan et al.[20] this competition was incorporated with a density-dependent birth rate, , as followswhere is the rate of growth under unlimited resources, the size of reduction caused by the expansion of competing cells, and P is the size of the population[20]. Figure 6 shows the density-dependent birth rate, , as a function of the population size P. It suggests that the competition for resources is greater in the case of OT-I T cells. In the experiments the number of OT-I cells after one week (about ) is larger than the population of F5 T cells at day 18 (about ). Therefore, the population of F5 T cells never reaches its carrying capacity and the role of competition for resources can be neglected.

Figure 6

Density-dependent birth rate, , as a function of the population size, P. The parameter , with units of , represents the rate of growth under no competition and quantifies the level of reduction caused by the expansion of competing cells. Values for (shown in the inset) and are taken from[20], Table 1].

We estimate model parameters with the ABC-SMC algorithm[33]. Thus, the posterior distribution of the parameters is obtained by T sequential applications of the ABC algorithm, where the posterior obtained in each iteration is used as prior for the next one. This algorithm requires the definition of prior distributions for the first iteration, a distance function, a tolerance threshold for each iteration, and a perturbation kernel[33]. We assume all parameters are initially distributed according to a uniform prior distribution, as described in Table 1. When a prior distribution spans several orders of magnitude, the uniform distribution is taken over the exponent to efficiently explore parameter space. Given , the experimentally determined mean number of cells in generation g at time t, for , and its corresponding model prediction, for a particular choice of parameters , the distance function is defined aswhere is the set of time points and depends on the clonotype of interest, represents the standard deviation of the experimental data at time t and generation g, and G is the merged (and maximum) generation, . In practice, we define the first tolerance threshold, , in the ABC-SMC algorithm as the median value of the distances obtained from preliminary realisations, with the parameters sampled from the prior distributions in Table 1. The subsequent tolerance thresholds, , can be then defined as the median of the distance values obtained from the previous iterations of the algorithm. Finally, we use a uniform perturbation kernel to perturb the parameters during the sequence of iterations[33], and implement the algorithm for in the case of the multi-stage model and for the single stage one.

Table 1

Prior distributions for model parameters. Units for , and are inverse hours ().

Model parameters	Description	Prior distribution
C_0	Initial number of cells	C_0=10^x, x∼U(4,6)
N_0, N	Number of stages	U_discrete(1,50)
λ_0, λ	Birth rate	λ_0=10^y, λ=10^z, y,z∼U(-3,1)
α	Death rate slope	α=10^w, w∼U(-5,-1)

Before performing the Bayesian inference, we make some assumptions based on the experimental set up. Several studies have shown that the time to a first division is larger the time to subsequent divisions, since cells require time to become activated before they divide[13,15,16]. Thus, we assume that all generations but 0 are comprised of the same number of stages N, whereas generation 0 is characterised by stages. Similarly, cells in generation 0 proceed to divide with birth rate , whilst all the other generations have a birth rate . Therefore, in contrast to the inference in[29], the number of stages and N are free parameters in the model. On the other hand, we propose that the per cell death rate in a given generation is linear on the number of cell divisions that the cell has undergone[35,45]. We writewhere is a parameter to estimate. These linear death rates encode the fact that cells are more likely to die when they have already undergone several divisions[35,45]. Finally, the initial number of cells, , is considered a parameter to be estimated, since the actual number of transferred cells which make it to the lymph nodes or spleen cannot be measured.

The predictions obtained for each model, and for each clonotype (OT-I or F5), are shown in Fig. 7. We run the model with the parameters being sampled from the estimated posterior distributions and compute the median of all the simulations, which corresponds to the solid magenta (multi-stage model) and turquoise (exponential model) lines in Fig. 7. The bands around median predictions represent 95 confidence intervals. Data points are plotted with the standard deviation from the multiple experimental replicates. As shown in Fig. 7, the calibrated MS-G model successfully captures the dynamics of the proliferating T cell populations (OT-I and F5), whereas the single stage model significantly underestimates the expected number of cells beyond generation 1, particularly in the case of OT-I T cells. We made use of the corrected version of the Akaike Information Criterion ()[46,47] to quantify the better fit of the MS-G model accounting for the extra parameters and N in the multi-stage model. The values of for each model and clonotype are listed in Table 2. Despite the two extra parameters, the values of corresponding to the multi-stage model are significantly lower for both clonotypes. Overall, the MS-G model is able to explain the data from the OT-I transfer experiment better, since this data set is less noisy than the F5 set.

Figure 7

Exponential (solid turquoise line) and multi-stage (solid magenta line) model predictions compared to the data sets (orange dots) for OT-I (A) and F5 (B) T cells. Bars on data points represent their standard deviation. The expected number of cells in each generation is plotted as a function of time. These predictions represent the median value of simulations with the accepted parameter values from the posterior distributions. Shaded areas represent 95 confidence intervals.

Table 2

values for the exponential and multi-stage models calibrated with CFSE data of murine T lymphocytes.

Mathematical model	Cell type	Value of AIC_C
Multi-stage	OT-I T cells	50.4
Exponential	OT-I T cells	283
Multi-stage	F5 T cells	206
Exponential	F5 T cells	317

The marginal posterior distributions for each parameter are shown in green and blue in Figs. 8 and 9, for the multi-stage and exponential models, respectively, and the (uniform) prior distributions are plotted in red. Summary statistics of these posterior distributions are shown in Tables 3, 4, 5 and 6. Cell death is governed by the parameter , and is estimated to be low for both models and clonotypes, suggesting that cell death does not have a significant impact on the dynamics during lymphopenia, which is in fact dominated by cell division. This result is in agreement with Hogan et al.[20], where the death rate is assumed to be zero. The initial number of cells can be estimated with relative success, and does not seem to depend heavily on the model considered. On the other hand, cell division is governed by parameters , with in the exponential model. We note that in both models, and represent the mean time to the first and subsequent divisions, respectively. Although all division-related parameters can be estimated from the data, for both models and clonotypes, a correlation between the division rate and the number of stages is seen in the scatter plots of Fig. 10. Instead of plotting the marginal posterior distributions for these parameters, one can consider the posterior distribution for the mean times and (see Fig. 10). The fact that is never chosen as an accepted parameter value in the posterior distribution for the multi-stage model and the OT-I clonotype already suggests that a multi-stage representation of cell division is preferred for this clonotype. On the other hand for the F5 clonotype the marginal distribution for N shows a non-zero frequency for the value 1, but larger values of N are also represented in its posterior distribution. The mean time to both first and subsequent divisions, and , are significantly longer for the F5 clonotype than the OT-I. In fact, our results estimate that F5 T cells divide slowly compared to OT-I cells, requiring on average 192 h to carry out a first division (59 h taken by OT-I T cells), as shown in Fig. 10 for the multi-stage model. The time to subsequent divisions is represented by the blue histograms. Interestingly, our estimation of the mean time to first division of OT-I cells, on average 59 h, is close to the value obtained by Hogan et al.[20] (52 h when considering the best fit parameter estimates). In the case of F5 cells, we predict an average of 192 h to undergo their first division, whereas Hogan et al. obtained a value of 137 h. We note that the value 137 h is within the range covered by our predicted posterior distribution.

Figure 8

Posterior distributions (green and blue) for the parameters in the multi-stage (A) and exponential (B) model for OT-I T cells. In the exponential model, the number of stages for all generations is equal to 1, i.e., . Prior distributions are shown in red.

Figure 9

Posterior distributions (green and blue) for the parameters in the multi-stage (A) and exponential (B) model for F5 T cells. In the exponential model, the number of stages for all generations is equal to 1, i.e., . Prior distributions are shown in red.

Table 3

Summary statistics of OT-I clonotype posterior distributions for the multi-stage model.

Parameter	Minimum	Maximum	Mean	Median	Standard deviation
C_0	3.00×10^4	1.03×10^5	4.67×10^4	4.54×10^4	8.03×10^3
N_0	1	7	2.83	3	1.23
N	2	34	6.59	5	4.30
λ_0	1.98×10^-2	1.08×10^-1	4.64×10^-2	4.56×10^-2	1.45×10^-2
λ	2.80×10^-2	8.08×10^-1	1.48×10^-1	1.20×10^-1	1.01×10^-1
α	1.00×10^-5	5.97×10^-3	5.06×10^-4	1.76×10^-4	7.47×10^-4

Table 4

Summary statistics of F5 clonotype posterior distributions for the multi-stage model.

Parameter	Minimum	Maximum	Mean	Median	Standard deviation
C_0	4.74×10^4	1.85×10^5	9.26×10^4	8.94×10^4	2.22×10^4
N_0	1	10	3.01	3	1.53
N	1	35	2.42	2	2.57
λ_0	2.68×10^-3	7.20×10^-2	1.70×10^-2	1.47×10^-2	1.07×10^-2
λ	2.06×10^-3	5.88×10^-1	2.20×10^-2	9.54×10^-3	3.90×10^-2
α	1.00×10^-5	6.21×10^-3	1.35×10^-3	8.19×10^-4	1.40×10^-3

Table 5

Summary statistics for the posterior distributions of the exponential model for the OT-I clonotype.

Parameter	Minimum	Maximum	Mean	Median	Standard deviation
C_0	1.43×10^4	5.46×10^5	2.40×10^4	2.34×10^4	5.28×10^3
λ_0	1.56×10^-2	6.19×10^-2	2.40×10^-2	2.28×10^-2	5.69×10^-3
λ	7.45×10^-3	1.29×10^-2	9.88×10^-3	9.78×10^-3	1.10×10^-3
α	1.00×10^-5	7.23×10^-3	4.76×10^-4	1.32×10^-4	8.17×10^-4

Table 6

Summary statistics for the posterior distributions of the exponential model for the F5 clonotype.

Parameter	Minimum	Maximum	Mean	Median	Standard deviation
C_0	4.54×10^3	1.10×10^5	7.33×10^4	7.22×10^4	1.39×10^4
λ_0	2.50×10^-3	7.18×10^-3	4.24×10^-3	4.13×10^-3	9.28×10^-4
λ	1.79×10^-3	3.92×10^-3	2.77×10^-3	2.76×10^-3	4.66×10^-4
α	1.00×10^-5	1.52×10^-3	1.72×10^-4	7.30×10^-5	2.23×10^-4

Figure 10

Joint posterior distributions (left two plots) of the number of stages , N and the birth rates , . Marginal posterior distributions (right two plots) for the mean time to first and subsequent divisions estimated from the multi-stage model (third column) and the exponential model (fourth column). Panel A for OT-I T cells and B for F5 T cells.

Our results indicate that OT-I T lymphocytes require on average 59 h for their first division, and a bit less, 46 h, for subsequent divisions (see upper left plot of Fig. 10). Based on our Bayesian approach, we conclude that a multi-stage model with a constant division rate after the first division event, is a suitable description of lymphopenia-induced proliferation[5,29,40]. The MS-G model estimates that F5 cells take on average slightly less than 200 h to divide, both for the first or subsequent division rounds, as shown in the lower left plot of Fig. 10. This difference can be explained by the different characteristics of OT-I and F5 T cells, and was previously observed[20]. The posterior distributions of the expected time to subsequent divisions in the MS-G model, , and in the exponential one , shown as blue histograms in Fig. 10, indicate that the exponential model predicts a longer division time than the multi-stage model for both clonotypes. This can be justified by the implementation of the ABC-SMC algorithm. Indeed, when parameterising the exponential model, the algorithm tries to keep the distance between the model predictions and the experimental observations low. This leads to the choice of parameter sets which limit cell proliferation, as shorter division times in the exponential model would lead to an increase in cell numbers not observed in the data set, and thus larger distance values. This is why the estimated birth rates in the exponential model are lower than the ones in the multi-stage representation. As a result, the exponential model predicts a greater average division time than the multi-stage model for both clonotypes. Finally, our results indicate that for both clonotypes the exponential model (see Fig. 10) found a shorter time to first division than to subsequent ones, contradicting previous findings[13,15,16], which support longer first division times. This is related to the fact that, overall, the exponential model is not able to capture the observed cell dynamics for neither of the clonotypes, as can be seen in Fig. 7.

Discussion

We analyse a multi-stage model of cell proliferation and death, tracking cell generations, in a framework that retains the benefits of a Markov process. With particular choices of rates, the models are equivalent to others in the literature[13,27,29,40-42]. In the case study of “Case study: lymphopenia-induced proliferation” section, the MS-G model performs better than the exponential model of time to division. The model implemented here provides a flexible framework for estimating the birth and death rates that describe the dynamics of lymphocyte populations[48,49]. The representation retains the advantages of a Markovian approach, including analytical tractability in some cases, and computational efficiency of numerical simulations with the Gillespie algorithm[50,51]. The expected number of cells in each generation satisfy a set of linear differential equations. Further comparison of this and of published models[13,27,29,40-42] with different experimental datasets is the aim of future work.

It has been observed[13,15,16] that immune cells typically need longer to divide for the first time, whereas later divisions require shorter times[5]. It is possible to assume that divided and undivided cells have different probability densities of time to cell division in exponential and Smith-Martin models[7,35,45]. With the multi-stage model introduced here, the separation need not be explicit because it is incorporated in the generation-dependent parameters. A longer mean time to first division, , than mean time to subsequent divisions, , is a natural part of the framework. Extension of the mathematical analysis in “Analytical results” section to the case , and possibly generation-dependent death rate , , would be desirable.

Our calculations rely on the assumption that cells are independent of each other. In particular, no fate correlation is assumed between daughter cells and their progenitors, or between siblings. However, data sets from time-lapse microscopy of B and T cell families[10-12,14,16,17,25] show that division and death times for siblings are correlated, and “division destiny” is a familial characteristic[24]. A further potential extension of the MS-G model is the introduction of a population carrying capacity. In the model as described in “Analytical results” section, the mean number of cells over time either increases without bound, dies out or reaches a steady-state, depending on the relation between division-related parameters (birth rate and number of stages in the cell cycle), and the death rate. Competition for resources may be modelled using density-dependent birth and/or death rates[20,52], or by rates that depend on the time-dependent availability of resources[53].