A mathematical model of the effects of aging on naive T-cell population and diversity
aa r X i v : . [ q - b i o . P E ] J a n Bulletin of Mathematical Biology manuscript No. (will be inserted by the editor)
A mathematical model of the effects of aging on naiveT-cell population and diversity
Stephanie Lewkiewicz · Yao-li Chuang · Tom Chou
Received: date / Accepted: date
Abstract
The human adaptive immune response is known to weaken in advancedage, 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-cellcompartment has been implicated as a likely cause, but the underlying mechanismsdriving this immunosenescence have not been quantitatively modeled and system-atically analyzed. T-cell receptor diversity, or the extent of pathogen-derived anti-gen responsiveness of the T-cell pool, is known to diminish with age, but inherentexperimental difficulties preclude accurate analysis on the full organismal level. Inthis paper, we formulate a mechanistic mathematical model of T-cell populationdynamics on the immunoclonal subpopulation level, which provides quantitativeestimates of diversity. We define different estimates for diversity that depend onthe individual number of cells in a specific immunoclone. We show that diversitydecreases with age primarily due to diminished thymic output of new T-cells andthe resulting overall loss of small immunoclones.
Keywords immunosenescence · T-cell · aging · diversity · thymus Stephanie LewkiewiczE-mail: [email protected] of Mathematics, UCLA, Los Angeles, CA 90095-1555, USAYao-li ChuangE-mail: [email protected] of Mathematics, CalState-Northridge, Northridge, CA 91330-8313, USADepartment of Biomathematics, UCLA, Los Angeles, CA 90095-1766, USATom ChouE-mail: [email protected] of Biomathematics, UCLA, Los Angeles, CA 90095-1766, USADepartment of Mathematics, UCLA, Los Angeles, CA 90095-1555, USA Stephanie Lewkiewicz et al.
Immunosenescence underlies poor health outcomes in the aging population, includ-ing diminished vaccine efficacy (Poland et al. 2010; McElhaney and Dutz 2008;Fleming and Elliot 2008), increased susceptibility to disease (including irregularpresentation, intensified symptoms, longer recovery times, and increased mortal-ity) (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 originatesfrom extensive fundamental changes to the size and functionality of immune cellpools, and the structure of lymphatic tissues in which they develop and oper-ate (Salam et al. 2013).Among the many changes associated with immunosenescence (Globerson and Effros2000), 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 ofwhich retain immunological record of previous infections. The human immune com-partment maintains ∼ T-cells in total, of which ∼ are naive (Jenkins et al.2009; Trepel 1974). During aging, the population of naive T-cells declines in over-all size, while the population of memory T-cells undergoes extensive proliferation,thereby reversing the balance of naive and memory T-cells that had persistedat younger ages (Globerson and Effros 2000; Fagnoni et al. 2000). The expan-sion of memory T-cells further enhances immunological memory of previously-encountered antigens, reinforcing existent immune protection. The remaining naivepool experiences loss of T-cell receptor (TCR) “structural diversity” (Goronzy et al.2007, 2015b)–the number of distinct TCR complexes present across the entirenaive pool. The diversity of T-cell clones, or “immunoclones”, characterized bythe number of distinct TCR complexes among the cell population, provides theextent of antigen specificity. Unique TCR complexes are generated during T-celldevelopment in the thymus, via recombination of genes encoding the V and J do-mains of the TCR α chain and the V, J, and D domains of the TCR β chain,along with additional insertion and deletion of nucleotide fragments (Murphy2012). Combinatorially, a possible Ω ∼ − unique TCR complexesmay be assembled via this rearrangement process (Laydon et al. 2015), but only Ω ∼ (0 . × Ω of those rearrangements are functionally viable (Yates 2014), asdetermined by positive and negative selection tests in the thymus, which screenfor appropriate reactivity to self-peptide/MHC molecules. Each TCR is activatedby at least one peptide fragment presented via MHC molecules on the surfaceof an antigen-presenting cell, thus loss of naive TCR structural diversity lim-its 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, suchas diminished binding affinity and proliferative capacity after antigenic stimula-tion (Moro-Garc´ıa et al. 2013). While these effects have been studied mostly usingmurine models to date (Appay and Sauce 2014), they are not yet well understoodin humans and are beyond the scope of this paper.The total abundance of naive T-cells, which inhabit both blood and lymphatictissue, can be reliably estimated from measurements in small samples (Westermann and Pabst1990; Bains et al. 2009a). Recently, Westera et al. (2015) estimated an ∼
52% de-crease in the naive T-cell population in aging. In contrast, accurate estimationof full-organism TCR structural diversity is currently impeded by experimental itle Suppressed Due to Excessive Length 3 imprecision and the inability to extrapolate small sample data to the full organ-ism (Laydon et al. 2015). Experimentation typically entails DNA sequencing of theTCR α or–more commonly– β chain, in particular the complimentarity-determiningregion 3 (CDR3), which is the site of TCR binding to antigenic peptide and mostsignificant basis for diversity (Murphy 2012).Increasingly sophisticated deep sequencing methods have improved estimatesfor the lower bound on TCR diversity but direct estimation of TCR diversityremains a challenge due to various experimental complications, such as the inabil-ity to detect rare clonotypes, sequencing errors, and inaccurate measurement ofclonotype frequencies resulting from inconsistencies in polymerase chain reaction(PCR) amplification (Laydon et al. 2015). Predicting full-organism TCR diversityfrom a small sample is typically formulated as an “unseen species problem”, andone of many canonical solutions to such a problem is employed in conjunctionwith experimental data (Chao 1984; Chao and Lee 1992; Colwell and Coddington1994), but the true relationship between sample and full diversity is fundamentallyelusive.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 β chain diversity of ∼ × that persisted indonors through age 60, before dropping by two orders of magnitude to ∼ × at age 70. More recently, Britanova et al. (2014) collected samples from donors ofall ages and observed an approximately linear decrease in TCR β CDR3 diversityfrom ∼ × in youth (6 −
25 years) to ∼ . × in advanced age (61 − et al. (2014) obtained a particularly high lower bound estimate of ∼ unique TCR β sequences in youth (20 −
35 years), which declined two- to five-foldin advanced age (70 −
85 years).Note that only the TCR β chain is sequenced in these experiments. Sequencingof both the α and β chains would potentially produce a more accurate measure ofTCR diversity, but the same experimental limitations preclude complete analysis.The measurement of diversity is further complicated by the potentially large dis-parity between structural diversity and “functional diversity”–that is, the numberof antigens to which the T-cell pool is capable of responding. Due to the potentialfor crossreactivity, in which one TCR might respond to many structurally similarpeptide fragments, it is possible that actual TCR diversity is much higher thanstructural diversity indicates. It has been speculated that one TCR might respondto as many as 10 different peptide epitopes (Mason 1998).To obtain lifetime estimates of TCR structural diversity, and develop an in-formed context for discussion of functional diversity, we introduce a mechanisticmathematical model of the generation and replenishment of the lymphocyte poolfrom birth through the end of life. Although experimental assessments of full-system information remain challenging, measurements for the dynamics of eachcomponent related to the T-cell population can be found throughout the litera-ture. Our mathematical approach combines the knowledge of these individual com-ponents to study their interplay, leading to an understanding of the full-systemdynamics. 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 formulationis able to efficiently track the total number of distinct T-cell clones, allowing for afull-system assessment of TCR structural diversity. Stephanie Lewkiewicz et al.
We develop our mathematical model by first constructing the equation governingthe total population size of the naive T-cell pool in Sec. 2.1, through which wequantitatively constrain the primary parameters of our model using experimentalmeasurements found in previous literature. The model that describes the evolutionof immunoclones is derived in Sec. 2.2, allowing us to define and estimate thediversity of the T-cell population in Sec. 2.3. In Sec. 2.4, we inspect the impactof sampling on the estimate of immunoclone diversity, as in practice it is onlypossible to extract a small fraction of the entire T-cell population from a body.2.1 Total T-cell population modelThere are three fundamental immunological mechanisms that sustain the naiveT-cell pool: 1) export of mature naive T-cells from the thymus, 2) peripheralproliferation, and 3) cell removal from the naive pool due to death or phenotypicchanges. These basic mechanisms constitute a birth-death-immigration processdescribed by the ordinary differential equation,d N ( t )d t = γ ( t ) + pN ( t ) − µ ( N ) N ( t ) , (1)where N ( t ) denotes the total T-cell count, γ > p > µ ( N ) > γ ( t ) = γ e − at , where γ > a > et al. (2013), most models incorporate IL-7 regu-lation implicitly in the form of carrying capacity, assuming quick equilibration itle Suppressed Due to Excessive Length 5 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 bothcell proliferation and cell death rates, considering the cytokine’s dual role underlymphoreplete and lymphopenic conditions described above. Our model assumescell-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 deathrate of the form µ ( N ) = µ + µ N N + K , (2)where the first term, µ >
0, is the basal rate of cellular death. The second onedescribes the IL-7-mediated regulation of cell death, with µ > N → ∞ . The quantity K is analogous to a“carrying capacity” and dictates the population at which signalling induced deathstarts to limit the population. The constant rate of cellular proliferation underhealthy conditions is supported by recent studies of Westera et al. (2015), showingnearly identical naive proliferation rates at young and old ages during moderateage-related non-lymphopenic loss of naive cells. IL-7 induced proliferation canarise in unhealthy lymphopenic conditions typically found in severe disease of theimmune system (Brass et al. 2014), cytotoxic drug use (Gergely 1999), radiationtreatment (Grossman et al. 2015), or other abnormal situations. These scenariosare, however, beyond the scope of our analysis.Our model has six adjustable parameters, γ , a , p , µ , µ and K . The firstfour are biologically inherent to the mechanism of T-cell homeostasis, and havebeen measured experimentally in humans and rodents. The last two have to beconstrained via parameter sweeps to match relevant experimental observations.Fig. 1(a) illustrates four qualitatively distinct evolution trajectories of N ( t ) thatmay arise from simulations of the model in the presence of a decaying thymicexport rate γ ( t ) (gray dash-dotted curve). To non-dimensionalize Eqs. 1, 2, weuse a − to rescale t and K to rescale N . The qualitative behavior of our modelis thus controlled by three independent parameters: γ a − K − , ( p − µ ) a − , and µ ( p − µ ) − . The black dashed curve arises when µ ( p − µ ) − <
1. In this casecell proliferation always exceeds cell death, leading to unbounded expansion ofthe naive T-cell population. This scenario is unrealistic, except perhaps duringa period of lymphopenia. For µ ( p − µ ) − ≥
1, cell death is able to balancecell proliferation at a homeostatic carrying capacity N = N ss ( γ = 0), defined by µ ( N ss ( γ = 0)) = p , as γ →
0. As illustrated by the green dotted curve, N ( t ) risesand asymptotically converges towards N ss ( γ = 0) provided that γ a − K − ≪ N ss ( γ = 0) primarily by homeostatic pro-liferation. The model’s behavior makes a transition from proliferation-driven to“thymus-driven” if we increase γ a − K − . As shown by the blue solid curve, N ( t ), driven by increased thymic export, overshoots and approaches N ss ( γ = 0)from above as γ ( t ) → p − µ ) a − ≤
0. In this case cell death always exceeds cell proliferation as γ ( t ) →
0, and N ( t ) → N ss ( γ = 0) = 0. As stated earlier, in this paper we fo-cus on scenarios of healthy aging (lymphoreplete) conditions, which immediatelyrules out the scenarios of unbounded growth (black dashed curve) and complete Stephanie Lewkiewicz et al. collapse of the T-cell population (the red dot-dashed curve), effectively constrain-ing our parameters to physiologically reasonable values µ ( p − µ ) − ≥ p − µ ) a − > years N ( t ) ✕ (cid:1) a - K - (cid:2) ( p - (cid:0) ) -1 thymus-driven unbounded growthproliferation-driventhymus-drivenpopulation collapse K =10 K =10 (a) (b) proliferation-driven K =10 (cid:3) ( N ) N y N o (cid:4) ( N ) = - 52% -1 -1.0-0.8-0.6-0.4 - γ ( t ) ✕ thymic export rate Fig. 1 Qualitative behavior of the total T-cell population model (Eqs. 1,2). (a) The total T-cell population N ( t ) as a function of time (in years) for fourqualitatively distinct scenarios. Unbounded growth arises when µ ( p − µ ) − <
1. and the T-cell population collapses when ( p − µ ) a − <
0. Outside of thesetwo regimes, N ( t ) converges asymptotically to a positive steady state as γ ( t ) →
0. If γ a − K − ≪ N ( t ) is driven primarily by homeostatic proliferation andincreases monotonically towards the constant plateau. Increasing γ a − K − leadsto a transition from proliferation-driven scenario to thymus-driven populations,in which N ( t ) reaches a peak value before converging to the steady state. Thedecaying thymic export rate γ ( t ) is alongside of the N ( t ) curves as a reference. Toquantify the decrease in cell counts with age, we define ¯ N y as the average of N ( t )between ages 20 and 30, and ¯ N o between 70 and 80; then ∆ (cid:0) ¯ N (cid:1) = (cid:0) ¯ N o − ¯ N y (cid:1) / ¯ N y is the relative change in cell counts. The parameter values used are γ = 1 . × , a = 0 . K = 10 and p = 0 . µ = 0 . µ = 0 .
004 for unboundedgrowth, p = 0 . µ = 0 .
18 and µ = 0 .
04 for the collapse scenario, p = 0 . µ = 0 .
17, and µ = 0 . p = 0 . µ =0 .
17, and µ = 0 .
04 for the thymus-driven case. The initial value is N (1) = 10 at t = 1 year. (b) ∆ (cid:0) ¯ N (cid:1) as a function of γ a − K − and µ ( p − µ ) − . When γ a − K − and µ ( p − µ ) − are small, N ( t ) is driven primarily by proliferationand keeps increasing well into old age, leading to positive ∆ ( ¯ N ) values. Conversely,for large γ a − K − and µ ( p − µ ) − , thymic export dominates and N ( t ) peaksat early ages, resulting in negative ∆ ( ¯ N ). The black dotted curve corresponds to ∆ ( ¯ N ) = −
52% as previously reported by Westera et al. for human adults. At fixed µ ( p − µ ) − = 4, we are able to reproduce this curve by setting γ a − K − ≃ K = 10 for our choice of parameter values). The value of ∆ ( ¯ N )increases with decreasing γ a − K − and become positive when γ a − K − . p − µ ) a − = 0 . a = 0 . p has itle Suppressed Due to Excessive Length 7 been measured by Westera et al. (2015) as 0 .
05% day − , or equivalently p = 0 . − . The basal death rate µ can be estimated from the lifespan of T-cells.Based on data from Vrisekoop et al. (2008), De Boer and Perelson (2013) obtainan average naive CD4 + T-cell lifespan of ∼ + lifespan of ∼ . + :CD8 + ratio of 2:1, the averagecombined naive T-cell clearance rate is µ = . year − = 0 .
17 year − . Thymicinvolution with age can be quantified by measuring the decrease in thymic epithe-lial volume (Steinmann 1986), based on which Murray et al. (2003) showed thatthymic output decreases by an average of 4 .
3% per year between ages 0 and 100,implying a decay factor of a = | ln(0 . | ≃ . −
25 years old) at ∼ . × trainedcells daily, or equivalently 5 . × per year (Westera et al. 2015). Assuming thatthis rate is γ ( t ) at t = 25 years, we can back-calculate γ = (5 . × ) × (cid:0) . (cid:1) ≈ . × cell exports / year. Note that these values of p , µ , and a satisfy theconstraint ( p − µ ) a − > µ and K are not available inthe literature, further inspection of Fig. 1(a) reveals that µ and K determinewhether thymic export or homeostatic proliferation dominates the evolution of N ( t ). Through the dimensionless parameters, γ a − K − and µ ( p − µ ) − , thetime at which N ( t ) peaks and how fast it declines from the peak vary with changesto the values of µ and K . Recently, Westera et al. (2015) reported a 52% decreasein total naive T-cell counts between young human adults and elderly individuals,which we can use to quantitatively constrain µ and K . Let us define individu-als 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 unre-lated to age averages out over large sample sizes in clinical data, we may evaluate¯ N y = R N ( t )d t and ¯ N o = R N ( t )d t as the average naive T-cell countsrespectively for the young and the elderly, as illustrated by the shaded areas underthe thymus-domination curve in Fig. 1(a). The relative change in the naive T-cellcount between young and elderly adults can thus be evaluated as ∆ ( ¯ N ) = ( ¯ N o − ¯ N y )¯ N y . (3)Fig. 1(b) plots ∆ ( ¯ N ) as a function of γ a − K − and µ ( p − µ ) − , with a =0 .
044 year − for converting the dimensionless time to years to compute ¯ N y and¯ N o . When γ a − K − . µ ( p − µ ) − . ∆ ( ¯ N ) >
0. Note that the home-ostatic carrying capacity when γ ( t ) = 0 is N ss ( γ = 0) = K (cid:0) µ ( p − µ ) − − (cid:1) − .A small γ a − K − value represents a relatively low thymic export rate, and thecarrying capacity increases rapidly as µ ( p − µ ) − →
1, both of which make itchallenging for thymic output to fill up the T-cell pool to carrying capacity before γ ( t ) considerably decays within t ∼ a − . As a result, N ( t ) does not reach a peakvalue at a young age and continues increasing into old age. The ≈
52% decreasein naive T-cell counts reported by Westera et al. (2015) is depicted by the blackdotted curve. If we set µ ( p − µ ) − = 4, our model can be calibrated to repro-duce this decrease in the cell count by choosing K = 10 ( γ a − K − ≃
41 with γ = 1 . × and a = 0 . K = 10 yields γ a − K − ≃ . ∆ ( ¯ N ) ≃ . K = 10 Stephanie Lewkiewicz et al. results in a moderate decrease in the cell count ( ∆ ( ¯ N ) ≃ − . K = 10 and µ ( p − µ ) − = 4, or equivalently µ = 0 .
04 giventhat p = 0 .
18 and µ = 0 .
17, so that the age-related decline of N ( t ) in our modelis consistent with Westera et al. (2015). (a) ✕
420 0 40 80 120 years N ( t ) N ss -4 -3 -2 -1 p - μ c e ll - popu l a t i on t i m e sc a l e ( | λ - | ) μ = 4×10 -2 (b) thymic expor t (22.7 years) μ = 4×10 -3 μ = 4×10 -1 Fig. 2 Comparison of Thymic Export and Cell-Population EvolutionTime Scales. (a) Plots of N ( t ) and N ss show discrepancy. The γ ( t )-dependencemakes N ss decline monotonically with the exponentially decaying thymic export,and N ss approaches a small positive value as γ ( t ) →
0. The solution N ( t ) evolvestowards N ss but never catches up with it because of a slower evolution time scale.(b) Comparison of timescales of thymic atrophy and cell-population evolution.Thymic atrophy is the faster mechanism for most choices of the system’s param-eters. Increasing µ shortens the time scale of clone evolution, indicating thatthe steady state solution can be a reasonable approximation to the fully time-dependent solution at very large µ and very small p − µ . Here, varying N ss within the range [10 , ] yields almost identical results, and the values of γ and K , chosen within the reasonable parameter regime, do not affect the resultssignificantly. Parameter values used are γ = 1 . × , a = 0 . p = 0 . µ = 0 . K = 10 , Ω = 10 . For (a) µ = 0 .
04, and the initial condition is N (1) = 10 .Note that there exist two intrinsic timescales in Eq. 1; thymic export decays ata rate a , while the homeostatic time scale is controlled by p , µ , and µ . If home-ostasis is much faster than thymic involution, the solution of N ( t ) will quicklyconverge to the quasisteady state solution as γ ( t ) evolves. We compare these twosolutions in Fig. 2(a), where the quasisteady-state solution is obtained by solv-ing for the steady-state solution N ss of Eq. 1 with fixed γ ( t ) at each time t , and N ss ( γ ( t )) (black dashed curve) decreases monotonically with age due to the contin-uous decline of γ ( t ). In contrast, N ( t ) (blue solid curve) slowly rises from the initialconditions N (1) = 10 and does not approach the quasisteady-state level untilage ≈
20 years. The trajectory of N ( t ) then overshoots the declining N ss ( γ ( t )),reaches a peak value, and reverts course to go after N ss ( γ ( t )). However, N ( t ) nevercatches up with N ss ( γ ( t )) before the latter reaches a steady state of very low cellcounts. That N ( t ) keeps lagging behind N ss ( γ ( t )) indicates that the timescale forthe full model solution to converge to the steady state is slower than the evolution itle Suppressed Due to Excessive Length 9 of the nonautonomous term γ ( t ). The results here suggest that steady-state solu-tions cannot adequately describe the temporal evolution of the T-cell populationin the biologically relevant range of parameter values that we have implemented.It is necessary to numerically compute the time-dependent solutions for the fullnonautonomous equation.Indeed, we find a disparity in the rates at which thymic export decays andthe steady state solutions evolve. The latter is provided by the inverse of theeigenvalue of Eq. 1 linearized around N = N ss ( γ ( t )). The eigenvalue takes theform λ = p − ( µ + µ ((3 N K + N ) / (( K + N ) )). Simulations in Fig. 2(b)show that for the biologically relevant parameter values we have implemented, thecell-population evolution timescale, | λ | − (red solid curve), is generally longerthan the timescale of thymic involution ( a − ≃ . a = 0 .
044 as denotedby the horizontal black dotted line). Hence the nonautonomous solutions N ( t ) areexpected to lag behind the thymus-driven steady-state solutions N ss . For N ( t ) tobe reasonably approximated by N ss , the cell population has to evolve much fasterthan thymic involution, corresponding to the regime of very large µ , as indicatedby the blue dash-dotted curve, where cell death is extremely sensitive to the cellpopulation size.2.2 Clonotype Abundance DistributionsQuantification of the populations of individual clonotypes would require analysisof 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 maybe 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,d n i d t = γ ( t ) Ω + pn i − µ ( N ) n i , (4)where γ ( t ) /Ω represents thymic export of naive T-cells of each clonotype (thetotal thymic export rate normalized by the total number of viable TCR combina-tions Ω ), and N ( t ) = P 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 sys-tem (Lythe et al. 2016; Desponds et al. 2015). Indeed, the total numbers of T-cellclonotypes Ω in rodent or human bodies are prohibitively large for direct numer-ical simulations of the full system using Eq. 4. It is thus common to reduce thefull system to a more manageable size with the assumption that the phenomenaunder investigation are scale-invariant. However, it is sometimes difficult to assertwhether a certain property really does not change in a re-scaled system, as non-linear phenomena, such as the Allee effects, often arise in population dynamicsand cast doubt on the scalability of the system. Moreover, some properties, suchas the thymic export rate γ ( t ), are naturally scale dependent. It is not alwaysclear how these quantities should be re-scaled in a reduced system, and they haveusually been omitted by simplification arguments in previous models, which limitsthe applicability of these models. In particular, thymic involution is known to be associated with the age-relatedloss 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 moremanageable full-system model, we consider a formulation that tracks how theexpected number of clones of a given size changes with time. By focusing on clonecount rather than the explicit cell count of each distinct clonotype, we are able toeffectively reduce the number of tracked variables and thus the dimension of themodel. This representation was used by Ewens in population genetics (1972), byGoyal 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 ˆ c k ( t ) to be thenumber of clones represented by exactly k naive T-cells in the organism at time t :ˆ c k ( t ) = Ω X i =1 δ n i ( t ) ,k , (5)where the Kronecker delta function δ x,y = 1 when x = y and 0 otherwise. Bylumping clonotypes of the same cell count into one single variable ˆ c k , this alterna-tive formulation can efficiently describe changes to the TCR clone diversity in thefull system, albeit at the expense of the ability to distinguish each specific clono-type (Morris et al. 2014; Mora and Walczak 2016). Individual clone informationis lost, and n i ( t ) cannot be recovered from ˆ c k ( t ) after the transformation in Eq. 5.Nonetheless, the amount of computation can be significantly reduced by truncat-ing ˆ c k ( t ) at a reasonably large k , as few large clones exist in realistic scenarios, andˆ c k ( t ) for large k is negligible. Letting c ( t ) ≡ h ˆ c ( t ) i denote the expected numberof all possible (thymus-allowed) clonotypes unrepresented in the periphery at time t , and c k ( t ) ≡ h ˆ c k ( t ) i the expected number of clones of size k at time t , a closedset of equations governing the evolution of c k ( t ) can be derived from Eq. 4 in themean-field limit,d c k ( t )d t = γ ( t ) Ω [ c k − − c k ] + p [( k − c k − − kc k ] + µ ( N ) [( k + 1) c k +1 − kc k ] , (6)where N ( t ) = P ∞ i n i ( t ) = P ∞ ℓ =1 ℓc ℓ ( t ). The expected values c k ( t ) are also calledspecies abundances in the ecology literature. The number of unrepresented clonesis c = Ω − P ∞ k =1 c k , and summing Eq. 6 multiplied by k over k = 1 , , · · · recoversEq. 1. The mean-field assumption is articulated in terms such as µ ( P ℓ ℓ ˆ c ℓ )ˆ c k thatinvolve higher-order products of ˆ c k rather than correlations of products of ˆ c k .We have found (unpublished) that this mean-field approximation breaks downonly when γ/µ < /Ω ≪ N ∼ K and all c k ∼ c N . Thus, wereasonably assume that γ ( t ) > µ/Ω allowing the use of the mean-field equations6. In Eq. 6, the terms in the forms of ( γ ( t ) /Ω ) c k , pkc k , and µ ( N ) kc k respectivelyrepresent the effect of thymic export, homeostatic proliferation and cell deathon a T-cell clone already represented by k cells in the peripheral blood. Addingone 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 itle Suppressed Due to Excessive Length 11 clone from the c k -compartment to the c k − -compartment. We approximate theproliferation 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 µ ( N ) kc k , where the IL-7 regulated death rate µ ( N ) is given byEq. 2, reduce c k and increase c k − .For a healthy aging adult, the TCR repertoire is mostly comprised of smallclones 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 ≫
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) weexamine the effect of the truncation clone size M , showing sufficient convergenceof c at t = 40 and 70 to fixed values when M &
30, which indicates that furtherinclusion of clones beyond c has little effects on the solution for t .
70 years. Fornumerical simulations of Eq. 6 in this paper, we set M = 200 to ensure minimaltruncation errors. M
10 40 70 10045678 (a) × c (40) c (70) × (b) c ( t ) × c ( t ) years × c ( t ) Fig. 3 Simulations of Eq. 6. (a) Effect of numerical truncation. We plot c (40)and c (70) as functions of M for 10 ≤ M ≤ M &
30. (b) Temporal evolution of c k ( t ). We plot c ( t ), c ( t ), and c ( t ). Each c k ( t ) curve rises to a peak value and subsequently decreases. As k increases, c k ( t ) decreases in magnitude, and the time at which it reaches the peakvalue is pushed back. Parameter values: γ = 1 . × , a = 0 . p = 0 . µ = 0 . µ = 0 . K = 10 , Ω = 10 . Initial values c (1) = 10 , c (1) = Ω − , c k (1) = 0 for all k ≥ 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 ( t ) peaks around t .
20 years, and thereare many fewer clones of exactly two copies at old ages than at young ages. Incontrast, c ( t ) peaks around age 55, and the numbers of clones that have exactly19 copies are roughly the same between old and young ages, whereas the numberof clones that have exactly 59 copies ( c ( t )) keeps increasing into old ages. The relatively earlier decline of c k ( t ) with smaller k is expected, consideringthat rare clones are introduced into the peripheral circulation primarily by thethymus, which started to involute after birth. With increasing k , the influence ofthymic export on c k ( t ) decreases, whereas the dependence on homeostatic prolif-eration increases. Recalling that the rate of thymic involution is faster than thetime scale for homeostasis to drive the clonal population towards equilibrium, thefast decline of the rare clone population leaves room for larger clones to expand.To accompany the steady state N ss , we compute analogous fixed- γ steady statevalues of the full system, c ss k , in Appendix B. The steady states satisfy c ss k → γ → ≤ k ≤ M . We further show that in spite of the fact that c ss k → N = lim M →∞ P Mk =1 kc ss k > M → ∞ , qualitatively consistent with Eq. 1. Moreover, we prove in Appendix Cthat solutions c k ( t ) of the full nonautonomous system satisfy c k ( t ) → k ≤ M , with arbitrarily large M , as t → ∞ . This result is completely independentof the assumed functional forms of the proliferation and death rates, suggestingthat manipulation of homeostatic regulatory mechanisms cannot prevent the ex-tinction of small T-cell clones caused by decaying γ ( t ). We thus conclude thatthymic involution dictates the age-related decline of the TCR diversity of thenaive compartment.2.3 Diversity of the Naive T-cell RepertoireBy 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 structuraldiversity over a lifetime. Expected TCR structural diversity or “richness” is thetotal number of distinct clones present in the immune compartment, for which wedefine a threshold TCR richness diversity, R q ( t ) = X k ≥ q c k ( t ) , (7)where q ∈ N is a lower threshold, so that the quantity R q ( t ) represents the numberof clones of size at least q present in the immune compartment at time t . Situationsin which such a q -dependent threshold arise may include consideration of immunesurveillance, 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 ( t )is approximately t ∼
16. Higher q lead to older peak ages of R q ( t ), consistent withthe 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 criterionas with total cell counts and compute window-averaged values of R q ( t ) betweenages 20 and 30 for the young and between ages 70 and 80 for the elderly. Bydefining ¯ R y ( q ) ≡ R R q ( t )d t , ¯ R o ( q ) ≡ R R q ( t )d t , we quantify the loss ofrichness by computing its relative change. ∆ ( ¯ R q ) ≡ ( ¯ R o ( q ) − ¯ R y ( q ))¯ R y ( q ) . (8) itle Suppressed Due to Excessive Length 13 years0 40 80 120 R q ( t ) (a) × q = 1 q = 2 q = 3 µ ∆ ( R q ) -1.0-0.50 (b) q ∆ ( R q ) -1.0-0.500.5 Lifetime DecreaseLifetime Increase (c) K = 10 K = 10 µ ∆ ( R ) -0.9-0.7-0.5-0.3 (d) K = 10 K = 10 K = 10 Fig. 4 Simulation of Threshold Richness Diversity. (a) R q ( t ) as a functionof t , for q = 1, 2, 3. R q peaks at later times as q increases. (b) ∆ ( ¯ R q ( t )) for varying q , µ . Higher µ correspond to more severe loss of T-cell clones in advanced age.(c) ∆ ( ¯ R q ) for varying q , K . Small values of q result in a lifetime decrease to R q ,but larger values result in a lifetime increase. This is due to the fact that R q peaks at later times as q increases. (d) ∆ ( ¯ R ) for varying µ , K . Initial values c (1) = Ω − , c (1) = 10 c k (1) = 0 for k ≥
2. Parameter values, whennot varying: Ω = 10 , K = 10 , p = 0 . µ = 0 . µ = 0 . a = 0 . γ = 1 . × .Using the same parameter values as in Fig. 4(a), we plot ∆ ( ¯ R q ) with respectto µ and q in Fig. 4 (b),(c). In Fig. 4(b), ∆ ( ¯ R q ) decreases monotonically withincreasing µ , suggesting that upregulated death rate exacerbates the age-relatedloss of richness, and the impact is more significant for larger q . Fig. 4(c) shows thatwhen K = 10 , ∆ ( ¯ R q ) < q ≤
4. This decreasing trend of R q generally agreeswith the loss of diversity observed in recent experiments where measurements wereavailable across multiple ages (Qi et al. 2014; Britanova et al. 2014). For q = 5 , ∆ ( ¯ R q ) ≈
0, and R q is nearly unchanged between youth and advanced age. For q ≥ ∆ ( ¯ R q ) >
0, indicating higher R q at older ages. Generally, the lifetimedecrease in R q ( t ) occurs with small q , whereas for large q , the trend is reversed, inagreement with our discussion of Fig. 3(b) and Fig. 4(a) regarding peak ages. Thisphenomenon indicates that loss of diversity is primarily due to the extinction ofrare 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 lifetimeincrease to R q ( t ) at higher q .Recent TCR- β sequencing studies have attempted to estimate the change inthe repertoire richness of the naive T-cells with age. Despite the difference inorders 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 ∼ × clonotypes in youth(ages 6 − ∼ . × in aged individuals (ages 61 − etal. (2014), in which a two-to-five-fold decline (i.e., a 50% – 80% drop) betweenyouth (ages 20 −
35) and advanced age (ages 70 −
84) was observed. These resultsare quantitatively consistent with our computation of ∆ ( ¯ R ) for K = 10 – 10 . and 0 . ≤ µ ≤ .
05 in Fig. 4(d), whereas the decline of R q for q ≥ ∆ ( ¯ R )changes between ∼ −
66% and ∼ −
76% for 0 . ≤ µ ≤ .
05 and K = 10 . Incontrast, Fig. 1(b) shows that for the same parameter range, ∆ ( ¯ N ) varies from ∼ −
30% to ∼ − et al. (2005). Increases to the cellular death rate do not cause as much additionalclonal extinction as they do additional cellular extinction, as many surviving clonesare too large to wipe out by the death of a few cells.2.4 Sampling StatisticsConsidering 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 detectedclone sizes in small samples and the true clone sizes in the full organism. Asbefore, denote by N the total number of naive T-cells in the human’s immunecompartment, and Y ≤ N the number of cells collected during sampling fromamong 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 Nk the mean numberof 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 Yk the mean number of clones of size k inthe sampling of Y cells taken from the N total cells. Then the expectation of c Yk ,denoted by E [ c Yk ], is, E [ c Yk ] = R X j =1 jP (cid:16) c Yk = j (cid:17) , (9)where P (cid:0) c Yk = j (cid:1) represents the probability that there are precisely j clones ofsize k in the sampling. Then E [ c Yk ] may be expressed explicitly in terms of the c Nk as: itle Suppressed Due to Excessive Length 15 E [ c Yk ] = R X l = k (cid:0) NY (cid:1) c Nl lk ! N − lY − k ! . (10)(See Appendix D for the detailed proof.) The collection of expressions given byEq. 10 for k = 1 , , · · · , R , yields a linear system of equations solvable for c Nk ,using sampled data for the quantities E [ c Yk ]. More specifically, if we define thevectors b E := ( E [ c Y ] , E [ c Y ] , · · · , E [ c YR ] , ) and E := ( c N , c N , · · · , c NR ), Eq. 10 can bewritten as b E = AE , where A is a constant matrix that has non-zero elements onlyin the upper triangle, with non-zero diagonal entry ( NY ) (cid:0) N − kY − k (cid:1) in position ( k, k ).The equation can always be solved uniquely for E given b E . Thus the full sizedistribution E can be uniquely reconstructed from the expected mean sample sizedistribution b E measured experimentally, provided that the latter can be reliablyestimated through a sufficient number of repeated samplings.In Fig. 5(a), we use Eq. 10 to compute E [ c Yk ] from simulated c Nk , comparingthe predicted sampling results for varying choices of Y . The results indicate thateach decrease by one order of magnitude to the sample size results in a decreaseby roughly the same order of magnitude to the predicted diversity. Thus, diversitypredictions vary with sample size, and small samples do not result in accuratemeasurements of diversity.In Fig. 5(b) we examine how sampling may affect the diagnosis of the age-related TCR richness decline ∆ (cid:0) ¯ R q (cid:1) defined in the previous subsection. We find ∆ (cid:0) ¯ R q (cid:1) , which is negative, increasing with decreasing sampling fraction f , reveal-ing that sampling causes an underestimate of the richness decline. As previouslydiscussed, the decline of TCR richness at old ages is primarily due to the extinc-tion of small clones. Since small clones often evade detection during sampling,their extinction is largely unaccounted for, leading to lessened reduction of therichness measure. When f is very small, most of the small clones have escapeddetection; thus decreasing f further does not change ∆ (cid:0) ¯ R q (cid:1) . Moreover, we notethat ∆ (cid:0) ¯ R (cid:1) , which is the most straightforward measure for age-related loss ofTCR richness, changes from −
73% for the full sample, to −
59% for a samplingfraction f ≤ − , which is close to the value of ∆ (cid:0) ¯ R (cid:1) for the full sample. Thisreaffirms our discussion in the previous subsection that a threshold q > ∆ (cid:0) ¯ R (cid:1) , clones fewer thenthree copies largely evade detection, yielding a result equivalent to ∆ (cid:0) ¯ R (cid:1) of thefull sample, which underestimates the actual decrease of the TCR richness. 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 expo-nentially decaying lifetime thymic export, a constant rate of cellular proliferation,and variable cellular death rate that adjusts to present cell counts and availabilityof survival resources. It depicts the generation of the naive T-cell pool in early lifevia thymic export, and long-term maintenance of the population via peripheral f years0 40 80 12010 (a) True (b) -1 -2 -3 -4 -(cid:5) (cid:6)(cid:7) -1 (cid:8)(cid:9)(cid:10)(cid:11)(cid:12)(cid:13)(cid:14)(cid:15) (cid:16) ( R q ) q = 1 q = 2 q = 3 q = 4 q = 5 R ( t ) f = 10 -1 f = 10 -2 f = 10 -3 Fig. 5 Comparison of Actual and Sampled Richness. (a) True lifetime R ,as well as the expected R that result from extracting 10%, 1%, and 0.1% of thetotal cell count for sampling. ( Y = f × N , with f = 10 − , 10 − , 10 − .) Eachdecrease to the sample size by one order of magnitude results in a decrease to theexpected R by approximately one order of magnitude. (b) The ratio of age-relatedTCR richness decline ∆ (cid:0) ¯ R q (cid:1) as a function of sampling fraction f for clone sizethresholds q = 1 – 5. As f decreases, the value of ∆ (cid:0) ¯ R q (cid:1) increases, indicating alower estimate of the TCR richness decline. When f is very small, ∆ (cid:0) ¯ R q (cid:1) becomesinsensitive to further decreases to f . Parameter values used: γ = 1 . × , a = 0 . p = 0 . µ = 0 . K = 10 , Ω = 10 , µ = 0 .
04. Initial values c (0) = Ω , c k (0) = 0 for k ≥ c k ( t ) (with 1 ≤ k ≤ M ) deplete as t → ∞ , inde-pendent of essentially any restrictive assumptions about the homeostatic prolifer-ative mechanism in the periphery. Concretely, for any choice of proliferation anddeath rates p ( N ) , µ ( N ), that satisfy p (0) , µ (0) > γ ( t ) = γ e − at with γ , a >
0, there exists a sufficiently small δ > c k ( t ) → t → ∞ for all 1 ≤ k ≤ M , provided that P | c k (1) | ≤ δ . Although this resultonly guarantees that trajectories c k ( t ) started sufficiently close to zero convergeto zero, simulation indicates that the basin of attraction to this “zero state” isactually quite large. In fact, for the typical initial conditions used throughout thispaper, simulation suggests convergence of all compartments c k to zero in infinitetime. Although it takes an extremely long time to deplete all c k compartments for1 ≤ k ≤ M , the initial phase of this process can still cause significant loss of T-celldiversity in aging individuals within a human lifespan. Most importantly, we find itle Suppressed Due to Excessive Length 17 that the T-cell loss driven by exponentially-diminishing thymic export alone isrobust against any assumptions about the homeostatic proliferative mechanism inthe periphery, as this outcome is universal for all functional forms of p ( N ) , µ ( N );even a particularly strong homeostatic mechanism (say, one with p (0) ≫ µ (0))cannot rescue a plunging diversity. This, in turn, suggests that in searching fortreatments of age-induced loss of diversity, efforts should be directed at the thy-mus, in particular to maintaining thymic productivity into advanced age.Moreover, we compare the real-time simulations and the quasisteady-state so-lutions of the total cell count, as well as the number of distinct clones, over thecourse of age-related thymic output erosion. We find that our simulation resultskeep lagging behind the quasi steady state solutions, suggesting that the erosiontime scale of thymic output is faster than the time scale for the population dy-namics to relax towards a steady state. Mathematically, this result reveals that theevolution of the T-cell population within the human lifespan is a rather dynamicalphenomenon, which may not be well-described by quasistatic solutions, requiringevaluation of the fully nonautonomous system. Biologically, our results indicatethat the loss of T-cell diversity is a delayed response to thymic involution, andassessment of thymic function may predict the health of the immune system.Although peripheral division cannot salvage the T-cell population on a longtime scale, higher basal proliferation rates may at least delay the erosion of the T-cell compartment, sustaining acceptable effectiveness of the immune system withinthe human lifespan (Naylor et al. 2005). We assumed a constant lifetime rate ofcellular proliferation, but alternative research suggests that proliferation rates mayincrease with age (Naylor et al. 2005). In light of this finding, we briefly inspectthe effect of increased proliferation rates at advanced ages on cellular and clonalloss by modifying p ( N ) and µ ( N ) in Eq. 6. For simplicity, we take the deathrate to be constant ( µ ( N ) = µ > p ( N, t ) = p ( t )(1 − N/K ), where a discrete increase in the proliferation rate is incorporated in p ( t ) = p (1 + rH ( t − T )), with p > H ( t ) the Heaviside function, with T the age at which the rate increases. Theconstant r specifies the increase to the proliferation rate. (Full simulation detailsare given in the caption of Fig. 6.) By varying r , simulation under these alternatehypotheses indicates that increased basal proliferation rates do lead to notablyhigher 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 unlikelyto rescue the eroding naive T-cell diversity, despite the increased cell count. Ifdiversity loss is the main cause of immunosenescence (still a debatable topic inthe medical community), peripheral proliferation may not be the sensible targetof treatments.The increased N (70) and nearly unchanged R (70) in Fig. 6 imply that thedecline of T-cell diversity at old age may appear more dramatic if the diversityis measured in terms of frequency of distinct TCR sequences among the cyclingcells, which corroborates the explanation that an increase of proliferation rate atold age leads to a sharp decrease of T-cell diversity (Naylor et al. 2005). Previousmodels have shown that even sharper decline of T-cell diversity can be induced byfitness selection, where certain clonotypes increase their fitness at old age possiblydue 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 theproliferative mechanism are unable to replenish the declining TCR diversity in thenaive T-cell pool, it is possible that the impact is less severe than the decayingrichness would have indicated, considering that most of the extinct clones areoriginally small clones, which may be much less effective than larger clones. Inthis regard, the viability of treating immunosenescence by expanding peripheralproliferation depends on the elucidation of the T-cell pool’s effectiveness clonesize –that is, the size a clone must have attained to effectively guarantee activationof the clone when its cognate antigen infiltrates the organism. The effectivenessclone size is intrinsically linked to true functional TCR diversity; if we can identifya threshold integer q ∗ , such that clones of size at least q ∗ are reliably activatedin 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 forprecisely those clones actively participating in the adaptive immune mechanism.The larger the “correct” choice of q ∗ is, the more effective treatments to boostcellular proliferation in the periphery will be. Our model directly yields the numberof clones of a particular size, making it straightforward to include or exclude clonesbelow a certain cell count, should such a threshold exist and be identified.The effective clone size is also significant to the question of whether diversityloss is the driving factor in immunosenescence. Using the parameter values that wefound in the literature, R q ( t ) decreases for q ≤ q = 5 ,
6, and increases for q ≥
7. The extinction of small clonesallows the surviving clones to expand in size, leading the richness of large clones toincrease at old ages. If the minimal size for a T-cell clone to effectively respond toantigens is large, the diversity of such “effective” clones may actually increase withage, strengthening the immune response. Therefore, either the minimal clone sizerequired for effective immune response is low, or the weakened immune responseat old ages is caused primarily by other mechanisms. For example, functional defi-ciencies acquired by naive T-cells in aging are one possible alternative cause of theweakened immune response. Such functional deficiencies have been studied heav-ily in mouse models, but research in humans is still lacking (Appay and Sauce2014). Diminished naive T-cell effector responsiveness and proliferative capacityhave been observed in aged mice (Moro-Garc´ıa et al. 2013). It is possible thatsimilar changes occur in humans. Conversely, experiments on mice have directlyshown that loss of TCR diversity does have an actively detrimental effect on im-mune responsiveness (Yagger et al. 2008), supporting the notion that loss of TCRdiversity as a significant contributor to immunosenescence.Our model illustrates the feasibility of several different scenarios, in which lossof diversity contributes to immunosenescence on drastically different levels. Thereis clearly a strong need to investigate the effects of both age-related structuraldiversity loss and T-cell functionality loss in human subjects in vivo , to betterunderstand the causes of immunosenescence. Moreover, our model indicates thatthe effectiveness clone size and crossreactivity in vivo are valuable pieces of missinginformation, the elucidation of which would allow for the identification of effectiveoptions to treat immunosenescence. itle Suppressed Due to Excessive Length 19 proliferation increase, r ∆ ( N ) -1.0-0.500.51.0 (a) T = 30 T = 70 proliferation increase, r ∆ ( R ) -0.8-0.7-0.6-0.5-0.4 (b) T = 30 T = 50 T = 70 Fig. 6 Total Cell Count and Richness with Rise in Proliferation.
Simu-lation of Eq. 6 with exponentially decaying thymic export, and peripheral home-ostasis described by time-varying logistic growth. We use the thymic export rate γ ( t ) = γ e − at , peripheral death rate µ ( N ) = µ >
0, and peripheral proliferationrate p ( N, t ) = p ( t )(1 − ( N/K )), with p ( t ) = p (1 + rH ( t − T )). Here, H ( t ) repre-sents the Heaviside function with jump at t = 0. The constant r determines themagnitude of the increase to the basal proliferation rate, and T represents the timeat which the jump occurs. We take the jump to occur at varying ages. (a) ∆ ( ¯ N )with jump at ages T = 30 and 70, for varying r . (Curve corresponding to T = 50is omitted due to close similarity to T = 30 curve.) Raising the basal proliferationrate diminishes cellular loss in advanced age, with sufficiently high values of r pro-ducing a lifetime increase in total cell counts. The positive steady state solutionof the autonomous total cell ODE, d N/ d t = γ + p (1 − N/K ) − µ N , is given by N ∗ = ( K/ − µ /p + p (1 − µ /p ) + 4 γ /Kp ), and can be seen to satisfy ∂N ∗ /∂p > γ < Kµ , suggesting that increases to the basal proliferationrate are likely to increase the total cell count. (b) ∆ ( ¯ R ) with T = 30 , , and 70,for varying r . Increases to the basal proliferation rate do mitigate diversity loss,but the effect is minor and potentially insignificant. Increases to the basal prolif-eration rate increase c k +1 due to a decrease in c k , preserving additional diversity,but the lifetime diversity loss is still observed, even when proliferation rates arehigh enough to generate a lifetime increase to the total cell count. Fixed param-eter values: γ = 1 . × , a = 0 . p = 0 . µ = 0 . K = 3 × , Ω = 10 . Initial values: c (1) = Ω − , c (1) = 10 c k (1) = 0 for k ≥
1. Eq. 6is truncated at k = 200. We have simulated the time evolution of the functions c k ( t ), which represent thenumber of naive T-cell clones of size k present in a human’s immune compartmentat time t . We determined that under essentially any realistic assumptions abouthomeostatic proliferation and death, all clones deplete in infinite time if thymicexport is assumed to decay exponentially. This implicates thymic export as a fun-damental cause of age-associated diversity loss. We simulated our model under the assumption that a carrying capacity is regulated by homeostatic proliferation anddeath through N -dependent rates. We found that the manipulation of homeostaticproliferation and death rates, which may notably raise the carrying capacity andthus 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 indicatesthat boosting the proliferation rate is unlikely an effective solution. However, ifonly clones of large size are sufficiently effective in the immune response, boostingproliferation rates might raise average clone sizes and help to mitigate the effectsof lost diversity. We simulated “threshold richness diversity”, R q ( t ), which countsthe total number of clones of size q or larger. We found that by increasing q , thetrajectory 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 effec-tive richness would actually increase with age, suggesting that it is important toidentify the minimal effective clone size in order to determine whether the loss ofTCR diversity is the primary driving mechanism of the immune dysfunction seenin advanced age. Lastly, we derived a one-to-one mapping between the full-samplediversity c Nk of N cells and the expected measurement of diversity E [ c Yk ] in sam-ples of Y cells. We found that the probability of detecting small clones shranksignificantly with small sample sizes, which could potentially skew small-samplestatistics. In particular, we show that small samples tend to underestimate theage-related loss of T-cell richness diversity. Our formulation provides a rigorousmethod for accurately inferring the statistical distribution of clonal sizes fromsmall-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).
A Implementation of Numerical Truncation
The most straightforward way to truncate Eq. 6 at k = M is to neglect the exchange termsbetween c M and c M +1 , assuming a negligible contribution for k > M and essentially imposinga “no-flux” boundary condition. This leads to the following equation for the boundary term c M ( t ): d c M ( t )d t = γ ( t ) Ω c M − + p ( M − c M − − µ ( N ) Mc M . (11)This formulation, however, introduces a truncation error in Eq. 1 if we define N = P Mk =1 kc k .The neglected terms leave a small loss of total cell count in d N/ d t . An alternative implemen-tation of the truncation is adding these small loss terms to the boundary equation:d c M ( t )d t = γ ( t ) Ω (cid:16) c M − + c M M (cid:17) + p ( M − c M − + pc M − µ ( N ) Mc M , (12)thus preserving the total cell count N . However, for Eq. 12 the truncation error shows up inthe total number of clonal types Ω = P Mk =0 c k , as the terms added to Eq. 12 to preserve N artificially introduce new clonal types into the model. In contrast, Ω is preserved with theitle Suppressed Due to Excessive Length 21implementation of Eq. 11. If M → ∞ , the truncation errors for both implementations go to zeroat ∼ /M , and the two implementations become equivalent. Assuming sufficiently large M ,the truncation errors can be negligible in the context of γ ( t ) >
0, or have minimal cumulativeeffects within a limited duration, such as a human lifetime, on which our investigations in thispaper have primarily focused.In this paper, we adopt, for simplicity, Eq. 11 to numerically truncate Eq. 6. Note that thischoice 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 carryingcapacity are far more sophisticated than simply eliminating the proliferation of c M . Not onlyshould the proliferation rate of c M go to zero, the proliferation rate of the other c k shouldalso have a k dependence. The k dependence may be weak for small k , but as k → M , theproliferation rate should attenuate significantly. The probability that c k → M will proliferateshould be very small, as it is highly likely that there exist other smaller clones to push thetotal cell count up to the carrying capacity, prohibiting further proliferation. The k -dependentproliferation 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 numericallynegligible and not biologically significant. B Steady States of the Autonomous Equations
If we fix γ ( t ) = γ , Eqs. 1, 6, and 11 become autonomous and admit the following steady statesolution, c ss1 = γ γ Ω M X i =1 i ! µ ( N ss ) i − i − Y j =1 h γ Ω + jp i + µ ( N ss ) − , (13) c ss k = c ss1 k ! µ ( N ss ) k − k − Y n =1 h γ Ω + np i! , (14)where N ss is the total population at steady state, given by the unique positive root of thecubic, c ( N ; γ ) = ( p − ( µ + µ )) N + γ N + ( p − µ ) K N + γ K . (15)When γ = 0, c ( N ; 0) has three real roots, N = 0, ± p (( p − µ ) K ) / ( µ + µ − p ). Thepositive steady state solution, which we denote by N ss (0), is stable, and the zero solutionunstable, under the parameter restrictions described in Section 2.1. We now demonstrate thateven though Eqs. 13, 14 indicate that each c ss k → γ →
0, the quantity lim M →∞ P Mk =1 kc ss k converges to a positive value qualitatively consistent with N ss (0) as γ → Proposition B:
The steady state solutions c ss k , as given in Eqs. 13, 14, satisfy,lim γ → lim M →∞ M X k =1 kc ss k > Proof.
We seek to derive upper and lower bounds, U ( γ ), L ( γ ), which satisfy, L ( γ ) ≤ lim M →∞ M X k =1 kc ss k ≤ U ( γ ) , for γ >
0, and lim γ → U ( γ ) ≥ lim γ → L ( γ ) >
0. We first establish two small results,which will be used later on:2 Stephanie Lewkiewicz et al.
Proposition B1:
For µ = µ ( N ss ( γ )), lim γ → µ d γ > Proof.
Recalling that µ = µ ( N ss ( γ )) = µ + µ ( N ss ( γ ) / ( N ss ( γ ) + K )), we have:d µ d γ = d µ d N ss d N ss d γ = 2 µ K N ss ( N + K ) (cid:20) − ( N + K )3( p − ( µ + µ )) N + 2 γ N ss + ( p − µ ) K (cid:21) = − µ K N ss ( N + K ) [3( p − ( µ + µ )) N + 2 γ N ss + ( p − µ ) K ]where we computed the derivative d N ss d γ implicitly from the expression c ( N ss ( γ ); γ ) = 0.From the explicit form N ss (0) = p ( p − µ ) K / (( µ + µ ) − p ), we have:lim γ −→ d µ d γ = − µ K N SS (0)( N SS (0) + K ) [3( p − ( µ + µ )) N SS (0) + ( p − µ ) K ]= − µ K N SS (0)( N SS (0) + K ) [ − p − µ ) K ] > Proposition B2:
For f ( p/µ ( N ss ( γ )); γ ) = γ pΩ (cid:16) − pµ ( N ss ( γ )) (cid:17) − γ pΩ − , lim γ → f ( p/µ ( N ss ( γ )); γ ) > Proof.
We write the function f ( p/µ ( N ss ( γ )); γ ) as a product of two functions as follows: f ( p/µ ( N ss ( γ )); γ ) = γ pΩ (cid:18) − pµ ( N ss ( γ )) (cid:19) − γ pΩ − = (cid:18) − pµ ( N ss ( γ )) (cid:19) − γ pΩ · γ pΩ (cid:18) − pµ ( N ss ( γ )) (cid:19) − = A ( γ ) · B ( γ )We define A = lim γ → A ( γ ) and B = lim γ → B ( γ ), and compute A and B :itle Suppressed Due to Excessive Length 23ln( A ) = lim γ −→ − γ pΩ ln (cid:18) − pµ ( N ss ( γ )) (cid:19) = − pΩ lim γ −→ ln (cid:16) − pµ ( N ss ( γ )) (cid:17) γ − = − pΩ lim γ −→ (cid:16) − pµ ( N ss ( γ )) (cid:17) − γ (cid:16) − pµ ( N ss ( γ )) (cid:17) − γ − = 1 pΩ lim γ −→ γ (cid:20) − pµ ( N ss ( γ )) (cid:21) − (cid:20) pµ ( N ss ( γ )) − d µ d γ (cid:21) = 1 pΩ lim γ −→ γ p d µ d γ µ ( N ss ( γ )) − pµ ( N ss ( γ )) = 1 Ω lim γ −→ γ µ d γ + γ
20 d µ d γ (2 µ − p ) d µ d γ = 1 Ω γ lim γ → µ d γ + γ lim γ → µ d γ p lim γ → µ d γ , where we used that µ ( N ss ( γ )) → p as γ →
0. From Proposition B1, lim γ −→ µ d γ >
0, anda similar computation shows that lim γ → µ d γ ∈ R . Thus, ln( A ) ∈ R , and A >
0. Now, B = lim γ → γ pΩ (cid:18) − pµ ( N ss ( γ )) (cid:19) − = lim γ → ( γ /pΩ ) (cid:16) − pµ ( N ss ( γ )) (cid:17) = lim γ → (1 /pΩ ) pµ ( N ss ( γ )) − µ d γ = lim γ → µ ( N ss ( γ )) p Ω d µ d γ > . Thus, lim γ → γ pΩ (cid:16) − pµ ( N ss ) (cid:17) − γ pΩ − = A B > c ss1 , to simplify calculations. From the nonnegativity of the parameters and coefficientfunctions, and the form in Eq. 13, c ss1 ≤ γ /µ , independent of M . To derive an M -independentlower bound on c ss1 , we observe that the sum in the denominator of Eq. 13 satisfies, γ Ω M X i =1 i ! µ ( N ss ( γ )) i − i − Y j =1 h γ Ω + jp i ≤ M X i =1 i − µ ( N ss ( γ )) i − i − Y j =0 h γ Ω + jp i = p M X i =1 i − i − Y j =0 (cid:20) γ pΩ + j (cid:21) (cid:18) pµ ( N ss ( γ )) (cid:19) i − M -th Taylor polynomial, S M,γ , for the function f ( x ; γ ) = γ pΩ (1 − x ) − γ pΩ − expanded around x = 0 and evaluated at x = pµ ( N ss ( γ )) . Thefunction f ( x ; γ ) is analytic in x away from x = 1, and in particular, the S M,γ increasemonotonically to f ( p/µ ( N ss ( γ )); γ ). It follows that,1 p M X i =1 i ! µ ( N ss ( γ )) i − i − Y j =0 h γ Ω + jp i ≤ S M,γ ≤ f (cid:18) pµ ( N ss ( γ )) ; γ (cid:19) := f γ and thus that c ss1 ≥ γ / ( pf γ + µ + µ ). After using the c ss1 bounds in the expression for c ss k ,we have: pΩpf γ + µ + µ S M,γ ≤ M X k =1 kc ss k ≤ pΩµ S M,γ −→ lim M →∞ pΩpf γ + µ + µ S M,γ ≤ lim M →∞ M X k =1 kc ss k ≤ lim M →∞ pΩµ S M,γ −→ pΩpf γ + µ + µ f γ ≤ lim M →∞ M X k =1 kc ss k ≤ pΩµ f γ Now we let L ( γ ) = pΩpf γ + µ + µ f γ and U ( γ ) = pΩµ f γ . From Proposition B2, lim γ → f γ >
0, so lim γ → L ( γ ) , lim γ → U ( γ ) >
0, and Proposition B follows.
C Convergence and Stability of c k when γ ( t ) → In this section we will prove that solutions c k to our ODE system initialized sufficiently close to −→ −→ t → ∞ . Denote by (P) the “perturbed” ODE system given by Eqs. 6, 11,with γ ( t ) = γ e − at , and by (U) the “unperturbed” ODE system resulting from the alternatechoice γ ( t ) ≡
0. For the sake of generality, we omit previous assumptions about the form ofthe functions p ( N ) , µ ( N ), except that p (0) , µ (0) >
0. Additionally, in this section, we regardthe term N that appears in the ODEs as P k ≥ kc k instead of its own variable, and thusdo not explicitly include Eq. 1 in our analysis as in Appendix B. Note that the residual N − P k ≥ kc k → M → ∞ . We begin by noting that the unperturbed system (U) hassteady-state c Uk ( t ) ≡ k ≥
1. To analyze the stability of this steady state, we consider thelinearization of (U) around this steady state, which is represented by the M × M matrix wecall L U ( L U = ( l ij ) ≤ i,j ≤ M ). The components l ij of L U are given explicitly by: l ij = − j ( p (0) + µ (0)) , if i = j ≤ M − − Mµ (0) , if i = j = Mjµ (0) , if i = j −
1; 2 ≤ j ≤ Mjp (0) , if i = j + 1; 1 ≤ j ≤ M − , otherwise (16)Although the matrix is tridiagonal, it is high-dimensional, and thus its eigenvalues cannot becomputed analytically. However, we may nevertheless demonstrate that all eigenvalues possessstrictly 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 λ ∈ C satisfying ℜ ( λ ) ≥
0, then λ = 0. We then verify that λ = 0 is never an eigenvalue of L U , by directlydemonstrating that L U has linearly independent rows.itle Suppressed Due to Excessive Length 25 Proposition C:
All eigenvalues λ ∈ C of the matrix L U satisfy ℜ λ < , so that the zero-solution of (U) is asymptotically stable. We first apply Gershgorin’s circle theorem to the columns of the matrix L U to concludethat all eigenvalues λ ∈ C of the truncated system (finite M ) are contained within the followingunion of disks: M − [ i =1 { λ ∈ C : | λ + i ( p (0) + µ (0)) | ≤ i ( p (0) + µ (0)) } ! [ { λ ∈ C : | λ + Mµ (0) | ≤ Mµ (0) } , (17)where we have used the fact that { λ ∈ C : | λ + D | ≤ D } ⊂ { λ ∈ C : | λ + ( D + ǫ ) | ≤ D + ǫ } for D, ǫ >
0. Given the assumption that p (0) , µ (0) >
0, each of these disks is tangent to theline ℜ λ = 0 at λ = 0, and otherwise lies entirely in the half plane ℜ λ <
0. Thus, L U can onlypossess an eigenvalue λ satisfying ℜ λ = 0 if λ = 0 is itself an eigenvalue. We next verify that λ = 0 is never an eigenvalue of L U directly, by establishing the linear independence of therows of L U .Let us assume that there exist scalars a , a , . . . , a M , such that P Mj =1 a j ( l ij −
0) = 0 forall 1 ≤ i ≤ M . Hence a normalized vector a = ( a , a , . . . , a M ) represents the eigenvectorof the zero eigenvalue. For i = 1, we find that 2 a µ (0) − a ( p (0) + µ (0)) = 0, so that a =2 − µ (0) − ( p (0) + µ (0)) a . By moving on to larger i , we can recursively derive a i = Θ i a for all 2 ≤ i ≤ M with a proportional constant coefficient Θ i . Moreover, P Mi =1 P Mj =1 a j l ij = − a µ (0) = 0, leading to a = 0 given that µ (0) >
0. If a = 0, a ≡
0, and a non-zero eigenvectordoes not exist, implying that zero is not among the eigenvalues of the M × M matrix L U . Wethus conclude that all eigenvalues λ of the matrix L U satisfy ℜ ( λ ) <
0, and the zero-solutionof (U) is asymptotically stable for Eq. 6 truncated using Eq. 11 at an arbitrarily large M . Notethat the proof in Eq. 17 does not hold if we use the alternative truncation formula Eq. 12. Byforcing all cells to remain below the truncation threshold M , it is not possible for all c k togo to zero with a finite M . For the alternative truncation, the stable steady state solution is c k = 2 N ss / ( M ( M + 1)), which nevertheless goes to zero as M → ∞ .We next proceed to demonstrate that the uniform asymptotic stability of the zero-solution( c Uk ( t ) ≡ k ≥
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 γ ( t ), implies that solutions ofthe perturbed system (P) also converge to zero in magnitude, in a sense to be made moreprecise later on. Here let us simplify our notation by writing (U) as d c / d t = f ( c ), where c ≡ ( c , c , . . . , c M ). The autonomous term f ( c ) consists of cell proliferation and death. Cor-respondingly we express (P) as d c / d t = f ( c ) + g ( t, c ), where the nonautonomous term g ( t, c )describes thymic export that depends explicitly on the argument t . We appeal to results ofStrauss and Yorke in (1967), in particular their Theorem 4.6, which we may invoke to provethat the solution of the perturbed system c P ( t ) → c U ( t ) ≡
0) of the unperturbed system (U) is uniformly asymptoticallystable.2. The autonomous term f ( c ) is C .3. There exists r > | c | ≤ r , then | g ( t, c ) | ≤ η ( t ) for all t ≥ G ( t ) := R t +1 t η ( s ) ds → t → ∞ . (Here, we use the norm | c | = P Mi =1 | c i | .)We now verify Conditions 1–3 above. Condition 1 follows immediately from the previous dis-cussion, and the fact that for an autonomous system, asymptotic stability and uniform asymp-totic stability are equivalent. Condition 2 is trivial. To verify Condition 3, we must constructa suitable function η ( t ), using the definition of the function g ( t, c ):6 Stephanie Lewkiewicz et al. | g ( t, c ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) γ e − at Ω Ω − M X j =1 c j − c (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + M − X j =2 (cid:12)(cid:12)(cid:12)(cid:12) γ e − at Ω ( c j − c j +1 ) (cid:12)(cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12)(cid:12) γ e − at Ω c M − (cid:12)(cid:12)(cid:12)(cid:12) (18) ≤ γ e − at Ω | Ω | + M X i =1 | c i | ! + | c | ! + M − X j =2 γ e − at Ω ( | c j | + | c j +1 | ) + γ e − at Ω | c M − | (19) ≤ γ e − at Ω Ω + 3 M − X i =1 | c i | ! (20) ≤ γ e − at Ω ( Ω + 3 | c | ) (21)= γ e − at (cid:18) Ω | c | (cid:19) (22)Thus, | g ( t, c ) | ≤ γ e − at (cid:0) Ω | c | (cid:1) , and for a given choice of r >
0, we may define η r ( t ) := γ e − at (cid:0) rΩ (cid:1) . From the exponential form of η r ( t ), it is clear that lim t →∞ R t +1 t η r ( s ) ds = 0.Moreover, not only does there exist a single choice of r > η r ( t ), but any choice of r produces a suitable η r ( t ).From Theorem 4.6 in (Strauss and Yorke 1967), we may conclude that for any T ≥ δ > t ≥ T and | c P ( t ) | ≤ δ , then the solution of the perturbedproblem, c P ( t ), passing through ( t , c P ( t )) converges to zero in magnitude as t → ∞ . Herethe proof of convergence holds for any sufficiently smooth function γ ( t ) →
0. Given Eq. 6truncated at an arbitrarily large threshold M , all c k decline with the decaying thymic exportas t → ∞ . While the total cell count is preserved by proliferation driving all cells above thetruncation threshold and out of the truncated system through truncation errors, the mean-fieldapproximation breaks down at the limit γ ( t ) /µ → /Ω ≪
1, and Eq. 6 no longer accuratelydescribes the real biology. Nonetheless, our analysis here describes the decline of the number ofT-cell clones with decaying γ ( t ) as t → ∞ , before the mean-field approximation breaks down. D Computation of Expected Sample Clonal Size Distribution
In this section, we detail the derivation of Eq. 10, the explicit expression for E [ c Yk ]. We beginwith Eq. 9, E [ c k ] = R X j =1 jP (cid:16) c Yk = j (cid:17) . (23)Each term P (cid:0) c Yk = j (cid:1) in Eq. 23 can itself be expanded as a sum over all the ways to choosethe j clones that are of size k . For a sample containing exactly Z clones of size k , we introducethe following Z -tuple notation, for Z ∈ N : I Z := { i Z = ( i , i , . . . , i Z ) : i j ∈ { , , . . . , R } , i j < i j +1 for all j } . (24)where 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 = y i = · · · = y i Z = k , but noother sample clone consists of k cells. Note that in i Z , clones are listed in numerical order, dueto the assumption i j < i j +1 , in order to avoid repetition (e.g., in I , ( i , i ) should be indistinctfrom ( i , i ), and this pair should not be counted twice, as the significance is in which clonenumbers are listed at all, and not the order in which they are written.) With this, let P ( i Z , k )denote the probability that there are precisely Z clones of size k in the sample, and that theiritle Suppressed Due to Excessive Length 27clone numbers are listed in the vector i Z . Additionally, for s ∈ N , denote by I Z,s ⊂ I Z thecollection of all i Z ∈ I Z such that i z ∗ = s for some z ∗ ∈ { , , · · · , Z } . Essentially, we areimposing the assumption that the s -th clone specifically belongs somewhere in the list i Z , s .Explicitly, we may write I Z,s as: I Z,s = { i Z , s = ( i , . . . , i z ∗ − , i z ∗ = s, i z ∗ +1 , . . . , i Z ) : i j ∈ { , , . . . , R } , i j < i j +1 for all j } . (25)We define P ( i Z , s , k ) as the probability that there are precisely Z clones of size k , with clonenumbers listed in i Z , s , recalling that the s -th clone is in the list. We may further simplify Eq. 23with this notation, rearranging sums by strategically regrouping clone size distributions thatshare a common size k clone. E [ c k ] = R X j =1 jP ( c Yk = j ) , (26)= R X j =1 j X i j ∈ I j P ( i j , k ) , (27)= R X s =1 R X j =1 X i j , s ∈ I j,s P ( i j , s , k ) , (28)= R X s =1 P ( y s = k ) , (29)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: P ( y s = k ) = 1 (cid:0) NY (cid:1) (cid:16) n s k (cid:17)(cid:16) N − n s Y − k (cid:17) . (30)Inserting Eq. 30 into Eq. 29, we obtain a simple expression for the expected sample clone sizedistribution: E [ c k ] = R X s =1 (cid:0) NY (cid:1) (cid:16) n s k (cid:17)(cid:16) N − n s Y − k (cid:17) . (31)We can further simplify Eq. 31 by recognizing that the term (cid:0) n s k (cid:1) is nonzero only if n s ≥ k .We can thus rewrite Eq. 31 in terms of the true clone size distribution { c Nl } Rl =1 as: E [ c k ] = R X l = k (cid:0) NY (cid:1) c Nl (cid:16) lk (cid:17)(cid:16) N − lY − k (cid:17) . (32) References
Appay V, Sauce D (2014) Naive T cells: The crux of cellular immune aging? ExperimentalGerontology 54:90–93Bains I, Antia R, Callard R, Yates AJ (2009a) Quantifying the development of the peripheralnaive CD4+ T-cell pool in humans. Immunobiology 113(22):5480–54878 Stephanie Lewkiewicz et al.Bains I, Thi´ebaut R, Yates AJ, Callard R (2009b) Quantifying thymic export: Combiningmodels of naive T cell proliferation and TCR excision circle dynamics gives an explicitmeasure of thymic output. The Journal of Immunology 183(7):4329–4336Berzins SP, Boyd R, Miller JF (1998) The role of the thymus and recent thymic migrantsin the maintenance of the adult peripheral lymphocyte pool. The Journal of ExperimentalMedicine 187(11):1839–1848Bradley LM, Haynes L, Swain SL (2005) IL-7: maintaining T-cell memory and achieving home-ostasis. Trends in Immunology 26(3):172–176Brass D, McKay P, Scott F (2014) Investigating an incidental finding of lymphopenia. BritishMedical Journal 348:1–3Britanova OV, Putintseva EV, Shugay M, Merzlyak EM, Turchaninova MA, Staroverov DB,Bolotin DA, Lukyanov S, Bogdanova EA, Mamedov IZ, Lebedev YB, Chudakov DM (2014)Age-related decrease in TCR repertoire diversity measured with deep and normalized se-quence profiling. The Journal of Immunology 192(6):2689–2698Chao A (1984) Nonparametric estimation of the number of classes in a population. Scandana-vian Journal of Statistics 11(4):265–270Chao A, Lee SM (1992) Estimating the number of classes via sample coverage. Journal of theAmerican Statistical Association 87(417):210–217Colwell RK, Coddington JA (1994) Estimating terrestrial biodiversity through extrapolation.Philosophical Transactions of the Royal Society B 345(1311):101–118de Boer RJ, Perelson AS (2013) Quantifying T lymphocyte turnover. Journal of TheoreticalBiology 327:45–87Desponds J, Mora T, Walczak A (2015) Fluctuating fitness shapes the clone-size distributionof immune repertoires. Proceedings of the National Academy of Sciences 113(2):274–279Desponds J, Mayer A, Mora T, Walczak AM (2017) Population dynamics of immune reper-toires. ArXiv e-prints
Ewens W (1972) The sampling theory of selectively neutral alleles. Theoretical PopulationBiology 3(1):87–112Fagnoni FF, Vescovini R, Passeri G, Bologna G, Pedrazzoni M, Lavagetto G, Casti A,Franceschi C, Passeri M, Sansoni P (2000) Shortage of circulating naive CD8+ T cellsprovides new insights on immunodeficiency in aging. Blood 95(9):2860–2868Fleming DM, Elliot AJ (2008) The impact of influenza on health and health care utilisationof elderly people. Vaccine 32(1):S1–S9Fry TJ, Mackall CL (2005) The many faces of IL-7: from lymphopoesis to peripheral T cellmaintenance. The Journal of Immunology 174(11):6571–6576Gergely P (1999) Drug-Induced Lymphopenia. Drug Safety 21(2):91–100Ginaldi L, Loreto MF, Corsi MP, Modesti M, de Martinis M (2001) Immunosenescence andinfectious diseases. Microbes and Infection 3(10):851–857Globerson A, Effros RB (2000) Aging of lymphocytes and lymphocytes in the aged. Immunol-ogy Today 21(10):515–521Goronzy JJ, Lee WW, Weyland CM (2007) Aging and T-cell diversity. Experimental Geron-tology 42(5):400–406Goronzy JJ, Fang F, Cavanagh MM, Qi Q, Weyand CM (2015a) Na¨ıve T cell maintenanceand function in human aging. Journal of Immunology 194(9):4073–4080Goronzy JJ, Qi Q, Olshen RA, Weyland CM (2015b) High-throughput sequencing insightsinto T-cell receptor diversity in aging. Genome Medicine 7:1–3Goyal S, Kim S, Chen ISY, Chou T (2015) Mechanisms of blood homeostasis: lineage trackingand a neutral model of cell populations in rhesus macaques. BMC Biology 13(85):1–14Grossman SA, Ellsworth S, Campian J, Wild AT, Herman JM, Laheru D, Brock M, Bal-manoukian A, Ye X (2015) Survival in patients with severe lymphopenia following treatmentwith radiation and chemotherapy for newly diagnosed solid tumors. Journal of the NationalComprehensive Cancer Network 13(10):1225–1231Gruver A, Hudson L, Sempowski J (2007) Immunosenescence of aging. The Journal of Pathol-ogy 211(2):144–156Hapuarachchi T, Lewis J, Callard RE (2013) A mechanistic mathematical model for naiveCD4 T cell homeostasis in healthy adults and children. Frontiers in Immunology 4(366):1–6Jenkins MK, Chu HH, McLachlan JB, Moon JJ (2009) On the composition of the preimmunerepertoire of T cells specific for peptide-major histocompatibility ligands. Annual Reviewof Immunology 28:275–294itle Suppressed Due to Excessive Length 29Johnson PLF, Yates AJ, Goronzy JJ, Antia R (2012) Peripheral selection rather than thymicinvolution explains sudden contraction in naive CD4 T-cell diversity with age. PNAS109(52):21432–21437Johnson PLF, Goronzy JJ, Atia R (2014) A population biological approach to understandingthe maintenance and loss of the T-cell repertoire during aging. Immunology 142(2):167–175Laydon DJ, Bangham CRM, Asquith B (2015) Estimating T-cell repertoire diversity: limita-tions of classical estimators and a new approach. Philosophical Transactions of the RoyalSociety B 370(1675):1–11Lythe G, Callard RE, Hoare RL, Molina-Par’ıs C (2016) How many TCR clonotypes does abody maintain? Journal of Theoretical Biology 389:214–224Mason D (1998) A very high level of crossreactivity is an essential feature of the T-cell receptor.Trends in Immunology 19(9):395–404McElhaney JA, Dutz JP (2008) Better influenza vaccines for older people: What will it take?The Journal of Infectious Diseases 198(5):632–634Mehr R, Perelson AS, Fridkis-Hareli M, Globerson A (1996) Feedback regulation of T celldevelopment: manifestations in aging. Mechanisms of Ageing and Development 91(3):195–210Mehr R, Perelson AS, Fridkis-Harelic M, Globersond A (1997) Regulatory feedback pathwaysin the thymus. Immunology Today 18(12):581–585Metcalf D (1963) The autonomous behavior of normal thymus grafts. Australian Journal ofExperimental Biology and Medical Sciences 41:437–444Mora T, Walczak A (2016) Quantifying lymphocyte receptor diversity. ArXiv e-prints1604.00487