## Abstract

To study the kinetics of lymphocytes, models have divided the cell population into subpopulations with different turnover rates. These have been called ‘kinetic heterogeneity models’ so as to distinguish them from ‘temporal heterogeneity models’, in which a cell population may have different turnover rates at different times, e.g. when resting versus when activated. We model labelling curves for temporally heterogeneous populations, and predict that they exhibit equal biphasic up- and downslopes. We show when cells divide only once upon activation, these slopes are dominated by the slowest exponent, yielding underestimates of the average turnover rate. When cells undergo more than one division, the labelling curves allow fitting of the two exponential slopes in the temporal heterogeneity model. The same data can also be described with a two-compartment kinetic heterogeneity model. In both instances, the average turnover rate is correctly estimated. Because both models assume a different cell biology but describe the data equally well, the parameters of either model have no simple biological interpretation, as each parameter could reflect a combination of parameters of another biological process. Thus, even if there are sufficient data to reliably estimate all exponentials, one can only accurately estimate an average turnover rate. We illustrate these issues by re-fitting labelling data from healthy and HIV-infected individuals.

## 1. Introduction

Despite great advances in immunological research during the last decades, estimating the kinetics of lymphocyte populations has proved quite difficult. There are widely divergent estimates of the production rates, division rates and lifespans of mouse and human lymphocyte populations [1]. As a consequence, it is poorly understood how memory is maintained, how the naive lymphocyte repertoire remains diverse and how homeostasis is brought about. Additionally, it remains poorly understood how viruses such as HIV, and diseases such as rheumatoid arthritis, disturb the normal population dynamics and deteriorate the functioning of the immune system.

Several *in vivo* labelling techniques have been developed that have enabled the generation of quantitative data on lymphocyte dynamics. Examples of labels are the fluorescent dye carboxy-fluorescein diacetate succinimidyl ester (CFSE), the base analogue 5-bromo-2′-deoxyuridine (BrdU) and the stable isotype deuterium. Major advantages of deuterium labelling are that it is unlikely to change the cellular dynamics, and that, unlike CFSE and BrdU, it can safely be used in humans, and hence provide insights into human lymphocyte turnover in health and disease. Animals and humans have been labelled with deuterated glucose (^{2}H_{2}-glucose) or heavy water (^{2}H_{2}O) for 1 or 2 days [2–6], 5 days [7] to one week [8] or several months [9–11], and have subsequently been followed to study the loss of label during a washout phase. After deuterium labelling, one sorts the cell population of interest, isolates the DNA from the cells, and uses mass spectrometry to determine the enrichment of deuterium in the DNA [12–15]. During the labelling period, some fraction of the hydrogen atoms in the body will be replaced by deuterium, and dividing cells will incorporate both deuterium and hydrogen in newly synthesized DNA molecules. Ultimately, the data correspond to a normalized fraction of labelled DNA, 0 ≤ *L*(*t*) < 1, which increases during the labelling phase, and which decreases after deuterium has been withdrawn. We will refer to these two phases of *L*(*t*) as the upslope and the downslope, respectively. Adding deuterium should not change the number of T cells, and thus one typically assumes steady state for the total cell population during the labelling experiment.

Although deuterium labelling is increasingly being used, the interpretation of the data in terms of the underlying kinetics of the cells has turned out to be notoriously difficult [1,8,10,16–18]. One major problem is that the estimated turnover rates differ markedly depending on whether ^{2}H_{2}-glucose or ^{2}H_{2}O has been used for the labelling. Another problem is that the upslopes and the downslopes depend on the period of labelling, with steeper slopes for shorter labelling periods [1] (figure 1*a*). For the downslopes, this is relatively well understood because Asquith *et al.* [16] demonstrated that whenever the population is kinetically heterogeneous, the labelled fraction will initially be enriched in cells with a higher than average turnover rate. For the upslopes, this is not understood. The initial upslope should reflect the average turnover rate [1,16], and there should be one unique labelling curve that monotonically approaches the asymptote, *L*(*∞*) = 1, where all cells are labelled. Thus, in situations where the estimated upslope seems to depend on the labelling period, one either has insufficient data to determine the one common slope, and/or one is fitting the true labelling curve with an incorrect model.

One of the major problems for interpreting labelling experiments is that we do not know which mathematical model should be employed to fit the data, and that truly realistic models have so many parameters that one would be overfitting the simple up- and downslopes of deuterium enrichment over time. The appropriate model may also depend on the cells types that have been labelled, as naive T cells, self-renewing memory T cells and clonally expanding effector T cells have very different dynamics. Even if one deals with a single population of sorted cells, the population may still be heterogeneous. Here, we will analyse two types of heterogeneity: heterogeneity in phenotype, i.e. what type of cell; and in history, i.e. what has recently happened to the cell [3], which have been called ‘kinetic heterogeneity’ and ‘temporal heterogeneity’, respectively [20]. Kinetic heterogeneity implies that there exist separate subpopulations of cells that have different turnover rates. Temporal heterogeneity is more subtle and implies that a cell may transiently have a different turnover rate. For example, lymphocytes after being activated may have an increased death rate, a phenomenon called activation-induced cell death.

Currently, there is reasonable consensus on what models to use for describing the labelling of a possibly heterogeneous population of cells that are turning over at different rates. Ganusov *et al.* [20] demonstrated that the several models that have previously been used for describing kinetically heterogeneous populations [8,16,21] can all be written as a model for *n* independent subpopulations, each turning over at a rate *d*_{i},
1.1where *α*_{i} is the fraction of cells with turnover rate *d*_{i}, and *t*_{end} defines the end of the labelling phase [20]. The first equation describes the kinetics of label acquisition from which one can derive the upslope. The second equation describes the kinetics of label loss after the end of the labelling period, from which one can derive the downslope. For the case of just one population, *n* = 1 and no asymptote, *α*_{1} = 1, this model is formally identical to the precursor product model that Hellerstein [22] used to describe deuterium data. For the case *n* = 1 and *α*_{1} < 1, the model has been shown [20] to be identical to the model of Asquith *et al.* [16], and to that used by Mohri *et al.* [8]. For *n* → *∞*, this model can be used to describe populations with a continuous distribution of turnover rates [20]. The average turnover rate of the model is defined as . The initial upslope reflects the average turnover rate, whereas the initial downslope is dominated by the subpopulation with the fastest turnover, and hence depends on the period of labelling. The difference between the up- and downslope diminishes when the length of the labelling periods increases [20].

Although equation (1.1) seems like an excellent general model for a kinetically heterogeneous population, it is not clear that it is the best model to describe the possibility of temporal dynamic heterogeneity in lymphocytes. Lymphocytes that have recently divided, and hence are enriched in deuterium, have gone through a process of activation and may transiently experience a faster death rate than a quiescent cell, by a process that is sometimes referred to as activation-induced cell death. This form of temporal heterogeneity is particularly important for lymphocytes that are undergoing a phase of clonal expansion [23], as is typical during an immune response, because clonal expansion is followed by a phase of contraction owing to rapid cell death by apoptosis [24]. Unfortunately, most models for clonal expansion have so many parameters that they would lead to overfitting labelling data [17,18,21,25], and therefore, we still lack a reliable model for cell labelling when temporal heterogeneity is involved.

Another example where temporal heterogeneity could be important is provided by the study of Choo *et al.* [26], who recently described the process by which CD8^{+} memory T cells maintain themselves in mice. Clones of memory cells were generated by infecting mice with lymphocytic choriomeningitis virus (LCMV), the cells were labelled with CFSE *ex vivo* and transferred to naive recipient mice. Drawing blood samples from the recipients at various later time points, it was shown that the clones maintained a stable size, and that cells completed single divisions stochastically with an average interval of 50 days [26]. Although it is not known whether the two daughter cells that result from a rare stochastic division of a cell recruited from an otherwise quiescent population have a transiently increased death rate (which would account for temporal heterogeneity), these data—for the first time—allow us to study what deuterium labelling in a self-renewing population of lymphocytes should look like with a mathematical model that we believe is correct because it captures the underlying biology. Here we present this model, study its properties and generalize it to other situations with temporal heterogeneity.

## 2. Model

Consider a population of resting cells that are recruited into division at an activation rate *a*. Cell division yields *c* daughter cells, where *c* = 2 corresponds to the single divisions that account for the renewal of an otherwise quiescent population of memory T cells [26], and *c* > 2 corresponds to a burst of divisions that characterize clonal expansion. Small non-overlapping bursts of clonal expansion may also occur during chronic infection [23], and may play a role in the normal maintenance of part of the CD4^{+} T cells expressing a memory phenotype [27]. It would obviously be best to model a burst of cell expansion by a cascade indexed by the number of divisions the cells have completed [23,28], but such a model typically has too many parameters to fit to limited experimental data. To keep the analyses tractable, we model clonal expansion as an instantaneous process, but we do not expect that our main results will change if clonal expansion were incorporated as a cascade. Moreover, because we assume steady state in our model, we will limit ourselves to small bursts (*c* = 32, corresponding to 5 divisions).

The activated daughter cells revert to a resting state at rate *r*, and resting and activated cells have death rates *d*_{R} and *d*_{A}, respectively. Whenever *d*_{A} > *d*_{R}, the model accounts for temporal heterogeneity. Without loss of generality, one can scale the total population size to one, so that *A* + *R* = 1, and write for the fractions of resting, *R*, and activated, *A*, cells
2.1

Since we expect that the system is at steady state and that deuterium labelling does not change the dynamics, we consider d*R*/d*t* = d*A*/d*t* = 0. Using *A* + *R* = 1, the steady state d*A*/d*t* = 0 gives that the fraction of dividing cells , and the steady state of d*R*/d*t* = 0 can be used to eliminate one parameter, e.g.
2.2

Substituting the latter into gives the fraction of activated cells , or 2.3

The average turnover of the population is defined as . This model is identical to that of Ribeiro *et al.* [18], who were modelling two subpopulations with different turnover rates. Equation (2.2) requires that (*c* − 1)*a* − *d*_{R} > 0, which implies that for *c* = 2 we require *d*_{R} < *a*. For biological reasons, we only consider *c* > 1 and *d*_{R} < *d*_{A}. Below we compare populations with different values of *c*, but with the same average turnover rate , which can be achieved by keeping (*c* − 1)*a* the same when *c* is changed, thus keeping *f* constant. Note that the average turnover rate, or the average death rate, , is not related to the average residence times in the two subpopulations, i.e. 1/(*a* + *d*_{R}) and 1/(*r* + *d*_{A}), respectively.

To model deuterium labelling, it is most convenient to write equations for the unlabelled fractions of resting and activated cells, *U*_{R} and *U*_{A}, respectively, during the labelling phase [18], i.e.
2.4where *U*_{R} + *L*_{R} = *R* and *U*_{A} + *L*_{A} = *A* and the total labelled fraction is defined as *L* = *L*_{R} + *L*_{A} = 1 − *U*_{R} − *U*_{A}. To obtain the downslope, it is preferable to write equations for the loss of the labelled fractions, and—conveniently—these are identical to those used to obtain the upslope [18] when *U* and *L* are interchanged, i.e.
2.5

### 2.1. Upslope

Because the model is linear, Ribeiro *et al.* [18] were able obtain the analytical solutions. For the labelling phase with the initial condition *U*_{A}(0) = *f* and *U*_{R}(0) = 1 − *f*, one obtains
2.6which, after some algebra, is identical to the upslope of the two-compartment variant of the kinetic heterogeneity model of equation (1.1). Here
2.7and
2.8and *r* is defined by equation (2.2). The upslope is a combination of *e*_{1} and *e*_{2}, where *e*_{1} < *e*_{2}. The initial upslope of equation (2.6), i.e. d*L*(*t*)/d*t* for *t* → 0, is the average turnover rate . Thus, although the structure of this solution is identical to that of equation (1.1) with *n* = 2, the parameters have completely different interpretations because the exponents *e*_{1} and *e*_{2} are complicated expressions that do not reflect the turnover rates, *d*_{R} and *d*_{A}, of the subpopulations *R* and *A*, and the parameter *α* is not the fraction of activated cells *f*.

### 2.2. Downslope

The general expression for the downslope is given by Ribeiro *et al.* [18], and is complicated because it involves the fraction of labelled resting cells, *L*_{R}(*t*_{end}), and the fraction of labelled activated cells, *L*_{A}(*t*_{end}). We therefore first consider the case where all cells are labelled at the end of labelling, i.e. *L*(*t*_{end}) = *L*_{R}(*t*_{end}) + *L*_{A}(*t*_{end}) = 1, because the solution can then be written as
2.9where the parameters are defined as in equation (2.7). The initial downslope reflects (because we consider the case where all cells are labelled at time *t*_{end}). For finite labelling periods *t*_{end}, it can numerically be shown that the general expression for the downslope [18] rapidly approaches the approximate solution
2.10which is indeed equal to equation (2.9) when all cells are labelled at *t*_{end}, i.e. when . Thus, the downslope of this temporal heterogeneity model strongly resembles the downslope of equation (1.1) for *n* = 2, albeit with different parameter interpretation, and for all practical purposes it will be difficult to distinguish between the models.

### 2.3. No death of quiescent cells

The death rate of resting cells, *d*_{R}, is the slowest timescale of this model, and for *c* = 2, we indeed require that *d*_{R} < *a*. Thus, it is enlightening to simplify the exponents by assuming *d*_{R} = 0. This yields
2.11where *c*′ = *c*/(*c* − 1). For *c* = 2, this simplifies further into
2.12which for low activation rates, *a*, approaches *e*_{1} ≃ *a*/2 and *e*_{2} ≃ *a*/2 + 2*d*_{A}. Using the Choo *et al.* [26] estimate of *a* = 0.02 d^{−1}, and giving recently divided cells a short lifespan, i.e. choosing *d*_{A} = 1 d^{−1}, the exponents would be *e*_{1} ≃ 0.01 d^{−1} and *e*_{2} ≃ 2.01 d^{−1}.

## 3. Results and discussion

Let us first study what one would expect for a deuterium labelling experiment of LCMV-specific CD8^{+} memory T cells, as described by Choo *et al.* [26]. The cells on average divide every 50 days and then do one division, i.e. we set *a* = 0.02 d^{−1} and *c* = 2. Considering long-lived resting cells, e.g. *d*_{R} = 0.001 d^{−1}, and short-lived activated cells, e.g. *d*_{A} = 1 d^{−1}, one calculates a reversion rate of *r* = 1.105 d^{−1}, and a steady-state fraction of activated cells of *f* = 0.01865, from equations (2.2) and (2.3), respectively. The average turnover rate of the cells is then *d*^{−1} (which is close to *fd*_{A} because we assume that resting cells are long-lived). We study deuterium labelling experiments in this system by numerically solving equations (2.4) and (2.5) for labelling periods of 10, 50 and 100 days (figure 2*a*). Although the analytical solutions of this model involve two exponentials, see equations (2.6) and (2.9), the up and downslopes shown in figure 2*a* are hardly biphasic. For the present case, the parameters of the solution of equation (2.6) are *α* = 0.995, *e*_{1} = 0.01 d^{−1} and *e*_{2} = 2.12 d^{−1}, from which one sees that the fast exponential describes only a minor fraction of the population (i.e. 1 − *α* ≃ 0.005). This explains why the system is hardly biphasic, and the slow exponential, *e*_{1}, largely determines the ‘observed’ upslope and downslopes. The downslope barely depends on the labelling period, and even after all cells have been labelled, the ‘observed’ downslope, *e*_{1} = 0.01 d^{−1}, is about half of the average turnover rate, as d^{−1}.

As shown by equations (2.7) and (2.8), the slowest exponential is determined by several parameters of the model, and is complicated to understand in biological terms. Because, in this example, there is a slow death rate of resting cells, one can obtain better intuition about the slowest exponential by letting *d*_{R} → 0 (see equation (2.12)). For the current parameters, the slowest exponential in equation (2.12) is *e*_{1} = 0.00995 d^{−1}, which is very close to the true *e*_{1} = 0.01 d^{−1}, and is close to its further approximation *e*_{1} ≃ *a*/2 for low values of *a* (as *a* = 0.02 d^{−1}). Thus, for the current parameters, the expected up- or downslopes from this model tend to look like a single exponential, reflecting half of the activation rate, which is here half of the true turnover rate. Summarizing, if one were to perform a deuterium labelling experiment on the LCMV-specific memory CD8^{+} T cells described by Choo *et al.* [26], and had somewhat imperfect experimental data that are fitted with equation (2.1), which we believe accurately summarizes the underlying biology, one would have no reason to conclude that the population is heterogeneous, and one would estimate an average turnover rate that is about half of the actual turnover rate.

If one were to fit these data not knowing that the underlying biology is described by the temporal heterogeneity model of equation (2.1), and instead used the general kinetic heterogeneity model assuming that memory T cells are composed of *n* subpopulations with different turnover rates, i.e. equation (1.1), one would obtain excellent fits to the data. Fitting the three labelling curves in figure 2*a* simultaneously to the one-compartment version of equation (1.1), i.e. *n* = 1 and *α*_{1} = 1, allowing just *t*_{end} to differ between the curves, describes the data with excellent quality, and one would estimate that d^{−1} (not shown). Note that for *n* = 1, equation (1.1) is equal to the simple precursor–product model [1,2,20]. Given the monophasic nature of the curves in figure 2*a*, the high quality of the fit with this simple model seems a natural result. As expected, the exponent that is picked up by this one compartment model is very similar to the slowest exponential, *e*_{1}, of the solution of the temporal heterogeneity model that underlies the data. As argued earlier, *e*_{1} is approximately half the activation rate, and half of the true turnover rate for these parameters. Hence, one would also make a twofold error when fitting these data perfectly with the wrong model.

The same data are obviously also well described by the two compartment variant of the kinetic heterogeneity model of equation (1.1) because adding two parameters to a model that already describes the data well can only make it better. Once sufficient data points are provided, the two compartment variant will fit the data in figure 2*a* significantly better than the *n* = 1 variant. Fitting the model to the 600 daily values of *L*(*t*) in figure 2*a*, we estimate that *α*_{1} = 0.995, *d*_{1} = 0.01 d^{−1}, *d*_{2} = 2.13 d^{−1} and an average turnover d^{−1}, which is very similar to the parameters of the true solution (see earlier text). This excellent fit was to be expected because the temporal heterogeneity model has the same solution for the upslope, equation (2.6), and a very similar solution for the downslope, equation (2.10), as those given in equation (2.1). The average turnover rate would therefore be estimated correctly, but the values of *α*_{1}, *d*_{1} and *d*_{2} would not describe the kinetics of any subpopulation in the data, which have turnover rates of *d*_{R} = 0.001 d^{−1} and *d*_{A} = 1 d^{−1}, respectively. One can sample the data in figure 2*a* every 3 days, and add 10 per cent noise by multiplying *L*(*t*) with a factor chosen from a normal distribution with mean *μ* = 1 and standard deviation *σ* = 0.1 to mimic a real labelling experiment. As expected, the information on the minor contributions (less than 0.5%) to the initial rapid up- and downslopes is lost from that sampled data, and fitting that data with equation (1.1) with *n* = 2 compartments gives excellent fits to the data, but no longer allows for estimating all three parameters correctly (not shown). The slowest exponent is the only parameter that is reliably estimated from the sampled data, i.e. with a 95% CI of 0.0102 ≤ *d*_{1} ≤ 0.0108 d^{−1}, where the ranges denote 95% CI based on 500 simulations bootstrapping the residuals [29]. For *α*, the 95% CI, 0.98 ≤ *α*_{1} ≤ 1, includes one, which would correspond to the *n* = 1 model. The fast exponent *d*_{2}, and the average turnover rate , cannot be estimated as they have large confidence limits, i.e. 0.074 ≤ *d*_{2} ≤ 85.7 d^{−1}, and d^{−1}.

### 3.1. Clonal expansion

Using CFSE labelling, Choo *et al.* [26] demonstrated that both CD8^{+} LCMV-specific memory T cells and CD8^{+} memory T cells in non-immunized mice typically perform only one division per renewal event, and that these events are on average 50 days apart. Conversely, using BrdU labelling, Ki67 staining and TCR sequencing, Younes *et al.* [27] show that both in conventional non-immunizied mice and in germ-free mice, CD4^{+} T cells expressing a memory phenotype have an almost fivefold faster turnover, and that division sometimes involves clonal bursts. We now expect that labelling cells performing single divisions with deuterium would give labelling curves that can be described well with a single exponential. Deuterium labelling curves from the whole CD45RO^{+} effector/memory population of CD4^{+} or CD8^{+} memory T cells in humans [10] and mice [30] tend to be biphasic, and hence reflect some form of heterogeneity (figure 1*b*). This heterogeneity could be temporal, kinetic or both, and it would seem very reasonable that effector/memory cells on average perform more than one division when they are activated, because some of them will be specific for persistent infections and environmental antigens.

We investigate the impact of clonal expansion by setting *c* > 2 without changing the average turnover rate and the fraction of activated cells. Choosing *c* = 32 we scale *a*, i.e. choose *a* = 0.02/31 = 0.0006452 d^{−1}, and compute from equation (2.2) that *r* = 0.08659 d^{−1} (*f*, *d*_{R}, *d*_{A} and remain the same). Simulating the same three labelling experiments in a system having exactly the same average turnover rate now gives strongly biphasic curves with a markedly lower accrual of label (figure 2*b*). The slower enrichment can intuitively be understood from the fact that most labelled cells die before they return to the resting state. Substituting the parameter values into equation (2.7), the true solution is now described by *α* = 0.983, *e*_{1} = 0.0016 d^{−1} and *e*_{2} = 1.09 d^{−1}, implying that the contribution of the fast exponential, (1 − *α*), has increased from 0.5 to about 2 per cent, and that both the slow and the fast exponents decreased. Note that the approximation of equation (2.12) no longer holds because now *a* < *d*_{R}, and *d*_{R} is no longer the slowest time scale. The labelling period has hardly any effect on the downslopes, which run reasonably parallel, but has a large effect on the initial contribution of the fast exponent to the downslope (see the magnitude of the initial fall in labelling in figure 2*b*). After a short labelling period, the label is not yet equally divided over the resting and activated compartments, and the activated compartment will contribute more to the total loss of label.

The initial upslopes in figure 2*a*,*b* are identical, i.e. for *t* → 0 the slope , but at later times the label accrual occurs slower when there is clonal expansion (figure 2*b*). Estimating turnover rates by fitting the whole labelling curve with a one-compartment model, or from the peak value at the end of the labelling phase, would thus be incorrect. Fitting these data with models other than the true equations (2.6) and (2.9) clearly requires a model with more than one exponential, and because equation (1.1) with *n* = 2 has an almost identical solution it almost perfectly describes the data with the same three parameters (i.e. *α*_{1} = 0.983, *d*_{1} = 0.0016 d^{−1} and *d*_{2} = 1.08 d^{−1}). Hence, having biphasic data, a kinetic heterogeneity model would get the average turnover rate correct, but the individual parameter estimates would not be related to the size and turnover of its subpopulations, as one might have expected.

### 3.2. Activation

Because the model was parameterized based on the Choo *et al.* [26] data on memory cells undergoing infrequent divisions, the fraction of activated cells in the model is small, i.e. *f* ≃ 0.018. To study labelling in systems with more activation, we increase the activation rate by 25-fold to *a* = 0.5 d^{−1} and again compare single divisions (i.e. *c* = 2) with clonal expansion (i.e. *c* = 32). Now the fraction of activated cells is *f* = 0.3329 and the average turnover is (figure 3). For *c* = 2, the parameters of the solutions are *α* = 0.94, *e*_{1} = 0.22 d^{−1} and *e*_{2} = 2.28 d^{−1}; for *c* = 32, they are *α* = 0.69, *e*_{1} = 0.017 d^{−1} and *e*_{2} = 1.03 d^{−1}. Owing to the fast activation, labelling accrual occurs rapidly. Even when a large fraction of the cells has recently divided, and transiently have a markedly increased death rate, the labelling curves show little evidence for heterogeneity when cells undertake one division at a time (figure 3*a*). The three labelling curves for *c* = 2 are well described with the *n* = 1 variant of equation (1.1) (not shown), but with an average turnover rate of *d*_{1} = 0.23 d^{−1}, which agrees with the *e*_{1} ≃ *a*/2 of equation (2.12) as *a* = 0.5, but which is only about 70 per cent of the true average turnover rate, .

The curves become biphasic when we allow for clonal expansion by setting *c* = 32 (figure 3*b*). The period of labelling again hardly affects the downslopes (which run approximately parallel for the three labelling periods), and largely affects the contribution of the fast exponent. Given the similar solutions, fitting the data with the two compartment model of equation (1.1) gives an excellent fit with very similar parameters i.e. *α* = 0.69, *d*_{1} = 0.017 d^{−1} and *d*_{2} = 1.03 d^{−1} (not shown). This excellent fit would again provide no information on the turnover of subpopulations, and would only get the average turnover rate right.

### 3.3. Memory T cell turnover in man

Healthy volunteers and HIV-infected patients have been labelled with deuterated glucose [2,8] and with deuterated water [9,10,19]. The resulting T cell labelling curves tend to be somewhat biphasic. Ribeiro *et al.* [17] interpreted the labelling data from a 7 day labelling experiment with deuterated glucose [8] using a variation of equation (2.1) as a model for kinetic heterogeneity with a slow and fast subpopulation of the total T cell pool. For instance, the resting and activated subpopulations could correspond to naive and memory T cells, respectively. The main difference is that instead of using a clonal expansion parameter, *c*, they explicitly allowed activated cells to proliferate by adding a *pA* term to equation (2.1). Since neither *c* nor *p* appears in equations (2.4) and (2.5) for the labelled fractions, one should be able to fit these data in exactly the same manner with the temporal heterogeneity model presented here (and, as a matter of fact, with the kinetic heterogeneity model of equation (1.1)).

Ribeiro *et al.* [17] fitted that data by assuming that on the timescale of the experiment, resting cells hardly die; i.e. they set *d*_{R} = 0. This seems a reasonable assumption because the peak enrichment in healthy volunteers after about nine weeks of labelling with heavy water is approximately 20 per cent, and fitting the up- and downslope with the kinetic heterogeneity model of Asquith *et al.* [16], the average lifespan of naive and memory T cells in volunteers was estimated to be several years and several months, respectively [10]. Adopting the *d*_{R} = 0 assumption, we computed the corresponding parameters of the temporal heterogeneity model from the parameter estimates obtained by fitting the two compartment model of Ribeiro *et al.* [17] to the deuterated glucose data of Mohri *et al.* [8] (table 1 and figure 4).

The labelling curves in HIV-1-infected patients in Mohri *et al.* [8], Ribeiro *et al.* [17], and the two examples in figure 4, are much more biphasic than those in healthy volunteers. Comparing the means over all individuals, the average turnover rate of CD4^{+} and CD8^{+} T cells increases 5.7-fold and 9.1-fold, respectively, in infected individuals compared with uninfected controls (table 1). Ribeiro *et al.* [17] found that in CD4^{+} T cells HIV infection increases the activation rate *a* (threefold), and the proliferation/death rate, *d*_{A} (2.4-fold). In CD8^{+} T cells, these parameters were not changed significantly by HIV infection (despite the threefold increase in the activation rate *a*), and the main difference was the size of activated compartment, *f* (see the asterisks in table 1). The increased death rate of CD4^{+} T cells could be attributed to the killing of CD4^{+} T cells by HIV-1, which would help us to explain that the fraction of activated CD4^{+} T cells, *f*, was not significantly increased (the twofold increase in *f* was borderline significant, and much less than the 4.3-fold increase of *f* in CD8^{+} T cells) [17].

Recomputing these parameter estimates into those of the temporal heterogeneity model, we find very similar results. The major new finding is the large increase in the clonal expansion of CD8^{+} T cells with HIV infection (table 1), i.e. *c* increases 14.2-fold, which seems a very natural result, given the specific immune responses to the virus. The absence of an increase in *c* in CD4^{+} T cells could be due to the fact that CD4^{+} T cells expand less during clonal expansion [31], and that clonally expanding CD4^{+} T cells are preferentially infected during HIV-1 infection [32,33]. Note the slow reversion rates of recently divided cells, *r*, with corresponding long expected residence times in the activated state of 1/*r* ≃ 50 days. This could be due to the fact that the clonal expansion phase in the temporal heterogeneity model is instantaneous, and that *r* would become faster if we were modelling clonal expansion with a cascade of equations [23,28]. Fitting a cascade model to labelling data is not straightforward because (i) one typically needs additional parameters, and (ii) one needs to estimate the optimal integer number of divisions that cells complete during clonal expansion.

Summarizing, reinterpreting the two compartment model in terms of temporal heterogeneity suggests that differences in activation and clonal expansion can equally well explain the observed differences in the labelling of volunteers and HIV-1-infected patients. The interpretation of the CD4^{+} T cell data remains similar. Since we now find evidence for a significant increase in clonal expansion of CD8^{+} T cells, it remains difficult to have a mechanistic understanding of the several-fold increase in the average turnover rate of both CD4^{+} and CD8^{+} T cells, and how an increase in the average turnover rate determines the depletion of CD4^{+} T cells, and not CD8^{+} T cells, in HIV-1-infected patients.

## 4. Conclusion

We have shown that a realistic model for temporal heterogeneity has solutions that are practically indistinguishable from those of a general model for kinetic heterogeneity. Because we typically do not know whether heterogeneity is temporal or kinetic, one should be careful with interpreting the compartment sizes, and the turnover rates estimated with kinetic heterogeneity models. The most reliable estimate that one can obtain from deuterium labelling is the average turnover rate. Interestingly, we have reached the same conclusion about the modelling of BrdU data [21,25]. In reality, there will be both kinetic and temporal heterogeneity, and a model can easily be written splitting each population of the kinetic heterogeneity model into the resting, *R*, and activated *A* fractions of equation (2.1). One would have to expand equation (2.1) with an index *i*, and allow each subpopulation, *R*_{i} + *A*_{i}, to comprise a certain fraction *α*_{i} of the total population, i.e. ∑*α*_{i}(*R*_{i} + *A*_{i}) = 1. A model with two kinetically different subpopulations would have four exponentials, which are too many to identify from the typical labelling curves that are currently available.

A second problem that we identified is that the average turnover rate will be underestimated if there is heterogeneity, especially if the data contain too little evidence of the fast exponentials because of their small contribution. In such a case, one would typically fit equation (1.1) to the data, and use the *F*-test to decide how many compartments, *n*, are required for an optimal fit to the data. Our results suggest that one should also test how the estimated average turnover rate, , depends on the number of compartments. If increasing *n* allows to converge to a particular value, this could be the best estimate for the average turnover rate, even though the *F*-test would argue that one is overfitting the data. One should therefore estimate and its confidence limits for several values of *n*, to see how the average turnover rate depends on the number of compartments.

The issues discussed here should also apply to other labelling techniques, such as CFSE and BrdU, although each of these have specificities that make a full discussion beyond the scope of the current study. In particular, the underestimation of the average cell cycle time by modelling of CFSE data has been discussed before: different death rates during the different phases of the cell cycle can give quite different estimates of the cell cycle time when compared with models assuming a constant death rate [34,35]. Likewise, data obtained with BrdU labelling have been interpreted with different models, and because the values of their parameters strongly depended on the choice of the model, it was proposed to just estimate an average turnover rate from BrdU data [21,25]. Thus, none of the current labelling methods provides sufficient information to accurately describe the kinetics of the possible subpopulations of the population of cells studied in a labelling experiment. To understand the kinetics of human T cells, where deuterium labelling remains the main source of data, we have shown that the estimation of the kinetics of memory T cell populations remains difficult because of its expected kinetic and temporal heterogeneity: memory T cell data obtained with deuterium labelling tend to provide too little information for estimating anything more than just the average turnover rate of this quite heterogeneous population.

## Acknowledgments

Portions of this work were done under the auspices of the U. S. Department of Energy under contract DE-AC52-06NA25396 and supported by NIH grants AI028433, P20-RR018754 and National Center for Research Resources and the Office of Research Infrastructure Programs (ORIP) through grant 8R01-OD011095-21 (ASP). R.J.D.B. thanks the Santa Fe Institute for their support and hospitality during a visit when this work was initiated.

- Received February 26, 2012.
- Accepted March 27, 2012.

- This journal is © 2012 The Royal Society