A mathematical model of the effects of aging on naive T-cell population and diversity

TL;DR

This paper presents a mechanistic mathematical model to analyze how aging affects naive T-cell population and diversity, highlighting thymic output decline as a key factor in immunosenescence.

Contribution

The paper introduces a novel quantitative model of T-cell dynamics that estimates diversity changes with age, addressing experimental challenges in full organismal analysis.

Findings

T-cell diversity decreases with age due to reduced thymic output.

Loss of small immunoclones significantly impacts overall T-cell diversity.

The model provides quantitative estimates of T-cell population dynamics.

Abstract

The human adaptive immune response is known to weaken in advanced age, resulting in increased severity of pathogen-born illness, poor vaccine efficacy, and a higher prevalence of cancer in the elderly. Age-related erosion of the T-cell compartment has been implicated as a likely cause, but the underlying mechanisms driving this immunosenescence have not been quantitatively modeled and systematically analyzed. T-cell receptor diversity, or the extent of pathogen-derived antigen responsiveness of the T-cell pool, is known to diminish with age, but inherent experimental difficulties preclude accurate analysis on the full organismal level. In this paper, we formulate a mechanistic mathematical model of T-cell population dynamics on the immunoclonal subpopulation level, which provides quantitative estimates of diversity. We define different estimates for diversity that depend on the individual number of cells in a specific immunoclone. We show that diversity decreases with age primarily due to diminished thymic output of new T-cells and the resulting overall loss of small immunoclones.

Equations

$$\frac{\mathrm{d}N(t)}{\mathrm{d}t}=\gamma (t)+pN(t)-\mu (N)N(t),$$

$$\mu (N)={\mu }_0+\frac{{\mu }_1{N}^2}{N^2+K^2},$$

$$\Delta (\bar{N})=\frac{({\bar{N}}_{\mathrm{o}}-{\bar{N}}_{\mathrm{y}})}{{\bar{N}}_{\mathrm{y}}}.$$

$$\frac{\mathrm{d}{n}_i}{\mathrm{d}t}=\frac{\gamma (t)}{\Omega }+p{n}_i-\mu (N)n_i,$$

$${\hat{c}}_k(t)=\sum _{i=1}^{\Omega }{\delta }_{n_i(t),k},$$

$$\frac{\mathrm{d}{c}_k(t)}{\mathrm{d}t}$$

$$R_q(t)=\sum _{k\ge q}{c}_k(t),$$

$$\Delta ({\bar{R}}_q)\equiv \frac{({\bar{R}}_{\mathrm{o}}(q)-{\bar{R}}_{\mathrm{y}}(q))}{{\bar{R}}_{\mathrm{y}}(q)}.$$

$$\mathbb{E}[c_k^Y]=\sum _{j=1}^{R}jP\left(c_k^Y=j\right),$$

$$\mathbb{E}[c_k^Y]$$

$$\frac{\mathrm{d}{c}_M(t)}{\mathrm{d}t}$$

$$c_1^{\mathrm{ss}}$$

$$c_k^{\mathrm{ss}}$$

$$c(N;{\gamma }_0)=\left(p_0-({\mu }_0+{\mu }_1)\right)N^3+{\gamma }_0{N}^2+(p_0-{\mu }_0)K^{2}N+{\gamma }_0{K}^2.$$

$$\underset{{\gamma }_0\to 0}{\lim }\underset{M\to \infty }{\lim }\sum _{k=1}^{M}k{c}_k^{\mathrm{ss}}>0$$

$$L({\gamma }_0)\le \underset{M\to \infty }{\lim }\sum _{k=1}^{M}k{c}_k^{\mathrm{ss}}\le U({\gamma }_0),$$

$$\frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}$$

$$=\frac{2{\mu }_1{K}^2{N}_{\mathrm{ss}}}{(N_{\mathrm{ss}}^2+K^2{)}^2}\left[\frac{-(N_{\mathrm{ss}}^2+K^2)}{3(p_0-({\mu }_0+{\mu }_1))N_{\mathrm{ss}}^2+2{\gamma }_0{N}_{\mathrm{ss}}+(p_0-{\mu }_0)K^2}\right]$$

$$=\frac{-2{\mu }_1{K}^2{N}_{\mathrm{ss}}}{(N_{\mathrm{ss}}^2+K^2)\left[3(p_0-({\mu }_0+{\mu }_1))N_{\mathrm{ss}}^2+2{\gamma }_0{N}_{\mathrm{ss}}+(p_0-{\mu }_0)K^2\right]}$$

$$\underset{{\gamma }_0⟶0}{\lim }\frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}$$

$$=\frac{-2{\mu }_1{K}^2{N}_{SS}(0)}{(N_{SS}(0{)}^2+K^2)\left[-2(p_0-{\mu }_0)K^2\right]}$$

$$>0$$

$$f(p/\mu (N_{\mathrm{ss}}({\gamma }_0));{\gamma }_0)$$

$$={\left(1-\frac{p}{\mu (N_{\mathrm{ss}}({\gamma }_0))}\right)}^{\frac{-{\gamma }_0}{p\Omega }}\cdot \frac{{\gamma }_0}{p\Omega }{\left(1-\frac{p}{\mu (N_{\mathrm{ss}}({\gamma }_0))}\right)}^{-1}$$

$$=A({\gamma }_0)\cdot B({\gamma }_0)$$

$$\ln (A_0)$$

$$=\frac{-1}{p\Omega }\underset{{\gamma }_0⟶0}{\lim }\frac{\ln \left(1-\frac{p}{\mu (N_{\mathrm{ss}}({\gamma }_0))}\right)}{{\gamma }_0^{-1}}$$

$$=\frac{-1}{p\Omega }\underset{{\gamma }_0⟶0}{\lim }\frac{{\left(1-\frac{p}{\mu (N_{\mathrm{ss}}({\gamma }_0))}\right)}^{-1}\frac{\mathrm{d}}{\mathrm{d}{\gamma }_0}\left(-\frac{p}{\mu (N_{\mathrm{ss}}({\gamma }_0))}\right)}{-{\gamma }_0^{-2}}$$

$$=\frac{1}{p\Omega }\underset{{\gamma }_0⟶0}{\lim }{\gamma }_0^2{\left[1-\frac{p}{\mu (N_{\mathrm{ss}}({\gamma }_0))}\right]}^{-1}\left[p\mu (N_{\mathrm{ss}}({\gamma }_0){)}^{-2}\frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}\right]$$

$$=\frac{1}{p\Omega }\underset{{\gamma }_0⟶0}{\lim }\left[\frac{{\gamma }_0^{2}p\frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}}{\mu (N_{\mathrm{ss}}({\gamma }_0){)}^2-p\mu (N_{\mathrm{ss}}({\gamma }_0))}\right]$$

$$=\frac{1}{\Omega }\underset{{\gamma }_0⟶0}{\lim }\left[\frac{2{\gamma }_0\frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}+{\gamma }_0^2\frac{{\mathrm{d}}^2\mu }{\mathrm{d}{\gamma }_0^2}}{(2\mu -p)\frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}}\right]$$

$$=\frac{1}{\Omega }\left[\frac{2{\gamma }_0{\lim }_{{\gamma }_0\to 0}\frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}+{\gamma }_0^2{\lim }_{{\gamma }_0\to 0}\frac{{\mathrm{d}}^2\mu }{\mathrm{d}{\gamma }_0^2}}{p{\lim }_{{\gamma }_0\to 0}\frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}}\right],$$

$$B_0$$

$$=\underset{{\gamma }_0\to 0}{\lim }\frac{({\gamma }_0/p\Omega )}{\left(1-\frac{p}{\mu (N_{\mathrm{ss}}({\gamma }_0))}\right)}$$

$$=\underset{{\gamma }_0\to 0}{\lim }\frac{(1/p\Omega )}{p\mu (N_{\mathrm{ss}}({\gamma }_0){)}^{-2}\frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}}$$

$$=\underset{{\gamma }_0\to 0}{\lim }\frac{\mu (N_{\mathrm{ss}}({\gamma }_0){)}^2}{p^2\Omega \frac{\mathrm{d}\mu }{\mathrm{d}{\gamma }_0}}$$

$$>0.$$

$$\frac{{\gamma }_0}{\Omega }\sum _{i=1}^M\frac{1}{i!\mu (N_{\mathrm{ss}}({\gamma }_0){)}^{i-1}}\left(\prod _{j=1}^{i-1}\left[\frac{{\gamma }_0}{\Omega }+jp\right]\right)$$

Full text

††journal: Bulletin of Mathematical Biology

1 Introduction

Immunosenescence underlies poor health outcomes in the aging population, including diminished vaccine efficacy (Poland et al. 2010; McElhaney and Dutz 2008; Fleming and Elliot 2008), increased susceptibility to disease (including irregular presentation, intensified symptoms, longer recovery times, and increased mortality) (Thomas-Crussels et al. 2012), as well as a heightened risk of cancer (Ginaldi et al. 2001). This degradative aging process of the human immune system originates from extensive fundamental changes to the size and functionality of immune cell pools, and the structure of lymphatic tissues in which they develop and operate (Salam et al. 2013).

Among the many changes associated with immunosenescence (Globerson and Effros 2000), the T-cell compartment is arguably the most damaged (Wick et al. 2000; Gruver et al. 2007). The T-cell pool is comprised of subpopulations of antigen-inexperienced naive cells, and antigen-experienced memory cells, the latter of which retain immunological record of previous infections. The human immune compartment maintains $\sim 10^{12}$ T-cells in total, of which $\sim 10^{11}$ are naive (Jenkins et al. 2009; Trepel 1974). During aging, the population of naive T-cells declines in overall size, while the population of memory T-cells undergoes extensive proliferation, thereby reversing the balance of naive and memory T-cells that had persisted at younger ages (Globerson and Effros 2000; Fagnoni et al. 2000). The expansion of memory T-cells further enhances immunological memory of previously-encountered antigens, reinforcing existent immune protection. The remaining naive pool experiences loss of T-cell receptor (TCR) “structural diversity” (Goronzy et al. 2007,  2015b)–the number of distinct TCR complexes present across the entire naive pool. The diversity of T-cell clones, or “immunoclones”, characterized by the number of distinct TCR complexes among the cell population, provides the extent of antigen specificity. Unique TCR complexes are generated during T-cell development in the thymus, via recombination of genes encoding the V and J domains of the TCR$\alpha$ chain and the V, J, and D domains of the TCR$\beta$ chain, along with additional insertion and deletion of nucleotide fragments (Murphy 2012). Combinatorially, a possible $\Omega_{0}\sim 10^{15}-10^{20}$ unique TCR complexes may be assembled via this rearrangement process (Laydon et al. 2015), but only $\Omega\sim(0.05)\times\Omega_{0}$ of those rearrangements are functionally viable (Yates 2014), as determined by positive and negative selection tests in the thymus, which screen for appropriate reactivity to self-peptide/MHC molecules. Each TCR is activated by at least one peptide fragment presented via MHC molecules on the surface of an antigen-presenting cell, thus loss of naive TCR structural diversity limits the number of new antigens to which the full naive T-cell pool can respond. Naive cells are also suspected to suffer major functional deficiencies in aging, such as diminished binding affinity and proliferative capacity after antigenic stimulation (Moro-García et al. 2013). While these effects have been studied mostly using murine models to date (Appay and Sauce 2014), they are not yet well understood in humans and are beyond the scope of this paper.

The total abundance of naive T-cells, which inhabit both blood and lymphatic tissue, can be reliably estimated from measurements in small samples (Westermann and Pabst 1990; Bains et al. 2009a). Recently, Westera et al. (2015) estimated an $\sim 52\%$ decrease in the naive T-cell population in aging. In contrast, accurate estimation of full-organism TCR structural diversity is currently impeded by experimental imprecision and the inability to extrapolate small sample data to the full organism (Laydon et al. 2015). Experimentation typically entails DNA sequencing of the TCR$\alpha$ or–more commonly–$\beta$ chain, in particular the complimentarity-determining region 3 (CDR3), which is the site of TCR binding to antigenic peptide and most significant basis for diversity (Murphy 2012).

Increasingly sophisticated deep sequencing methods have improved estimates for the lower bound on TCR diversity but direct estimation of TCR diversity remains a challenge due to various experimental complications, such as the inability to detect rare clonotypes, sequencing errors, and inaccurate measurement of clonotype frequencies resulting from inconsistencies in polymerase chain reaction (PCR) amplification (Laydon et al. 2015). Predicting full-organism TCR diversity from a small sample is typically formulated as an “unseen species problem”, and one of many canonical solutions to such a problem is employed in conjunction with experimental data (Chao 1984; Chao and Lee 1992; Colwell and Coddington 1994), but the true relationship between sample and full diversity is fundamentally elusive.

Despite variations across experimental measurements of TCR diversity, its age-related loss has been consistently observed. An early study conducted by Naylor et al. (2005) predicted a TCR$\beta$ chain diversity of $\sim 2\times 10^{7}$ that persisted in donors through age $60$, before dropping by two orders of magnitude to $\sim 2\times 10^{5}$ at age $70$. More recently, Britanova et al. (2014) collected samples from donors of all ages and observed an approximately linear decrease in TCR$\beta$ CDR3 diversity from $\sim 7\times 10^{6}$ in youth ($6-25$ years) to $\sim 2.4\times 10^{6}$ in advanced age ($61-66$ years). Qi et al. (2014) obtained a particularly high lower bound estimate of $\sim 10^{8}$ unique TCR$\beta$ sequences in youth ($20-35$ years), which declined two- to five-fold in advanced age ($70-85$ years).

Note that only the TCR$\beta$ chain is sequenced in these experiments. Sequencing of both the $\alpha$ and $\beta$ chains would potentially produce a more accurate measure of TCR diversity, but the same experimental limitations preclude complete analysis. The measurement of diversity is further complicated by the potentially large disparity between structural diversity and “functional diversity”–that is, the number of antigens to which the T-cell pool is capable of responding. Due to the potential for crossreactivity, in which one TCR might respond to many structurally similar peptide fragments, it is possible that actual TCR diversity is much higher than structural diversity indicates. It has been speculated that one TCR might respond to as many as $10^{6}$ different peptide epitopes (Mason 1998).

To obtain lifetime estimates of TCR structural diversity, and develop an informed context for discussion of functional diversity, we introduce a mechanistic mathematical model of the generation and replenishment of the lymphocyte pool from birth through the end of life. Although experimental assessments of full-system information remain challenging, measurements for the dynamics of each component related to the T-cell population can be found throughout the literature. Our mathematical approach combines the knowledge of these individual components to study their interplay, leading to an understanding of the full-system dynamics. By extending previous model studies of total cell counts (Mehr et al. 1996,  1997; Ribeiro and Perelson 2007; Bains et al. 2009a,  b; Hapuarachchi et al. 2013; Murray et al. 2003; Reynolds et al. 2013), our multi-component formulation is able to efficiently track the total number of distinct T-cell clones, allowing for a full-system assessment of TCR structural diversity.

2 Mathematical Models and Results

We develop our mathematical model by first constructing the equation governing the total population size of the naive T-cell pool in Sec. 2.1, through which we quantitatively constrain the primary parameters of our model using experimental measurements found in previous literature. The model that describes the evolution of immunoclones is derived in Sec. 2.2, allowing us to define and estimate the diversity of the T-cell population in Sec. 2.3. In Sec. 2.4, we inspect the impact of sampling on the estimate of immunoclone diversity, as in practice it is only possible to extract a small fraction of the entire T-cell population from a body.

2.1 Total T-cell population model

There are three fundamental immunological mechanisms that sustain the naive T-cell pool: 1) export of mature naive T-cells from the thymus, 2) peripheral proliferation, and 3) cell removal from the naive pool due to death or phenotypic changes. These basic mechanisms constitute a birth-death-immigration process described by the ordinary differential equation,

[TABLE]

where $N(t)$ denotes the total T-cell count, $\gamma>0$ denotes the rate of thymic output, $p>0$ denotes the rate of proliferation, and $\mu(N)>0$ denotes the rate of population-dependent regulated cellular death or loss of naive phenotype.

While more complex feedback mechanisms have been proposed  (Mehr et al. 1997), other experiments have shown that thymic export is independent of naive T-cell counts (Ribeiro and Perelson 2007; Berzins et al. 1998; Metcalf 1963), it is well-established that the export rate consistently decays throughout the human lifespan (Murray et al. 2003). The lifelong decline of thymic export is caused by thymic involution which leads to degradation of structural integrity and functional capacity of the thymus with age (Steinmann et al. 1985). The age dependence of the thymic export rate of newly-trained T-cells is often approximated by an exponentially decaying function, $\gamma(t)=\gamma_{0}e^{-at}$, where $\gamma_{0}>0$ is the maximum rate of thymic output that arises in early years, and $a>0$ is the rate of decrease in thymic output.

The immune systems of vertebrates maintain a healthy amount of naive T-cells through complex homeostatic mechanisms, which include controlled production and distribution of common gamma chain cytokines, particularly IL-7, to the naive pool (Fry and Mackall 2005). IL-7 is secreted by stromal and endothelial cells in the thymus, bone marrow, and lymphatic tissue, providing T-cells with necessary survival signals. In lymphoreplete conditions, competition for this limited resource regulates population size (Bradley et al. 2005; Tan et al. 2001; Vivien et al. 2001), but in lymphopenic conditions, high levels of IL-7 resulting from low T-cell counts can even stimulate cellular proliferation. While IL-7 concentration may be explicitly formulated in a mathematical model of the peripheral T-cell population, as in the work of Reynolds et al. (2013), most models incorporate IL-7 regulation implicitly in the form of carrying capacity, assuming quick equilibration in a state of competition for IL-7 in the presence of a given number of T-cells. Such simplification commonly leads to the dependence on total cell counts of both cell proliferation and cell death rates, considering the cytokine’s dual role under lymphoreplete and lymphopenic conditions described above. Our model assumes cell-count dependence only of the cell death rate, focusing on scenarios of healthy aging, i.e., lymphoreplete conditions. We thus assume an $N$-dependent cell death rate of the form

[TABLE]

where the first term, $\mu_{0}>0$, is the basal rate of cellular death. The second one describes the IL-7-mediated regulation of cell death, with $\mu_{1}>0$ representing the maximal increase to the death rate as $N\to\infty$. The quantity $K$ is analogous to a “carrying capacity” and dictates the population at which signalling induced death starts to limit the population. The constant rate of cellular proliferation under healthy conditions is supported by recent studies of Westera et al. (2015), showing nearly identical naive proliferation rates at young and old ages during moderate age-related non-lymphopenic loss of naive cells. IL-7 induced proliferation can arise in unhealthy lymphopenic conditions typically found in severe disease of the immune system (Brass et al. 2014), cytotoxic drug use (Gergely 1999), radiation treatment (Grossman et al. 2015), or other abnormal situations. These scenarios are, however, beyond the scope of our analysis.

Our model has six adjustable parameters, $\gamma_{0}$, $a$, $p$, $\mu_{0}$, $\mu_{1}$ and $K$. The first four are biologically inherent to the mechanism of T-cell homeostasis, and have been measured experimentally in humans and rodents. The last two have to be constrained via parameter sweeps to match relevant experimental observations. Fig. 1(a) illustrates four qualitatively distinct evolution trajectories of $N(t)$ that may arise from simulations of the model in the presence of a decaying thymic export rate $\gamma(t)$ (gray dash-dotted curve). To non-dimensionalize Eqs. 1, 2, we use $a^{-1}$ to rescale $t$ and $K$ to rescale $N$. The qualitative behavior of our model is thus controlled by three independent parameters: $\gamma_{0}a^{-1}K^{-1}$, $(p-\mu_{0})a^{-1}$, and $\mu_{1}\left(p-\mu_{0}\right)^{-1}$. The black dashed curve arises when $\mu_{1}\left(p-\mu_{0}\right)^{-1}<1$. In this case cell proliferation always exceeds cell death, leading to unbounded expansion of the naive T-cell population. This scenario is unrealistic, except perhaps during a period of lymphopenia. For $\mu_{1}\left(p-\mu_{0}\right)^{-1}\geq 1$, cell death is able to balance cell proliferation at a homeostatic carrying capacity $N=N_{\rm ss}(\gamma=0)$, defined by $\mu(N_{\rm ss}(\gamma=0))=p$, as $\gamma\to 0$. As illustrated by the green dotted curve, $N(t)$ rises and asymptotically converges towards $N_{\rm ss}(\gamma=0)$ provided that $\gamma_{0}a^{-1}K^{-1}\ll 1$. We refer to this scenario as being in the “proliferation-driven” regime, given that the cell population is driven to $N_{\rm ss}(\gamma=0)$ primarily by homeostatic proliferation. The model’s behavior makes a transition from proliferation-driven to “thymus-driven” if we increase $\gamma_{0}a^{-1}K^{-1}$. As shown by the blue solid curve, $N(t)$, driven by increased thymic export, overshoots and approaches $N_{\rm ss}(\gamma=0)$ from above as $\gamma(t)\to 0$ asymptotically. Finally, the red dash-dotted curve arises when $(p-\mu_{0})a^{-1}\leq 0$. In this case cell death always exceeds cell proliferation as $\gamma(t)\to 0$, and $N(t)\to N_{\rm ss}(\gamma=0)=0$. As stated earlier, in this paper we focus on scenarios of healthy aging (lymphoreplete) conditions, which immediately rules out the scenarios of unbounded growth (black dashed curve) and complete collapse of the T-cell population (the red dot-dashed curve), effectively constraining our parameters to physiologically reasonable values $\mu_{1}\left(p-\mu_{0}\right)^{-1}\geq 1$ and $(p-\mu_{0})a^{-1}>0$.

We can further quantitatively calibrate the parameter values using experimental measurements in the literature. The constant peripheral proliferation rate $p$ has been measured by Westera et al. (2015) as $0.05$% $\text{day}^{-1}$, or equivalently $p=0.18$ $\text{year}^{-1}$. The basal death rate $\mu_{0}$ can be estimated from the lifespan of T-cells. Based on data from Vrisekoop et al. (2008), De Boer and Perelson (2013) obtain an average naive $\text{CD4}^{+}$ T-cell lifespan of $\sim 5$ years and an average naive $\text{CD8}^{+}$ lifespan of $\sim 7.6$ years. Given the normal $\text{CD4}^{+}\text{:CD8}^{+}$ ratio of 2:1, the average combined naive T-cell clearance rate is $\mu_{0}=\frac{1}{5.9}$ $\text{year}^{-1}$ = $0.17$ $\text{year}^{-1}$. Thymic involution with age can be quantified by measuring the decrease in thymic epithelial volume (Steinmann 1986), based on which Murray et al. (2003) showed that thymic output decreases by an average of $4.3\%$ per year between ages [math] and $100$, implying a decay factor of $a=|\ln(0.957)|\simeq 0.044$. The rate of thymic export has recently been measured for young adults ($20-25$ years old) at $\sim 1.6\times 10^{7}$ trained cells daily, or equivalently $5.8\times 10^{9}$ per year (Westera et al. 2015). Assuming that this rate is $\gamma(t)$ at $t=25$ years, we can back-calculate $\gamma_{0}=(5.8\times 10^{9})\times\left(\frac{100}{33.3}\right)\approx 1.75\times 10^{10}\text{cell exports}/\text{year}$. Note that these values of $p$, $\mu_{0}$, and $a$ satisfy the constraint $\left(p-\mu_{0}\right)a^{-1}>0$ that prevents the T-cell population from completely collapsing.

While direct experimental measurements of $\mu_{1}$ and $K$ are not available in the literature, further inspection of Fig. 1(a) reveals that $\mu_{1}$ and $K$ determine whether thymic export or homeostatic proliferation dominates the evolution of $N(t)$. Through the dimensionless parameters, $\gamma_{0}a^{-1}K^{-1}$ and $\mu_{1}\left(p-\mu_{0}\right)^{-1}$, the time at which $N(t)$ peaks and how fast it declines from the peak vary with changes to the values of $\mu_{1}$ and $K$. Recently, Westera et al. (2015) reported a $52\%$ decrease in total naive T-cell counts between young human adults and elderly individuals, which we can use to quantitatively constrain $\mu_{1}$ and $K$. Let us define individuals of an age between $t=20$ and $30$ years as young adults, and those between $t=70$ and $80$ as the elderly. Assuming that interpersonal heterogeneity unrelated to age averages out over large sample sizes in clinical data, we may evaluate $\bar{N}_{\rm y}=\frac{1}{10}\int_{20}^{30}N(t){\rm d}t$ and $\bar{N}_{\rm o}=\frac{1}{10}\int_{70}^{80}N(t){\rm d}t$ as the average naive T-cell counts respectively for the young and the elderly, as illustrated by the shaded areas under the thymus-domination curve in Fig. 1(a). The relative change in the naive T-cell count between young and elderly adults can thus be evaluated as

[TABLE]

Fig. 1(b) plots $\Delta(\bar{N})$ as a function of $\gamma_{0}a^{-1}K^{-1}$ and $\mu_{1}\left(p-\mu_{0}\right)^{-1}$, with $a=0.044\textrm{ year}^{-1}$ for converting the dimensionless time to years to compute $\bar{N}_{\rm y}$ and $\bar{N}_{\rm o}$. When $\gamma_{0}a^{-1}K^{-1}\lesssim 1$ and $\mu_{1}\left(p-\mu_{0}\right)^{-1}\lesssim 2$, $\Delta(\bar{N})>0$. Note that the homeostatic carrying capacity when $\gamma(t)=0$ is $N_{\rm ss}(\gamma=0)=K\left(\mu_{1}(p-\mu_{0})^{-1}-1\right)^{-1}$. A small $\gamma_{0}a^{-1}K^{-1}$ value represents a relatively low thymic export rate, and the carrying capacity increases rapidly as $\mu_{1}\left(p-\mu_{0}\right)^{-1}\to 1$, both of which make it challenging for thymic output to fill up the T-cell pool to carrying capacity before $\gamma(t)$ considerably decays within $t\sim a^{-1}$. As a result, $N(t)$ does not reach a peak value at a young age and continues increasing into old age. The $\approx 52\%$ decrease in naive T-cell counts reported by Westera et al. (2015) is depicted by the black dotted curve. If we set $\mu_{1}\left(p-\mu_{0}\right)^{-1}=4$, our model can be calibrated to reproduce this decrease in the cell count by choosing $K=10^{10}$ ($\gamma_{0}a^{-1}K^{-1}\simeq 41$ with $\gamma_{0}=1.8\times 10^{10}$ and $a=0.044$). In contrast, $K=10^{12}$ yields $\gamma_{0}a^{-1}K^{-1}\simeq 0.41$, leading to an increase in the cell count ($\Delta(\bar{N})\simeq 0.63$). In between, $K=10^{11}$ results in a moderate decrease in the cell count ($\Delta(\bar{N})\simeq-0.33$). For the rest of the paper, we fix $K=10^{10}$ and $\mu_{1}\left(p-\mu_{0}\right)^{-1}=4$, or equivalently $\mu_{1}=0.04$ given that $p=0.18$ and $\mu_{0}=0.17$, so that the age-related decline of $N(t)$ in our model is consistent with Westera et al. (2015).

Note that there exist two intrinsic timescales in Eq. 1; thymic export decays at a rate $a$, while the homeostatic time scale is controlled by $p$, $\mu_{0}$, and $\mu_{1}$. If homeostasis is much faster than thymic involution, the solution of $N(t)$ will quickly converge to the quasisteady state solution as $\gamma(t)$ evolves. We compare these two solutions in Fig. 2(a), where the quasisteady-state solution is obtained by solving for the steady-state solution $N_{\rm ss}$ of Eq. 1 with fixed $\gamma(t)$ at each time $t$, and $N_{\rm ss}(\gamma(t))$ (black dashed curve) decreases monotonically with age due to the continuous decline of $\gamma(t)$. In contrast, $N(t)$ (blue solid curve) slowly rises from the initial conditions $N(1)=10^{11}$ and does not approach the quasisteady-state level until age $\approx 20$ years. The trajectory of $N(t)$ then overshoots the declining $N_{\rm ss}(\gamma(t))$, reaches a peak value, and reverts course to go after $N_{\rm ss}(\gamma(t))$. However, $N(t)$ never catches up with $N_{\rm ss}(\gamma(t))$ before the latter reaches a steady state of very low cell counts. That $N(t)$ keeps lagging behind $N_{\rm ss}(\gamma(t))$ indicates that the timescale for the full model solution to converge to the steady state is slower than the evolution of the nonautonomous term $\gamma(t)$. The results here suggest that steady-state solutions cannot adequately describe the temporal evolution of the T-cell population in the biologically relevant range of parameter values that we have implemented. It is necessary to numerically compute the time-dependent solutions for the full nonautonomous equation.

Indeed, we find a disparity in the rates at which thymic export decays and the steady state solutions evolve. The latter is provided by the inverse of the eigenvalue of Eq. 1 linearized around $N=N_{\rm ss}(\gamma(t))$. The eigenvalue takes the form $\lambda_{1}=p_{0}-(\mu_{0}+\mu_{1}((3N_{\rm ss}^{2}K^{2}+N_{\rm ss}^{4})/((K^{2}+N_{\rm ss}^{2})^{2}))$. Simulations in Fig. 2(b) show that for the biologically relevant parameter values we have implemented, the cell-population evolution timescale, $|\lambda_{1}|^{-1}$ (red solid curve), is generally longer than the timescale of thymic involution ($a^{-1}\simeq 22.7$ years for $a=0.044$ as denoted by the horizontal black dotted line). Hence the nonautonomous solutions $N(t)$ are expected to lag behind the thymus-driven steady-state solutions $N_{\rm ss}$. For $N(t)$ to be reasonably approximated by $N_{\rm ss}$, the cell population has to evolve much faster than thymic involution, corresponding to the regime of very large $\mu_{1}$, as indicated by the blue dash-dotted curve, where cell death is extremely sensitive to the cell population size.

2.2 Clonotype Abundance Distributions

Quantification of the populations of individual clonotypes would require analysis of models that track the population dynamics of naive T-cells of each TCR type. Assuming the same population dynamics for each T-cell clonotype $i$, which may be appropriate for certain scenarios, the evolution of the expected cell count $n_{i}(t)$ may be deduced from Eq. 1 and take the following generalized form,

[TABLE]

where $\gamma(t)/\Omega$ represents thymic export of naive T-cells of each clonotype (the total thymic export rate normalized by the total number of viable TCR combinations $\Omega$), and $N(t)=\sum_{i}n_{i}(t)$. Within the framework of these “neutral” models, basic qualitative behaviors of T-cell population dynamics have been investigated, particularly for scale-invariant properties that can be studied in a reduced system (Lythe et al. 2016; Desponds et al. 2015). Indeed, the total numbers of T-cell clonotypes $\Omega$ in rodent or human bodies are prohibitively large for direct numerical simulations of the full system using Eq. 4. It is thus common to reduce the full system to a more manageable size with the assumption that the phenomena under investigation are scale-invariant. However, it is sometimes difficult to assert whether a certain property really does not change in a re-scaled system, as nonlinear phenomena, such as the Allee effects, often arise in population dynamics and cast doubt on the scalability of the system. Moreover, some properties, such as the thymic export rate $\gamma(t)$, are naturally scale dependent. It is not always clear how these quantities should be re-scaled in a reduced system, and they have usually been omitted by simplification arguments in previous models, which limits the applicability of these models.

In particular, thymic involution is known to be associated with the age-related loss of T-cell diversity. Without the explicit inclusion of the thymic export rate, such loss of T-cell diversity cannot be properly investigated. To facilitate a more manageable full-system model, we consider a formulation that tracks how the expected number of clones of a given size changes with time. By focusing on clone count rather than the explicit cell count of each distinct clonotype, we are able to effectively reduce the number of tracked variables and thus the dimension of the model. This representation was used by Ewens in population genetics (1972), by Goyal et al. (2015) in the context of hematopoietic stem cell population dynamics, and by Desponds et al. in the context of T-cells (2017). We define $\hat{c}_{k}(t)$ to be the number of clones represented by exactly $k$ naive T-cells in the organism at time $t$:

[TABLE]

where the Kronecker delta function $\delta_{x,y}=1$ when $x=y$ and [math] otherwise. By lumping clonotypes of the same cell count into one single variable $\hat{c}_{k}$, this alternative formulation can efficiently describe changes to the TCR clone diversity in the full system, albeit at the expense of the ability to distinguish each specific clonotype (Morris et al. 2014; Mora and Walczak 2016). Individual clone information is lost, and $n_{i}(t)$ cannot be recovered from $\hat{c}_{k}(t)$ after the transformation in Eq. 5. Nonetheless, the amount of computation can be significantly reduced by truncating $\hat{c}_{k}(t)$ at a reasonably large $k$, as few large clones exist in realistic scenarios, and $\hat{c}_{k}(t)$ for large $k$ is negligible. Letting $c_{0}(t)\equiv\langle\hat{c}_{0}(t)\rangle$ denote the expected number of all possible (thymus-allowed) clonotypes unrepresented in the periphery at time $t$, and $c_{k}(t)\equiv\langle\hat{c}_{k}(t)\rangle$ the expected number of clones of size $k$ at time $t$, a closed set of equations governing the evolution of $c_{k}(t)$ can be derived from Eq. 4 in the mean-field limit,

[TABLE]

where $N(t)=\sum_{i}^{\infty}n_{i}(t)=\sum_{\ell=1}^{\infty}\ell c_{\ell}(t)$. The expected values $c_{k}(t)$ are also called species abundances in the ecology literature. The number of unrepresented clones is $c_{0}=\Omega-\sum_{k=1}^{\infty}c_{k}$, and summing Eq. 6 multiplied by $k$ over $k=1,2,\cdots$ recovers Eq. 1. The mean-field assumption is articulated in terms such as $\mu(\sum_{\ell}\ell\hat{c}_{\ell})\hat{c}_{k}$ that involve higher-order products of $\hat{c}_{k}$ rather than correlations of products of $\hat{c}_{k}$.

We have found (unpublished) that this mean-field approximation breaks down only when $\gamma/\mu<1/\Omega\ll 1$ for which the total population is proliferation driven and the quasistatic configuration is $N\sim K$ and all $c_{k}\sim 0$ except $c_{N}$. Thus, we reasonably assume that $\gamma(t)>\mu/\Omega$ allowing the use of the mean-field equations 6.

In Eq. 6, the terms in the forms of $(\gamma(t)/\Omega)c_{k}$, $pkc_{k}$, and $\mu(N)kc_{k}$ respectively represent the effect of thymic export, homeostatic proliferation and cell death on a T-cell clone already represented by $k$ cells in the peripheral blood. Adding one cell via thymic export or homeostatic proliferation moves one clone from the $c_{k}$-compartment to the $c_{k+1}$-compartment, while the death of one cell shifts one clone from the $c_{k}$-compartment to the $c_{k-1}$-compartment. We approximate the proliferation rate $p$ as a constant, at which rate all cells of all clones of size $k$ replicate via homeostatic proliferation. Proliferation reduces $c_{k}$ and increases $c_{k+1}$. Terms of the form $\mu(N)kc_{k}$, where the IL-7 regulated death rate $\mu(N)$ is given by Eq. 2, reduce $c_{k}$ and increase $c_{k-1}$.

For a healthy aging adult, the TCR repertoire is mostly comprised of small clones with the probability of finding large clones decreasing with clone size $k$. To numerically solve Eq. 6, we thus truncate the model at a maximum clone size $M\gg 1$, beyond which the probability of finding a clone is assumed negligible. For our implementation of the truncation, please see Appendix A. In Fig. 3(a) we examine the effect of the truncation clone size $M$, showing sufficient convergence of $c_{10}$ at $t=40$ and $70$ to fixed values when $M\gtrsim 30$, which indicates that further inclusion of clones beyond $c_{30}$ has little effects on the solution for $t\lesssim 70$ years. For numerical simulations of Eq. 6 in this paper, we set $M=200$ to ensure minimal truncation errors.

Fig. 3(b) shows the temporal evolution of $c_{k}(t)$ for $k=2$, $19$, and $59$. As $k$ increases, the overall magnitude of the $c_{k}(t)$ curve decreases, and the age at which $c_{k}(t)$ peaks increases. For example, $c_{2}(t)$ peaks around $t\lesssim 20$ years, and there are many fewer clones of exactly two copies at old ages than at young ages. In contrast, $c_{19}(t)$ peaks around age $55$, and the numbers of clones that have exactly $19$ copies are roughly the same between old and young ages, whereas the number of clones that have exactly $59$ copies ($c_{59}(t)$) keeps increasing into old ages.

The relatively earlier decline of $c_{k}(t)$ with smaller $k$ is expected, considering that rare clones are introduced into the peripheral circulation primarily by the thymus, which started to involute after birth. With increasing $k$, the influence of thymic export on $c_{k}(t)$ decreases, whereas the dependence on homeostatic proliferation increases. Recalling that the rate of thymic involution is faster than the time scale for homeostasis to drive the clonal population towards equilibrium, the fast decline of the rare clone population leaves room for larger clones to expand.

To accompany the steady state $N_{\rm ss}$, we compute analogous fixed-$\gamma_{0}$ steady state values of the full system, $c_{k}^{\rm ss}$, in Appendix B. The steady states satisfy $c_{k}^{\rm ss}\to 0$ as $\gamma_{0}\to 0$ for all $1\leq k\leq M$. We further show that in spite of the fact that $c_{k}^{\rm ss}\to 0$, Eq. 6 asymptotically yields a positive total cell count $N=\lim_{M\to\infty}\sum_{k=1}^{M}kc_{k}^{\rm ss}>0$ as $M\to\infty$, qualitatively consistent with Eq. 1. Moreover, we prove in Appendix C that solutions $c_{k}(t)$ of the full nonautonomous system satisfy $c_{k}(t)\to 0$ for all $k\leq M$, with arbitrarily large $M$, as $t\to\infty$. This result is completely independent of the assumed functional forms of the proliferation and death rates, suggesting that manipulation of homeostatic regulatory mechanisms cannot prevent the extinction of small T-cell clones caused by decaying $\gamma(t)$. We thus conclude that thymic involution dictates the age-related decline of the TCR diversity of the naive compartment.

2.3 Diversity of the Naive T-cell Repertoire

By computing the functions $c_{k}$ that track the number of clones consisting of $k$ cells, we should have sufficient information to evaluate the variation in TCR structural diversity over a lifetime. Expected TCR structural diversity or “richness” is the total number of distinct clones present in the immune compartment, for which we define a threshold TCR richness diversity,

[TABLE]

where $q\in\mathds{N}$ is a lower threshold, so that the quantity $R_{q}(t)$ represents the number of clones of size at least $q$ present in the immune compartment at time $t$. Situations in which such a $q$-dependent threshold arise may include consideration of immune surveillance, in which small clones may evade detection.

As shown in Fig. 4(a), $R_{q}(t)$ increases at young ages, peaks at a mature age, and declines afterwards. For our previous parameter values, the peak age of $R_{1}(t)$ is approximately $t\sim 16$. Higher $q$ lead to older peak ages of $R_{q}(t)$, consistent with the results in Fig. 3(b), in which the number of larger clones peaks at older ages.

To compare $R_{q}(t)$ between the elderly and young, we adopt the same criterion as with total cell counts and compute window-averaged values of $R_{q}(t)$ between ages $20$ and $30$ for the young and between ages $70$ and $80$ for the elderly. By defining $\bar{R}_{\rm y}(q)\equiv\frac{1}{10}\int_{20}^{30}R_{q}(t){\rm d}t$, $\bar{R}_{\rm o}(q)\equiv\frac{1}{10}\int_{70}^{80}R_{q}(t){\rm d}t$, we quantify the loss of richness by computing its relative change.

[TABLE]

Using the same parameter values as in Fig. 4(a), we plot $\Delta(\bar{R}_{q})$ with respect to $\mu_{1}$ and $q$ in Fig. 4 (b),(c). In Fig. 4(b), $\Delta(\bar{R}_{q})$ decreases monotonically with increasing $\mu_{1}$, suggesting that upregulated death rate exacerbates the age-related loss of richness, and the impact is more significant for larger $q$. Fig. 4(c) shows that when $K=10^{10}$, $\Delta(\bar{R}_{q})<0$ for $q\leq 4$. This decreasing trend of $R_{q}$ generally agrees with the loss of diversity observed in recent experiments where measurements were available across multiple ages (Qi et al. 2014; Britanova et al. 2014). For $q=5,6$, $\Delta(\bar{R}_{q})\approx 0$, and $R_{q}$ is nearly unchanged between youth and advanced age. For $q\geq 7$, $\Delta(\bar{R}_{q})>0$, indicating higher $R_{q}$ at older ages. Generally, the lifetime decrease in $R_{q}(t)$ occurs with small $q$, whereas for large $q$, the trend is reversed, in agreement with our discussion of Fig. 3(b) and Fig. 4(a) regarding peak ages. This phenomenon indicates that loss of diversity is primarily due to the extinction of rare clones, which is consistent with the observation made by Naylor et al. (2005). In contrast, the number of larger clones increases over time, leading to the lifetime increase to $R_{q}(t)$ at higher $q$.

Recent TCR-$\beta$ sequencing studies have attempted to estimate the change in the repertoire richness of the naive T-cells with age. Despite the difference in orders of magnitude regarding the total number of circulated naive T-cell clones, these studies agreed quantitatively in the ratio of the age-related loss of richness. For example, Britanova et al. (2014) estimated $\sim 7\times 10^{6}$ clonotypes in youth (ages $6-25$), and $\sim 2.4\times 10^{6}$ in aged individuals (ages $61-66$), a roughly $66$% drop from the youth figure. Similar measurements were also reported by Qi et al. (2014), in which a two-to-five-fold decline (i.e., a $50$% – $80$% drop) between youth (ages $20-35$) and advanced age (ages $70-84$) was observed. These results are quantitatively consistent with our computation of $\Delta(\bar{R}_{1})$ for $K=10^{10}$ – $10^{11.5}$ and $0.03\leq\mu_{1}\leq 0.05$ in Fig. 4(d), whereas the decline of $R_{q}$ for $q\geq 2$ is not as pronounced as in these experimental observations.

Also note that the loss of clonal richness is more severe than the decrease in the total cell count between young and aged individuals. In Fig. 4(a) $\Delta(\bar{R}_{1})$ changes between $\sim-66\%$ and $\sim-76\%$ for $0.03\leq\mu_{1}\leq 0.05$ and $K=10^{10}$. In contrast, Fig. 1(b) shows that for the same parameter range, $\Delta(\bar{N})$ varies from $\sim-30\%$ to $\sim-62\%$. However, the figures also reveal that richness is relatively less sensitive to changes to the cellular death rate, compared to the total cell count. This outcome reflects the fact that homeostatic cellular death is uniformly random across the entire naive T-cell population. The drop in richness is due to cell death within small clones that drives these clones to extinction, as observed by Naylor et al. (2005). Increases to the cellular death rate do not cause as much additional clonal extinction as they do additional cellular extinction, as many surviving clones are too large to wipe out by the death of a few cells.

2.4 Sampling Statistics

Considering that naive T-cell richness is often assessed via small blood samples, let us next use the same framework to examine the relation between the detected clone sizes in small samples and the true clone sizes in the full organism. As before, denote by $N$ the total number of naive T-cells in the human’s immune compartment, and $Y\leq N$ the number of cells collected during sampling from among the $N$ total. We assume that the $N$ total cells consist of $R$ distinct clones, which we number from $1$ to $R$. In this section, we denote by $c_{k}^{N}$ the mean number of clones of size $k$ from among the $N$ total cells in the full organism (denoted by $c_{k}$ in the previous simulations), and by $c_{k}^{Y}$ the mean number of clones of size $k$ in the sampling of $Y$ cells taken from the $N$ total cells. Then the expectation of $c_{k}^{Y}$, denoted by $\mathbb{E}[c_{k}^{Y}]$, is,

[TABLE]

where $P\left(c_{k}^{Y}=j\right)$ represents the probability that there are precisely $j$ clones of size $k$ in the sampling. Then $\mathbb{E}[c_{k}^{Y}]$ may be expressed explicitly in terms of the $c_{k}^{N}$ as:

[TABLE]

(See Appendix D for the detailed proof.) The collection of expressions given by Eq. 10 for $k=1,2,\cdots,R$, yields a linear system of equations solvable for $c_{k}^{N}$, using sampled data for the quantities $\mathbb{E}[c_{k}^{Y}]$. More specifically, if we define the vectors $\mathbf{\widehat{E}}:=(\mathbb{E}[c_{1}^{Y}],\mathbb{E}[c_{2}^{Y}],\cdots,\mathbb{E}[c_{R}^{Y}],)$ and $\mathbf{E}:=(c_{1}^{N},c_{2}^{N},\cdots,c_{R}^{N})$, Eq. 10 can be written as $\mathbf{\widehat{E}}=\mathbf{AE}$, where $\mathbf{A}$ is a constant matrix that has non-zero elements only in the upper triangle, with non-zero diagonal entry $\frac{1}{\binom{N}{Y}}\binom{N-k}{Y-k}$ in position $(k,k)$. The equation can always be solved uniquely for $\mathbf{E}$ given $\mathbf{\widehat{E}}$. Thus the full size distribution $\mathbf{E}$ can be uniquely reconstructed from the expected mean sample size distribution $\mathbf{\widehat{E}}$ measured experimentally, provided that the latter can be reliably estimated through a sufficient number of repeated samplings.

In Fig. 5(a), we use Eq. 10 to compute $\mathbb{E}[c_{k}^{Y}]$ from simulated $c_{k}^{N}$, comparing the predicted sampling results for varying choices of $Y$. The results indicate that each decrease by one order of magnitude to the sample size results in a decrease by roughly the same order of magnitude to the predicted diversity. Thus, diversity predictions vary with sample size, and small samples do not result in accurate measurements of diversity.

In Fig. 5(b) we examine how sampling may affect the diagnosis of the age-related TCR richness decline $\Delta\left(\bar{R}_{q}\right)$ defined in the previous subsection. We find $\Delta\left(\bar{R}_{q}\right)$, which is negative, increasing with decreasing sampling fraction $f$, revealing that sampling causes an underestimate of the richness decline. As previously discussed, the decline of TCR richness at old ages is primarily due to the extinction of small clones. Since small clones often evade detection during sampling, their extinction is largely unaccounted for, leading to lessened reduction of the richness measure. When $f$ is very small, most of the small clones have escaped detection; thus decreasing $f$ further does not change $\Delta\left(\bar{R}_{q}\right)$. Moreover, we note that $\Delta\left(\bar{R}_{1}\right)$, which is the most straightforward measure for age-related loss of TCR richness, changes from $-73\%$ for the full sample, to $-59\%$ for a sampling fraction $f\leq 10^{-3}$, which is close to the value of $\Delta\left(\bar{R}_{3}\right)$ for the full sample. This reaffirms our discussion in the previous subsection that a threshold $q>1$ may arise during the process of sampling. The results here indicate that when only a small fraction of a T-cell population is used to measure $\Delta\left(\bar{R}_{1}\right)$, clones fewer then three copies largely evade detection, yielding a result equivalent to $\Delta\left(\bar{R}_{3}\right)$ of the full sample, which underestimates the actual decrease of the TCR richness.

3 Discussion

We have formulated a model of lifetime human naive T-cell population dynamics, which traces T-cell lineages on the level of individual clones. It accounts for exponentially decaying lifetime thymic export, a constant rate of cellular proliferation, and variable cellular death rate that adjusts to present cell counts and availability of survival resources. It depicts the generation of the naive T-cell pool in early life via thymic export, and long-term maintenance of the population via peripheral turnover after thymic export has waned. Values of most of the model’s parameters can be found in previous literature, while the few exceptions are obtained by fitting some basic results of the model, such as age-related T-cell loss, to previous experiments. Our analysis serves two important purposes: to map the thymic machinery, identifying which components do and do not contribute to age-related cellular loss, and then to interpret the nuanced role of that cellular loss in immunosenescence.

First, we have found that if thymic export is assumed to decay exponentially to zero, then all compartments $c_{k}(t)$ (with $1\leq k\leq M$) deplete as $t\to\infty$, independent of essentially any restrictive assumptions about the homeostatic proliferative mechanism in the periphery. Concretely, for any choice of proliferation and death rates $p(N),\mu(N)$, that satisfy $p(0),\mu(0)>0$ and the choice $\gamma(t)=\gamma_{0}e^{-at}$ with $\gamma_{0},a>0$, there exists a sufficiently small $\delta>0$ guaranteeing $c_{k}(t)\to 0$ as $t\to\infty$ for all $1\leq k\leq M$, provided that $\sum|c_{k}(1)|\leq\delta$. Although this result only guarantees that trajectories $c_{k}(t)$ started sufficiently close to zero converge to zero, simulation indicates that the basin of attraction to this “zero state” is actually quite large. In fact, for the typical initial conditions used throughout this paper, simulation suggests convergence of all compartments $c_{k}$ to zero in infinite time. Although it takes an extremely long time to deplete all $c_{k}$ compartments for $1\leq k\leq M$, the initial phase of this process can still cause significant loss of T-cell diversity in aging individuals within a human lifespan. Most importantly, we find that the T-cell loss driven by exponentially-diminishing thymic export alone is robust against any assumptions about the homeostatic proliferative mechanism in the periphery, as this outcome is universal for all functional forms of $p(N),\mu(N)$; even a particularly strong homeostatic mechanism (say, one with $p(0)\gg\mu(0)$) cannot rescue a plunging diversity. This, in turn, suggests that in searching for treatments of age-induced loss of diversity, efforts should be directed at the thymus, in particular to maintaining thymic productivity into advanced age.

Moreover, we compare the real-time simulations and the quasisteady-state solutions of the total cell count, as well as the number of distinct clones, over the course of age-related thymic output erosion. We find that our simulation results keep lagging behind the quasi steady state solutions, suggesting that the erosion time scale of thymic output is faster than the time scale for the population dynamics to relax towards a steady state. Mathematically, this result reveals that the evolution of the T-cell population within the human lifespan is a rather dynamical phenomenon, which may not be well-described by quasistatic solutions, requiring evaluation of the fully nonautonomous system. Biologically, our results indicate that the loss of T-cell diversity is a delayed response to thymic involution, and assessment of thymic function may predict the health of the immune system.

Although peripheral division cannot salvage the T-cell population on a long time scale, higher basal proliferation rates may at least delay the erosion of the T-cell compartment, sustaining acceptable effectiveness of the immune system within the human lifespan (Naylor et al. 2005). We assumed a constant lifetime rate of cellular proliferation, but alternative research suggests that proliferation rates may increase with age (Naylor et al. 2005). In light of this finding, we briefly inspect the effect of increased proliferation rates at advanced ages on cellular and clonal loss by modifying $p(N)$ and $\mu(N)$ in Eq. 6. For simplicity, we take the death rate to be constant ($\mu(N)=\mu_{0}>0$), and adopt a logistic growth rate, $p(N,t)=p(t)(1-N/K)$, where a discrete increase in the proliferation rate is incorporated in $p(t)=p_{0}(1+rH(t-T))$, with $p_{0}>0$ the early-life basal cellular proliferation rate, and $H(t)$ the Heaviside function, with $T$ the age at which the rate increases. The constant $r$ specifies the increase to the proliferation rate. (Full simulation details are given in the caption of Fig. 6.) By varying $r$, simulation under these alternate hypotheses indicates that increased basal proliferation rates do lead to notably higher total cell counts (Fig. 6(a)), but have little effect on diversity (Fig. 6(b)). These results further affirm that expansion of peripheral proliferation is unlikely to rescue the eroding naive T-cell diversity, despite the increased cell count. If diversity loss is the main cause of immunosenescence (still a debatable topic in the medical community), peripheral proliferation may not be the sensible target of treatments.

The increased $N(70)$ and nearly unchanged $R_{1}(70)$ in Fig. 6 imply that the decline of T-cell diversity at old age may appear more dramatic if the diversity is measured in terms of frequency of distinct TCR sequences among the cycling cells, which corroborates the explanation that an increase of proliferation rate at old age leads to a sharp decrease of T-cell diversity (Naylor et al. 2005). Previous models have shown that even sharper decline of T-cell diversity can be induced by fitness selection, where certain clonotypes increase their fitness at old age possibly due to higher avidity to self-antigens (Johnson et al. 2012,  2014; Goronzy et al. 2015a).

Although the boosts to the total cell count through artificial expansion of the proliferative mechanism are unable to replenish the declining TCR diversity in the naive T-cell pool, it is possible that the impact is less severe than the decaying richness would have indicated, considering that most of the extinct clones are originally small clones, which may be much less effective than larger clones. In this regard, the viability of treating immunosenescence by expanding peripheral proliferation depends on the elucidation of the T-cell pool’s effectiveness clone size–that is, the size a clone must have attained to effectively guarantee activation of the clone when its cognate antigen infiltrates the organism. The effectiveness clone size is intrinsically linked to true functional TCR diversity; if we can identify a threshold integer $q^{*}$, such that clones of size at least $q^{*}$ are reliably activated in the presence of their cognate antigen(s), but that smaller clones are not, then $R_{q^{*}}(t)$ is naturally the most useful measure of diversity, because it accounts for precisely those clones actively participating in the adaptive immune mechanism. The larger the “correct” choice of $q^{*}$ is, the more effective treatments to boost cellular proliferation in the periphery will be. Our model directly yields the number of clones of a particular size, making it straightforward to include or exclude clones below a certain cell count, should such a threshold exist and be identified.

The effective clone size is also significant to the question of whether diversity loss is the driving factor in immunosenescence. Using the parameter values that we found in the literature, $R_{q}(t)$ decreases for $q\leq 4$ from youth to advanced age, stays nearly constant for $q=5,6$, and increases for $q\geq 7$. The extinction of small clones allows the surviving clones to expand in size, leading the richness of large clones to increase at old ages. If the minimal size for a T-cell clone to effectively respond to antigens is large, the diversity of such “effective” clones may actually increase with age, strengthening the immune response. Therefore, either the minimal clone size required for effective immune response is low, or the weakened immune response at old ages is caused primarily by other mechanisms. For example, functional deficiencies acquired by naive T-cells in aging are one possible alternative cause of the weakened immune response. Such functional deficiencies have been studied heavily in mouse models, but research in humans is still lacking (Appay and Sauce 2014). Diminished naive T-cell effector responsiveness and proliferative capacity have been observed in aged mice (Moro-García et al. 2013). It is possible that similar changes occur in humans. Conversely, experiments on mice have directly shown that loss of TCR diversity does have an actively detrimental effect on immune responsiveness (Yagger et al. 2008), supporting the notion that loss of TCR diversity as a significant contributor to immunosenescence.

Our model illustrates the feasibility of several different scenarios, in which loss of diversity contributes to immunosenescence on drastically different levels. There is clearly a strong need to investigate the effects of both age-related structural diversity loss and T-cell functionality loss in human subjects in vivo, to better understand the causes of immunosenescence. Moreover, our model indicates that the effectiveness clone size and crossreactivity in vivo are valuable pieces of missing information, the elucidation of which would allow for the identification of effective options to treat immunosenescence.

4 Summary and Conclusions

We have simulated the time evolution of the functions $c_{k}(t)$, which represent the number of naive T-cell clones of size $k$ present in a human’s immune compartment at time $t$. We determined that under essentially any realistic assumptions about homeostatic proliferation and death, all clones deplete in infinite time if thymic export is assumed to decay exponentially. This implicates thymic export as a fundamental cause of age-associated diversity loss. We simulated our model under the assumption that a carrying capacity is regulated by homeostatic proliferation and death through $N$-dependent rates. We found that the manipulation of homeostatic proliferation and death rates, which may notably raise the carrying capacity and thus the total cell count, was unable to save falling diversity as an individual ages. It affirms the vital role of thymic output in age-related diversity loss, and indicates that boosting the proliferation rate is unlikely an effective solution. However, if only clones of large size are sufficiently effective in the immune response, boosting proliferation rates might raise average clone sizes and help to mitigate the effects of lost diversity. We simulated “threshold richness diversity”, $R_{q}(t)$, which counts the total number of clones of size $q$ or larger. We found that by increasing $q$, the trajectory of $R_{q}(t)$ changes from decreasing to increasing over a human lifetime. From this trend, we concluded that if only large clones are effective, the effective richness would actually increase with age, suggesting that it is important to identify the minimal effective clone size in order to determine whether the loss of TCR diversity is the primary driving mechanism of the immune dysfunction seen in advanced age. Lastly, we derived a one-to-one mapping between the full-sample diversity $c_{k}^{N}$ of $N$ cells and the expected measurement of diversity $\mathbb{E}[c_{k}^{Y}]$ in samples of $Y$ cells. We found that the probability of detecting small clones shrank significantly with small sample sizes, which could potentially skew small-sample statistics. In particular, we show that small samples tend to underestimate the age-related loss of T-cell richness diversity. Our formulation provides a rigorous method for accurately inferring the statistical distribution of clonal sizes from small-sample measurements.

Acknowledgements

This work was supported by the NIH (SL, R56HL126544), the NSF (TC and SL, DMS-1516675 and DMS-1814364) and the Army Research Office (YLC, W1911NF14-1-0472).

Appendix A Implementation of Numerical Truncation

The most straightforward way to truncate Eq. 6 at $k=M$ is to neglect the exchange terms between $c_{M}$ and $c_{M+1}$, assuming a negligible contribution for $k>M$ and essentially imposing a “no-flux” boundary condition. This leads to the following equation for the boundary term $c_{M}(t)$:

[TABLE]

This formulation, however, introduces a truncation error in Eq. 1 if we define $N=\sum_{k=1}^{M}kc_{k}$. The neglected terms leave a small loss of total cell count in ${\rm d}N/{\rm d}t$. An alternative implementation of the truncation is adding these small loss terms to the boundary equation:

[TABLE]

thus preserving the total cell count $N$. However, for Eq. 12 the truncation error shows up in the total number of clonal types $\Omega=\sum_{k=0}^{M}c_{k}$, as the terms added to Eq. 12 to preserve $N$ artificially introduce new clonal types into the model. In contrast, $\Omega$ is preserved with the implementation of Eq. 11. If $M\to\infty$, the truncation errors for both implementations go to zero at $\sim 1/M$, and the two implementations become equivalent. Assuming sufficiently large $M$, the truncation errors can be negligible in the context of $\gamma(t)>0$, or have minimal cumulative effects within a limited duration, such as a human lifetime, on which our investigations in this paper have primarily focused.

In this paper, we adopt, for simplicity, Eq. 11 to numerically truncate Eq. 6. Note that this choice may seem “natural” if one regards $M$ as the carrying capacity, making it reasonable for $c_{M}$ to have zero proliferation rate. However, the full mechanisms associated with the carrying capacity are far more sophisticated than simply eliminating the proliferation of $c_{M}$. Not only should the proliferation rate of $c_{M}$ go to zero, the proliferation rate of the other $c_{k}$ should also have a $k$ dependence. The $k$ dependence may be weak for small $k$, but as $k\to M$, the proliferation rate should attenuate significantly. The probability that $c_{k\to M}$ will proliferate should be very small, as it is highly likely that there exist other smaller clones to push the total cell count up to the carrying capacity, prohibiting further proliferation. The $k$-dependent proliferation rate will yield a natural truncation threshold at the carrying capacity. However, such a sophisticated $k$-dependence of the proliferation rate is beyond the scope of this paper. Our assumption here is simply that the truncation errors introduced by Eq. 11 are numerically negligible and not biologically significant.

Appendix B Steady States of the Autonomous Equations

If we fix $\gamma(t)=\gamma_{0}$, Eqs. 1, 6, and 11 become autonomous and admit the following steady state solution,

[TABLE]

where $N_{\rm ss}$ is the total population at steady state, given by the unique positive root of the cubic,

[TABLE]

When $\gamma_{0}=0$, $c(N;0)$ has three real roots, $N=0$, $\pm\sqrt{((p-\mu_{0})K^{2})/(\mu_{0}+\mu_{1}-p)}$. The positive steady state solution, which we denote by $N_{\rm ss}(0)$, is stable, and the zero solution unstable, under the parameter restrictions described in Section 2.1. We now demonstrate that even though Eqs. 13, 14 indicate that each $c_{k}^{\rm ss}\to 0$ as $\gamma_{0}\to 0$, the quantity $\lim_{M\to\infty}\sum_{k=1}^{M}kc_{k}^{\rm ss}$ converges to a positive value qualitatively consistent with $N_{\rm ss}(0)$ as $\gamma_{0}\to 0$.

Proposition B: The steady state solutions $c_{k}^{\rm ss}$, as given in Eqs. 13, 14, satisfy,

[TABLE]

Proof. We seek to derive upper and lower bounds, $U(\gamma_{0})$, $L(\gamma_{0})$, which satisfy,

[TABLE]

for $\gamma_{0}>0$, and $\lim_{\gamma_{0}\to 0}U(\gamma_{0})\geq\lim_{\gamma_{0}\to 0}L(\gamma_{0})>0$. We first establish two small results, which will be used later on:

Proposition B1: For $\mu=\mu(N_{\rm ss}(\gamma_{0}))$, $\lim_{\gamma_{0}\to 0}\frac{{\rm d}\mu}{{\rm d}\gamma_{0}}>0$.

Proof. Recalling that $\mu=\mu(N_{\rm ss}(\gamma_{0}))=\mu_{0}+\mu_{1}(N_{\rm ss}(\gamma_{0})^{2}/(N_{\rm ss}(\gamma_{0})^{2}+K^{2}))$, we have:

[TABLE]

where we computed the derivative $\frac{{\rm d}N_{\rm ss}}{{\rm d}\gamma_{0}}$ implicitly from the expression $c(N_{\rm ss}(\gamma_{0});\gamma_{0})=0$. From the explicit form $N_{\rm ss}(0)=\sqrt{(p_{0}-\mu_{0})K^{2}/((\mu_{0}+\mu_{1})-p_{0})}$, we have:

[TABLE]

Proposition B2: For $f(p/\mu(N_{\rm ss}(\gamma_{0}));\gamma_{0})=\frac{\gamma_{0}}{p\Omega}\left(1-\frac{p}{\mu(N_{\rm ss}(\gamma_{0}))}\right)^{\frac{-\gamma_{0}}{p\Omega}-1}$, $\lim_{\gamma_{0}\to 0}f(p/\mu(N_{\rm ss}(\gamma_{0}));\gamma_{0})>0$.

Proof. We write the function $f(p/\mu(N_{\rm ss}(\gamma_{0}));\gamma_{0})$ as a product of two functions as follows:

[TABLE]

We define $A_{0}=\lim_{\gamma_{0}\to 0}A(\gamma_{0})$ and $B_{0}=\lim_{\gamma_{0}\to 0}B(\gamma_{0})$, and compute $A_{0}$ and $B_{0}$:

[TABLE]

where we used that $\mu(N_{\rm ss}(\gamma_{0}))\to p$ as $\gamma_{0}\to 0$. From Proposition B1, $\lim_{\gamma_{0}\longrightarrow 0}\frac{{\rm d}\mu}{{\rm d}\gamma_{0}}>0$, and a similar computation shows that $\lim_{\gamma_{0}\to 0}\frac{{\rm d}^{2}\mu}{{\rm d}\gamma_{0}^{2}}\in\mathbf{R}$. Thus, $\ln(A_{0})\in\mathbf{R}$, and $A_{0}>0$. Now,

[TABLE]

Thus, $\lim_{\gamma_{0}\to 0}\frac{\gamma_{0}}{p\Omega}\left(1-\frac{p}{\mu(N_{\rm ss})}\right)^{\frac{-\gamma_{0}}{p\Omega}-1}=A_{0}B_{0}>0$.

We now resume the proof of Proposition B. We first derive upper and lower bounds on the term $c_{1}^{\rm ss}$, to simplify calculations. From the nonnegativity of the parameters and coefficient functions, and the form in Eq. 13, $c_{1}^{\rm ss}\leq\gamma_{0}/\mu_{0}$, independent of $M$. To derive an $M$-independent lower bound on $c_{1}^{\rm ss}$, we observe that the sum in the denominator of Eq. 13 satisfies,

[TABLE]

and that the sum on the right above is the $M$-th Taylor polynomial, $S_{M,\gamma_{0}}$, for the function $f(x;\gamma_{0})=\frac{\gamma_{0}}{p\Omega}\left(1-x\right)^{\frac{-\gamma_{0}}{p\Omega}-1}$ expanded around $x=0$ and evaluated at $x=\frac{p}{\mu(N_{\rm ss}(\gamma_{0}))}$. The function $f(x;\gamma_{0})$ is analytic in $x$ away from $x=1$, and in particular, the $S_{M,\gamma_{0}}$ increase monotonically to $f(p/\mu(N_{\rm ss}(\gamma_{0}));\gamma_{0})$. It follows that,

[TABLE]

and thus that $c_{1}^{\rm ss}\geq\gamma_{0}/(pf_{\gamma_{0}}+\mu_{0}+\mu_{1})$. After using the $c_{1}^{\rm ss}$ bounds in the expression for $c_{k}^{\rm ss}$, we have:

[TABLE]

Now we let $L(\gamma_{0})=\frac{p\Omega}{pf_{\gamma_{0}}+\mu_{0}+\mu_{1}}f_{\gamma_{0}}$ and $U(\gamma_{0})=\frac{p\Omega}{\mu_{0}}f_{\gamma_{0}}$. From Proposition B2, $\lim_{\gamma_{0}\to 0}f_{\gamma_{0}}>0$, so $\lim_{\gamma_{0}\to 0}L(\gamma_{0}),\lim_{\gamma_{0}\to 0}U(\gamma_{0})>0$, and Proposition B follows.

Appendix C Convergence and Stability of $c_{k}$ when $\gamma(t)\to 0$

In this section we will prove that solutions $c_{k}$ to our ODE system initialized sufficiently close to $\mathbf{\overrightarrow{0}}$ converge to $\mathbf{\overrightarrow{0}}$ as $t\to\infty$. Denote by (P) the “perturbed” ODE system given by Eqs. 6, 11, with $\gamma(t)=\gamma_{0}e^{-at}$, and by (U) the “unperturbed” ODE system resulting from the alternate choice $\gamma(t)\equiv 0$. For the sake of generality, we omit previous assumptions about the form of the functions $p(N),\mu(N)$, except that $p(0),\mu(0)>0$. Additionally, in this section, we regard the term $N$ that appears in the ODEs as $\sum_{k\geq 1}kc_{k}$ instead of its own variable, and thus do not explicitly include Eq. 1 in our analysis as in Appendix B. Note that the residual $N-\sum_{k\geq 1}kc_{k}\to 0$ as $M\to\infty$. We begin by noting that the unperturbed system (U) has steady-state $c_{k}^{U}(t)\equiv 0$ for $k\geq 1$. To analyze the stability of this steady state, we consider the linearization of (U) around this steady state, which is represented by the $M\times M$ matrix we call $\mathbf{L_{U}}$ ($\mathbf{L_{U}}=(l_{ij})_{1\leq i,j\leq M}$). The components $l_{ij}$ of $\mathbf{L_{U}}$ are given explicitly by:

[TABLE]

Although the matrix is tridiagonal, it is high-dimensional, and thus its eigenvalues cannot be computed analytically. However, we may nevertheless demonstrate that all eigenvalues possess strictly negative real part, indicating that the zero solution is asymptotically stable. To do this, we use Gershgorin’s circle theorem to show that if there exists an eigenvalue $\lambda\in\mathbf{C}$ satisfying $\Re(\lambda)\geq 0$, then $\lambda=0$. We then verify that $\lambda=0$ is never an eigenvalue of $\mathbf{L_{U}}$, by directly demonstrating that $\mathbf{L_{U}}$ has linearly independent rows.

Proposition C: All eigenvalues $\lambda\in\mathbf{C}$ of the matrix $\mathbf{L_{U}}$ satisfy $\Re{\lambda}<0$, so that the zero-solution of (U) is asymptotically stable.

We first apply Gershgorin’s circle theorem to the columns of the matrix $\mathbf{L_{U}}$ to conclude that all eigenvalues $\lambda\in\mathbf{C}$ of the truncated system (finite $M$) are contained within the following union of disks:

[TABLE]

where we have used the fact that $\{\lambda\in\mathbf{C}:|\lambda+D|\leq D\}\subset\{\lambda\in\mathbf{C}:|\lambda+(D+\epsilon)|\leq D+\epsilon\}$ for $D,\epsilon>0$. Given the assumption that $p(0),\mu(0)>0$, each of these disks is tangent to the line $\Re\lambda=0$ at $\lambda=0$, and otherwise lies entirely in the half plane $\Re\lambda<0$. Thus, $\mathbf{L_{U}}$ can only possess an eigenvalue $\lambda$ satisfying $\Re\lambda=0$ if $\lambda=0$ is itself an eigenvalue. We next verify that $\lambda=0$ is never an eigenvalue of $\mathbf{L_{U}}$ directly, by establishing the linear independence of the rows of $\mathbf{L_{U}}$.

Let us assume that there exist scalars $a_{1},a_{2},\ldots,a_{M}$, such that $\sum_{j=1}^{M}a_{j}\left(l_{ij}-0\right)=0$ for all $1\leq i\leq M$. Hence a normalized vector $\mathbf{a}=\left(a_{1},a_{2},\ldots,a_{M}\right)$ represents the eigenvector of the zero eigenvalue. For $i=1$, we find that $2a_{2}\mu(0)-a_{1}(p(0)+\mu(0))=0$, so that $a_{2}=2^{-1}\mu(0)^{-1}(p(0)+\mu(0))a_{1}$. By moving on to larger $i$, we can recursively derive $a_{i}=\Theta_{i}a_{1}$ for all $2\leq i\leq M$ with a proportional constant coefficient $\Theta_{i}$. Moreover, $\sum_{i=1}^{M}\sum_{j=1}^{M}a_{j}l_{ij}=-a_{1}\mu(0)=0$, leading to $a_{1}=0$ given that $\mu(0)>0$. If $a_{1}=0$, $\mathbf{a}\equiv 0$, and a non-zero eigenvector does not exist, implying that zero is not among the eigenvalues of the $M\times M$ matrix $\mathbf{L_{U}}$. We thus conclude that all eigenvalues $\lambda$ of the matrix $\mathbf{L_{U}}$ satisfy $\Re(\lambda)<0$, and the zero-solution of (U) is asymptotically stable for Eq. 6 truncated using Eq. 11 at an arbitrarily large $M$. Note that the proof in Eq. 17 does not hold if we use the alternative truncation formula Eq. 12. By forcing all cells to remain below the truncation threshold $M$, it is not possible for all $c_{k}$ to go to zero with a finite $M$. For the alternative truncation, the stable steady state solution is $c_{k}=2N_{\rm ss}/(M(M+1))$, which nevertheless goes to zero as $M\to\infty$.

We next proceed to demonstrate that the uniform asymptotic stability of the zero-solution ($c_{k}^{U}(t)\equiv 0$ for $k\geq 1$) of the unperturbed system (U) confers a similar notion of “stability” on the perturbed system (P). In particular, the uniform asymptotic stability of the system (U), in conjunction with the exponential decay of the function $\gamma(t)$, implies that solutions of the perturbed system (P) also converge to zero in magnitude, in a sense to be made more precise later on. Here let us simplify our notation by writing (U) as $\mathrm{d}\mathbf{c}/\mathrm{d}t=\mathbf{f}\left(\mathbf{c}\right)$, where $\mathbf{c}\equiv\left(c_{1},c_{2},\ldots,c_{M}\right)$. The autonomous term $\mathbf{f}\left(\mathbf{c}\right)$ consists of cell proliferation and death. Correspondingly we express (P) as $\mathrm{d}\mathbf{c}/\mathrm{d}t=\mathbf{f}\left(\mathbf{c}\right)+\mathbf{g}\left(t,\mathbf{c}\right)$, where the nonautonomous term $\mathbf{g}\left(t,\mathbf{c}\right)$ describes thymic export that depends explicitly on the argument $t$. We appeal to results of Strauss and Yorke in (1967), in particular their Theorem 4.6, which we may invoke to prove that the solution of the perturbed system $\mathbf{c}^{P}(t)\to 0$ if the unperturbed and perturbed systems (U) and (P) satisfy the following conditions:

The zero solution ($\mathbf{c}^{U}(t)\equiv 0$) of the unperturbed system (U) is uniformly asymptotically stable. 2. 2.

The autonomous term $\mathbf{f}(\mathbf{c})$ is $C^{1}$. 3. 3.

There exists $r>0$ such that if $|\mathbf{c}|\leq r$, then $|\mathbf{g}(t,\mathbf{c})|\leq\eta(t)$ for all $t\geq 0$ where $G(t):=\int_{t}^{t+1}\eta(s)ds\to 0$ as $t\to\infty$. (Here, we use the norm $|\mathbf{c}|=\sum_{i=1}^{M}|c_{i}|$.)

We now verify Conditions 1–3 above. Condition 1 follows immediately from the previous discussion, and the fact that for an autonomous system, asymptotic stability and uniform asymptotic stability are equivalent. Condition 2 is trivial. To verify Condition 3, we must construct a suitable function $\eta(t)$, using the definition of the function $g(t,\mathbf{c})$:

[TABLE]

Thus, $|\mathbf{g}(t,\mathbf{c})|\leq\gamma_{0}e^{-at}\left(1+\frac{3}{\Omega}|\mathbf{c}|\right)$, and for a given choice of $r>0$, we may define $\eta_{r}(t):=\gamma_{0}e^{-at}\left(1+\frac{3r}{\Omega}\right)$. From the exponential form of $\eta_{r}(t)$, it is clear that $\lim_{t\to\infty}\int_{t}^{t+1}\eta_{r}(s)ds=0$. Moreover, not only does there exist a single choice of $r>0$ that produces a suitable $\eta_{r}(t)$, but any choice of $r$ produces a suitable $\eta_{r}(t)$.

From Theorem 4.6 in (Strauss and Yorke 1967), we may conclude that for any $T_{0}\geq 0$, there exists a $\delta_{0}>0$ such that if $t_{0}\geq T_{0}$ and $|\mathbf{c}^{P}(t_{0})|\leq\delta_{0}$, then the solution of the perturbed problem, $\mathbf{c}^{P}(t)$, passing through $(t_{0},\mathbf{c}^{P}(t_{0}))$ converges to zero in magnitude as $t\to\infty$. Here the proof of convergence holds for any sufficiently smooth function $\gamma(t)\to 0$. Given Eq. 6 truncated at an arbitrarily large threshold $M$, all $c_{k}$ decline with the decaying thymic export as $t\to\infty$. While the total cell count is preserved by proliferation driving all cells above the truncation threshold and out of the truncated system through truncation errors, the mean-field approximation breaks down at the limit $\gamma(t)/\mu\to 1/\Omega\ll 1$, and Eq. 6 no longer accurately describes the real biology. Nonetheless, our analysis here describes the decline of the number of T-cell clones with decaying $\gamma(t)$ as $t\to\infty$, before the mean-field approximation breaks down.

Appendix D Computation of Expected Sample Clonal Size Distribution

In this section, we detail the derivation of Eq. 10, the explicit expression for $\mathbb{E}[c_{k}^{Y}]$. We begin with Eq. 9,

[TABLE]

Each term $P\left(c_{k}^{Y}=j\right)$ in Eq. 23 can itself be expanded as a sum over all the ways to choose the $j$ clones that are of size $k$. For a sample containing exactly $Z$ clones of size $k$, we introduce the following $Z$-tuple notation, for $Z\in\mathds{N}$:

[TABLE]

where $\vec{i_{Z}}$ lists the indices of all the sample clones consisting of precisely $k$ cells. Additionally, let $y_{i}$ denote the size of the $i$-th ordered sample clone, so that $y_{i_{1}}=y_{i_{2}}=\cdots=y_{i_{Z}}=k$, but no other sample clone consists of $k$ cells. Note that in $\vec{i_{Z}}$, clones are listed in numerical order, due to the assumption $i_{j}<i_{j+1}$, in order to avoid repetition (e.g., in $I_{2}$, $(i_{1},i_{2})$ should be indistinct from $(i_{2},i_{1})$, and this pair should not be counted twice, as the significance is in which clone numbers are listed at all, and not the order in which they are written.) With this, let $P(\vec{i_{Z}},k)$ denote the probability that there are precisely $Z$ clones of size $k$ in the sample, and that their clone numbers are listed in the vector $\vec{i_{Z}}$. Additionally, for $s\in\mathbb{N}$, denote by $I_{Z,s}\subset I_{Z}$ the collection of all $\vec{i_{Z}}\in I_{Z}$ such that $i_{z^{*}}=s$ for some $z^{*}\in\{1,2,\cdots,Z\}$. Essentially, we are imposing the assumption that the $s$-th clone specifically belongs somewhere in the list $\vec{i_{Z,s}}$. Explicitly, we may write $I_{Z,s}$ as:

[TABLE]

We define $P(\vec{i_{Z,s}},k)$ as the probability that there are precisely $Z$ clones of size $k$, with clone numbers listed in $\vec{i_{Z,s}}$, recalling that the $s$-th clone is in the list. We may further simplify Eq. 23 with this notation, rearranging sums by strategically regrouping clone size distributions that share a common size $k$ clone.

[TABLE]

The terms of the final sum in Eq. 29 give the probability that the $s$-th clone is of size $k$, independent of any other information about the sampling. This probability is easy to compute, and given by:

[TABLE]

Inserting Eq. 30 into Eq. 29, we obtain a simple expression for the expected sample clone size distribution:

[TABLE]

We can further simplify Eq. 31 by recognizing that the term $\binom{n_{s}}{k}$ is nonzero only if $n_{s}\geq k$. We can thus rewrite Eq. 31 in terms of the true clone size distribution $\{c_{l}^{N}\}_{l=1}^{R}$ as:

[TABLE]