Dynamics of T cell receptor distributions following acute thymic atrophy and resumption

Naive human T cells are produced and developed in the thymus, which atrophies abruptly and severely in response to physical or psychological stress. To understand how an instance of stress affects the size and “diversity” of the peripheral naive T cell pool, we derive a mean-field autonomous ODE model of T cell replenishment that allows us to track the clone abundance distribution (the mean number of different TCRs each represented by a specific number of cells). We identify equilibrium solutions that arise at different rates of T cell production, and derive analytic approximations to the dominant eigenvalues and eigenvectors of the mathematical model linearized about these equilibria. From the forms of the eigenvalues and eigenvectors, we estimate rates at which counts of clones of different sizes converge to and depart from equilibrium values—that is, how the number of clones of different sizes ”adjusts” to the changing rate of T cell production. Under most physiological realizations of our model, the dominant eigenvalue (representing the slowest dynamics of the clone abundance distribution) scales as a power law in the thymic output for low output levels, but saturates at higher T cell production rates. Our analysis provides a framework for quantitatively understanding how the clone abundance distribution evolves under small changes in the overall T cell production rate.


Introduction
The thymus, a small organ located above the heart in humans, is a crucial component of the primary lymphoid architecture, as the site of T cell development [1,2]. The many different T cell subpopulations together guide and assist the action of other immune agents during infection [3], regulate the immune response [4], and retain memory of encountered pathogens [5]. As such, the thymus supplies the immune compartment with its most essential source of direction, support, and regulation. T cells are produced when lymphocyte progenitors derived from hematopoietic stem cells in the bone marrow migrate to the thymus and begin a process of role selection, maturation, and vetting, before being exported to the peripheral blood [6]. The most significant event during thymocyte development is the rearrangement of the α and β chains of the T cell receptor (TCR) [7] that occurs in the thymic cortex, which is populated with thymocyte progenitors. The particular rearrangement a T cell undergoes determines its antigen specificity; a naive T cell in the peripheral pool is activated when its TCR is bound by a cognate antigen, a pathogen-derived peptide fragment capable of stimulating that particular TCR [8]. The total number of distinct TCRs present across the full T cell pool is the "TCR diversity" [9], and this quantifies the breadth of the pool's antigen responsiveness [10]. Thymocytes also undergo negative selection to eliminate cells that react too strongly to self antigens presented by resident macrophages and dendritic cells. The small number of T cells that survive this process are functionally competent and thus exported to the peripheral blood to participate in the immune mechanism.
The thymus is known to experience both chronic and acute forms of atrophy [11], resulting from both normal biological processes and the presence of disease or stress. The most universal form of thymic atrophy is age-related involution, the process by which productive thymic tissue is gradually replaced with nonproductive fat [12]. Involution begins at puberty and continues indefinitely, and the resulting decline in T cell production has been implicated as a likely source of immune dysfunction in the elderly [13][14][15]. Acute atrophy can occur under a plethora of conditions associated with a state of disease or stress [16][17][18], including viral, bacterial and fungal infection [19][20][21], malnutrition [22], cancer treatment [23], bone marrow transplant [24], psychological stress and pregnancy [25][26][27]. Each condition facilitates thymic atrophy in (at least) one of several ways, including reducing thymic cellularity [11], decreasing thymocyte proliferation and increasing apoptosis [28], instigating premature export of underdeveloped thymocytes [29], and inducing morphological changes to the thymic microenvironment [30]. Such disturbances may consequently alter the size and composition of the peripheral T cell pool. Decreased lymphocyte prevalence in the periphery during acute involution has been documented [31][32][33][34], and Salmonella, which infects the thymus itself, has been shown to disrupt positive and negative selection, producing a skewed TCR repertoire [35]. Radiation and chemotherapy drugs, such as temozolmide, used to treat cancer can also be highly lymphotoxic, producing a lymphopenic state referred to as "treatment-related lymphopenia" (TRL) [36][37][38]. Viral infections, particularly HIV, and autoimmune disorders can induce lymphopenia by increasing peripheral cellular death and redistributing cells to inappropriate tissues, in addition to affecting production in the thymus [39]. Congenital thymic aplasia, as seen in complete DiGeorge syndrome, results in a lymphopenic state at birth [40].
The activation of the hypothalamic-pituitary-adrenal axis by stress stimuli and subsequent release of glucocorticoids, which are known to induce apoptosis in double-positive thymocytes [41] and inhibit their differentiation [28], is also likely a major underlying catalyst of this acute involution [17,42]. Evidence suggests that glucocorticoid release is actually necessary to affect thymic atrophy [18,28]. Several other chemical agents have been observed to participate in thymic atrophy, notably sex hormones [17], which have been shown to weaken thymocyte proliferation [25] and induce apoptosis [43], and the IL-6 cytokine family, which is demonstrably thymosuppressive [16]. Despite this apparent sensitivity to stress, the thymus is highly plastic, and generally recovers in size and functionality after removal of the stressor [11]. Studies of the thymus during and after chemotherapy treatment in cancer patients indicate a return of thymic volume and productivity during recovery from treatment [23]. A recuperating thymus may even surpass its pre-treatment volume, in a phenomenon known as "thymic rebound" [44,45]. Such thymic recovery has also been seen after infection [46] and traumatic injury [47]. However, recovery is demonstrably age-dependent, with thymi of older patients reconstituting the naive T cell compartment more weakly than those of younger patients [23]. Although acute thymic atrophy has been observed extensively in humans, much has yet to be learned about it, and clear treatment protocol is lacking [11].
To this end, we present a mechanistic mathematical model to predict changes in the size and diversity of the peripheral naive T cell compartment in response to various immunologically diseased conditions. We study how this pool's size and composition adjust to changes in the rate of thymic output. We compartmentalize the peripheral T cell pool by grouping clones-collections of T cells with the same TCR-according to their size. We then use a high-dimensional autonomous ODE system to follow the time evolution of the number of clones in each compartment. We assume that the size of the peripheral naive T cell pool is dictated by rates of thymic export of new T cells, along with homeostatic proliferation and death mechanisms. We assume a piecewise constant rate of thymic export, as the atrophy/recovery cycle is known to be a rapid process, and that the proliferative and death processes are subject to homeostatic regulation based on the total T cell pool size. We derive analytic approximations to the dominant eigenvalues and eigenvectors of the system linearized around its equilibria in both the presence and the absence of thymic activity. From this, we assess the rates of convergence of different T cell compartments to equilibria that result from a changing thymic export rate. We then compare the linearized and fully nonlinear models, and study several special cases. We also compute explicit representations of solutions to an infinite-dimensional extension of our model.

Mathematical model and analysis
We assume that the total naive T cell population N(t) in the immune compartment (the blood and lymphatic tissue) satisfies a general ODE of the form where γ ≥ 0 is the rate of naive T cell export from the thymus, and p(N), μ(N) ≥ 0 are regulated, N-dependent rates of proliferation and death of naive T cells in the peripheral bloodstream. To prevent unbounded growth, we take p(N) to be non-increasing and μ(N) to be non-decreasing as cell counts, N, increase. We assume that p(0) > μ(0), as the lymphopenic proliferation rate should be higher than the lymphopenic death rate [48][49][50][51]. At steady-state, when a healthy, homeostatic cell count N* is achieved, increasing functions μ(N)) admit multiple-typically two-fixed points. The N = 0 fixed point is unstable, while the one at N > 0 is stable.
In order to compute the peripheral naive T cell diversity, we couple Eq 1 with a system of ODEs that describes the time evolution of the size-segregated subpopulations of the peripheral naive T cell pool. Let c k (t) denote the number of clones that are of size k at time t ≥ 0. As formally shown in Song & Chou [52], the mean clone count is c k ∝ P (n i = k, t), the marginalized probability that any single clone i has population k. The master equation for P(n i , t) is difficult to solve with regulation terms. Here, we provide a heuristic derivation of the equations obeyed by c k (t) by using a mean-field approximation that the total population N is uncorrelated with any n i , (although we know N = ∑ i = 1 ∞ n i ). Under this mean-field approximation, the evolution of P(n i = k, t), and hence c k (t), obeys where k = 2, 3,..., M − 1, p(N) and μ(N) are approximated by p(N(t)) and μ(N(t)) (where N(t) is the solution to Eq 1), and the index M in Eq 4 is the hypothetical maximum size a clone can achieve. We take M to be finite for mathematical tractability and in accordance with evidence of intraclonal competition that restricts clone sizes and preserves a balanced TCR diversity [53]. Each of the ODEs in Eqs 2-4 describes how c k changes due to the effects of thymic export of new cells, and proliferation and death in the periphery. The constant Ω denotes the large total number of clonotypes that can potentially be assembled in and exported from the thymus. The basic model includes immigration, birth, and death of multiple species (i.e., clones) and can be developed in a fully stochastic setting; however, in that case, only steady-state solutions are available [54].
In Eq 3, the term γ Ω represents the rate at which cells of a given clonotype are exported to the periphery from the thymus, and thus γ Ω c k represents the export rate of cells of clonotypes already represented by size-k clones in the periphery. The addition of a new cell to a clone of k cells decreases by one the number of clones with k cells and increases by one the number of clones with k + 1 cells. Similarly, γ Ω c k − 1 represents the rate at which clones move from c k − 1 to c k due to thymic export. We assume that cells immigrate, proliferate, or die one cell at a time, forcing clones to move only among adjacent compartments. Thus, the term γ Ω c k − 1 − c k fully accounts for changes to c k due to thymic export. The term p(N)kc k rate at which clones move from c k to c k + 1 due to peripheral proliferation. Analogously, p(N)(k − 1)c k − 1 denotes the rate at which clones enter c k from c k − 1 due to proliferation, so that p(N)[(k − 1)c k − 1 − kc k ] accounts for changes to c k due to proliferation. The death term in Eq 3, given by μ(N) [(k + 1)c k + 1 − kc k ], is defined analogously. Eqs 2 and 4 represent "boundary conditions" for Eq 3. In Eq 2, the term Ω − ∑ i = 1 M c i gives the number of clonotypes unrepresented in the periphery, so that M c i provides the rate at which new clones enter the periphery from the thymus. Equation 2 also retains the terms from Eq 3 that account for loss of clones in c 1 due to thymic export, proliferation, and death, and the addition of clones into c 1 due to death of cells in c 2 . Finally, Eq 4 retains terms accounting for the introduction of clones into c M via thymic export to and proliferation within clones in c M − 1 , as well as loss of clones from c M due to cellular death. This represents a "boundary condition" that prevents clones from surpassing size M.
Summing Eqs 2-4, we find that so that the ODE satisfied by N(t) (Eq 1) and that satisfied by ∑ k = 1 M kc k (t) (Eq 5) differ by p(N)Mc M + γ 0 Ω c M , and thus N(t) ≠ ∑ k = 1 M kc k . Thus, the very few numbers of clones of large sizes k > M, whose population is accounted for in Eq 1, are not accounted for fully in Eqs 2-4. This is especially salient in the γ → 0 + limit, in which the N > 0 fixed point is completely "missed" by the c k = 0 solution to Eqs 2-4. Here, the γ → 0 + limit of the truncated system represents a singular limit where the only solution to Eqs 2-4, c k → 0 + , appears to violate the N = ∑ k = 1 ∞ kc k > 0 constraint at the stable fixed point.
Nonetheless, we can use the ODE system Eqs 2-4 to analyze the effects of changes in the thymic output rate γ provided we carefully use the solutions c k and N that are consistent with the N = 0 and N > 0 fixed points. We denote the normal level of thymic export in an adult of a given age by γ = γ 0 . To represent diminished thymic activity during atrophy, we take γ ≪ γ 0 , or even γ = 0, depending on the severity of the atrophy. As the thymus is highly plastic, the changes in γ throughout the process of atrophy and recovery tend to be rapid. With this in mind, we model such cycles of disease with a piecewise ODE system. Specifically, let us observe a human's response to disease-induced changes in thymic activity over some time interval I = [t 0 , t S + 1 ]. We assume that this individual's thymic export rate undergoes S abrupt changes, at times t 1 , t 2 ,..., t S , where t 0 < t 1 < t 2 < … < t S < t S + 1 . Letting I i = [t i , t i+1 ], so that I = ∪ i = 0 S I i , we assume that γ = γ i ≥ 0, on I i .
If the initial condition c k t 0 k = 1 represents the time evolution of the c k compartments after a transition to a thymic activity level γ i . This is the most general description of our model; in practice, we will typically take S = 1, representing a single abrupt change in γ(t). Further descriptions of the piecewise ODE formulation specific to certain disease patterns and particular initial conditions are included in the relevant sections below.
Linear analysis of Eqs 2-4 will provide information on how the clone counts c k (t) evolve from their steady-state values after a small perturbation in the system (through changes in γ). The dynamics of c k are not equivalent-but qualitatively related to-those of n i (t), the number of cells in clone i. For example, when large values of n i (t) increase, c k ≈ n i (t) decreases while c k ≈ n i + 1 (t) increases. Thus, increases in large n i convect c k forward, especially for larger k. As will be explicitly shown, since c k is typically monotonically decreasing in k, the dynamics of small(large) clone populations are correlated with the dynamics of c k for small(large) k.

Analysis for γ > 0 (functioning thymus)
We begin by studying the behavior of solutions of our ODE model under the assumption of a strictly positive thymic export rate, γ > 0. We perform an analysis of equilibrium solutions of Eqs 1-4 and their stability, and also compute an explicit solution in the infinite dimensional case that arises when M → ∞. At the beginning of sections 3.1 and 3.2 below, as well as in sections 4.1, 4.2 and 5, we focus on solutions over one individual interval in the piecewise formulation described in section 2 above. For simplicity when doing this, we omit the i notation that distinguishes the different subintervals, writing γ instead of γ i , etc. When the discussion returns to the full piecewise ODE, the i notation is reintroduced.

Analytic solution of the infinite dimensional system
We begin by computing analytic expressions for the solutions c k of Eqs 1-4. If we take M → ∞ and consider instead the infinite dimensional system, the c k compartments can be obtained through a generating function in a conjugate variable q, defined as Note that ∂ k Q/ ∂q k q = 0 = k!c k , which allows us to extract c k with k ≥ 0. In addition, the total population can also be recovered from the generating function via the expression ∂Q/ ∂q q = 1 = ∑ k = 0 ∞ kc k = N .
In order to derive an explicit form for Q(q, t), we assume that an explicit solution N = N(t) of Eq 1 can be found, so that we may write p and μ as functions of t (p = p(t), μ = μ(t)). By substituting Eqs 2 and 3 for dc k /dt, the time derivative of Q can be expressed as The above partial differential equation can be solved analytically along any characteristic curve q(t) defined by the solutions to dq dt = − (q − 1)(p(t)q − μ(t)): where q 0 = q(0) and Along each trajectory q(t), the generating function obeys dQ (1 − q(t))Q and can be expressed as By allowing all possible initial values q 0 , we can express the full solution as The solutions c k (t) can then be extracted from Q(q, t) by taking a series expansion of Q(q, t) and identifying coefficients with c k (t) according to Eq 6. These exact solutions c k (t) will be compared with our subsequent results derived from direct numerical solution of Eqs 1-4. By verifying that the time evolution of c k under a finite dimensional formulation as in Eqs 2-4 is sufficiently close to that of c k under the infinite dimensional formulation in Eq 11, we allow the infinite and finite dimensional systems to be used more or less interchangeably. The finite dimensional formulation has the advantage of not only admitting simple, explicit steady state solutions, but also rates of convergence to steady state.

Equilibrium solution and linearization
Returning to the truncated formulation (finite M), we now study the equilibrium solution that results when taking γ > 0 in Eqs 1-4. Denote such a generic equilibrium solution by c k * (γ) k = 1 M , N*(γ) . For a given N * (γ), the c k * (γ) have the form In the discussion below, we will write N*(γ) as N* and c k * (γ) as c k * for simplicity, unless For clarity, an example of the matrix L S with M = 4 is We now apply a simplifying assumption to the matrix L s to analytically compute its eigenvalues more easily. In general, using previous estimates for γ [55] and Ω [56,57], Proof. We begin by assuming that the terms λ k S = k p N* − μ N* are themselves eigenvalues of L S , and search for their corresponding eigenvectors, y k = y k 1 , y k 2 , ⋯, y k M , 0 .
The solution of this recurrence relation is then where we let ∏ j = 1 0 (i ± j) = 1, whenever such a term appears in the above sum. (This is verified in Appendix 6.) As previously mentioned, L S − λ k S I y k i = 0 for i = 1, 2, ⋯, M − 1.
We now compute L S − λ k S I y k M , obtaining Each term in the sum above has the form p The solutions y k i , as functions of i for fixed k, are characterized by patterns of oscillatory behavior, as shown in Figure 1, which depicts numerical solutions of true and approximate eigenvalues and eigenvectors of L S . We choose p N* and μ N* to satisfy p N* < μ N* , so that at homeostatic population levels, the death rate exceeds the proliferation rate, preventing exponential growth of the population.  Increases to μ(N*) relative to p(N*) also cause an intensified dampening of the oscillations in the eigenvector y k at lower components j, which creates the illusion of the oscillatory mass shifting to the left as μ(N*) increases for fixed p(N*), as in Figure 2b. At the same time, the entire eigenvalue spectrum becomes more negative as μ(N*) increases, as indicated in Figure 2a, so that increases to the death rate at homeostatic levels indicate much faster convergence to equilibrium of all c k compartments.

Behavior of the linearized and fully nonlinear systems
This section addresses the convergence behavior of solutions in the presence of a positive thymic export rate γ > 0. This situation represents a functioning thymus, with the possibility for many different levels of functionality, ranging from total health high γ γ 0 to dramatically diminished functionality (low γ). It could also represent a new thymic export rate after transplant of thymic tissue, as studied in the context of DiGeorge's syndrome in Ciupe et al. [40]. In this case, we determined that for each equilibrium solution of Eq 1, the system has an equilibrium solution given by Eqs 12 and 13. If the steady state solution N * (γ) > 0 of Eq 1 represents a stable fixed point, the corresponding equilibrium c k * (γ) will also be stable (typical regulated forms of proliferation and death, p(N) and μ(N), tend to result in one positive stable equilibrium solution in Eq 1). If for some i, γ i > 0, the can be discerned from the relationship between p N i t i − μ N i t i and p N* γ i − μ N* γ i , the disparity in proliferation and death rates at the starting and terminal population levels.
If these quantities differ significantly, solution trajectories are generally characterized by a transient period of fast convergence, which carries the trajectory close enough to the stable equilibrium that convergence rates from then on are dictated by the linearized eigenvalues. For example, assume that an abrupt drop in thymic productivity occurs at t i , so that γ i < γ i − 1 . As the naive T cell population evolves from N i t i to N* γ i , for which 0 > p N* γ i − μ N* γ i > p N i t i − μ N i t i , the T cell pool will experience a brief period of higher cellular death. As N(t) approaches N* γ i , the convergence rates correspond to the eigenvalues found from the linearized approximation.

Analysis for γ = 0 (full thymic cessation)
We now proceed to study the system after thymic export is shut off (γ = 0). As in section 3 above, we compute equilibrium solutions of the truncated system (finite M) that arise when γ = 0, and identify the rates of convergence of the different c k (t) to equilibrium under the linearized model. We also take M → ∞ and consider explicit solutions of the infinite-dimensional system.

Analytic solutions
In the γ = 0 case, the solution c k (t) for M → ∞ can be readily expressed using the method of characteristics. By using the generating function Q(q, t) defined in Eq 11 and taking the k-th order derivative of Q with respect to q at q = 0, we find j (18) and N(t) = A −1 (t)∑ k = 0 ∞ kc k (0) . Note that for depleted initial conditions c 0 (0) = Ω and c k (0) = 0 for k ≥ 1, Eq 18 leads to c 0 (t) = Ω and c k (t) = 0 for k ≥ 1 at all times.
Indeed, the T cell pool is expected to remain empty since there is no thymic export. To compare the infinite dimensional system with the truncated system, we also compute solutions, y k , of the truncated Eqs 2-4 (not pictured). Figure 3b depicts the relative error, c k − y k /c k . As we see, the error is several orders of magnitude smaller than c k , y k themselves at each of the times t = 30, 60, 90, indicating that the numerical solution of the truncated system Eqs 2-4 is accurately described by the exact method-of-characteristics solution and that the infinite-and finite-dimensional systems may be used more or less interchangeably.

Equilibrium solutions and linearization
We now investigate the solution c k (t) near the fixed points that arise when we take γ = 0 in assumed to be non-increasing and non-decreasing, respectively. Therefore, we have that c k * = 0 and N* = N > 0 is a stable equilibrium solution.
If γ = 0, the solution c k (t) k = 1 M , N(t), will evolve away from the unstable equilibrium c k * = 0, N* = 0 and towards the equilibrium c k * = 0, N* = N > 0. In this instance, the pool of low-population clones is eradicated due to lack of thymic productivity, and the high lymphopenic proliferation rate pushes existent clones past the truncation threshold M, where they are no longer accounted for in c k but are accounted for in N, causing N(t) → N* despite the fact that c k (t) → 0 for all k. As before, we wish to explore further the rates at which individual functions c k diverge from the unstable fixed point towards the stable one under the linearized and fully nonlinear models. To this end, we study the eigenvalues of the linearization matrix, L U , evaluated at the two equilibria.
First, consider the eigenvalues of L U evaluated at the unstable equilibrium. In this case, we assume p(0) > μ(0), as described earlier. Without thymic export, new clones are not generated in the periphery, and existent clones expand due to the high proliferation rate. Under the dynamics described by Eqs 2-4, clones quickly expand beyond the small-k compartments and get "caught" at the boundary at size M, before depleting due to the slow death-induced passage of single cell clones through the boundary at k = 1. According to Eq 1, the total cell population reaches a natural homeostatic level through peripheral maintenance alone. To investigate the rates at which these processes occur under the  , we let z 0 be defined by solutions to the recurrence relation, for i = 1, 2, …, M − 2. It can be directly verified by induction that the solution to this recurrence relation is The eigenvalue λ 1 U of second smallest magnitude is well separated from λ 0 U , and it encodes information about how the number of small clones evolves. Similar analysis of an eigenvalue λ 1 U and eigenvector z 1 approximating λ 1 U and z 1 indicates that c k empties much more rapidly for small k than it does for large k, as small clones expand in size and race to the boundary at k = M. In particular, all c k except those with k ~ M, which had been preserved by the slow eigenvalue λ 0 U , empty at nearly the same rate, λ 1 U ≈ λ 1 U = (μ(0) − p(0)) . i + 1 be given by the solutions to the following recurrence relation: It is worth nothing that if we were to instead choose

Author Manuscript
Author Manuscript

Author Manuscript
Author Manuscript as M → ∞. Thus, for M ≫ 1, L U − λ 1 U I z 1 ≈ 0, and we find that z 1 is "almost" an eigenvector of L U corresponding to the approximate eigenvalue λ 1 U . □ We now wish to identify which of the c k empty at the rate determined by the second approximate eigenvalue, λ 1 U . As it turns out, the eigenvector z 1 corresponding to this eigenvalue is "nearly" constant, and thus all c k empty at essentially the same rate. We see this by identifying that even though we are only concerned with a finite number (M) of terms of the sequence generated by the recurrence relation in Eq 23, as the index M becomes infinitely large, the sequence z 1 i i = 1 M exhibits "Cauchy-like" behavior, mimicking "convergence" to a limiting value. We make this more precise in the following proposition: allowing us to conclude that z 1 m − z 1 n < ε for all M ≥ m, n ≥ cM. □ By taking 0 < c, ε ≪ 1, we find that "most" components of  The analysis conducted in section 3 does not apply, as it relied on the assumption p(N*) < μ(N*), which no longer holds. However, numerical computation indicates that in the case p(N*) = μ(N*), which applies here, the eigenvalue spectrum and associated eigenvectors qualitatively resemble the λ k S , y k studied analytically in section 3.

Behavior of the linearized and fully nonlinear systems
We now interpret these results in the context of the particular diseased states to which they naturally apply. We first identified an unstable equilibrium solution, c k * = 0, N* = 0, and studied the linearization of the system around this equilibrium. Under the linearized model, if γ i = 0 for some i, the eigenvalue/eigenvectors pairs suggest that solutions diverge away from this equilibrium, with c k i for small k evolving at a rate λ 1 U = (μ(0) − p(0)), and c k i for k − M evolving at the very slow rate given by the small-magnitude eigenvalue, λ 0 U .
The total population N i (t) evolves at the rate (p(0) − μ(0)). This situation represents the repopulation of the T cell pool from a small number of cells via peripheral proliferation in a highly pathological state involving both complete thymic inactivity (e.g. thymectomy or total functional cessation) and near full lymphopenia (as may result from treatment regimens for cancer, etc.).
We then identified a stable equilibrium solution, c k * = 0, N* = N > 0. As N is asymptotically stable, N i (t) N after diverging from N* = 0. Under the linearized model, the eigenvalue/ eigenvector pairs suggest that c k i (t) 0 slowly for small k, and c k i (t) 0 much more quickly for large k.
As before, the validity of the eigenvalues in providing accurate convergence rates of c k i , N i to and from equilibria depends on the initial condition c k i t i in the full nonlinear model. If the human is in a state of immune health for t < t i , so that γ i−1 > 0 and the initial conditions c k i t i , N i t i > 0 satisfy c k i t i c k * γ i − 1 , N i t i N* γ i − 1 , we expect that p N i t i < μ N i t i . The higher rate of death than proliferation at t i may cause a transient period of quick collapse, with N i (t) decreasing to N . As N i (t) N, convergence occurs at the rates dictated by the linearized eigenvalues. If γ i−1 > 0 but c k i t i , N i t i 0, so that the thymus is functioning to some extent but the T cell pool has been eradicated, trajectories first diverge away from the unstable zero equilibrium at rates given by the linearized eigenvalues. As p N i t i − μ N i t i p(N) − μ(N) = 0, the motion of trajectories transitions from being dictated by the eigenvalues of the unstable equilibrium to those of the stable equilibrium.
In summary, we show in Proposition 4.1 that c k for larger k is sustained by a near-zero eigenvalue λ 0 U ≃ 0, as the solution evolves away from the unstable empty state when γ = 0.
Furthermore, in Propositions 4.2 and 4.3 we identify a series of negative eigenvalues, and a uniform eigenvector corresponding to the slowest decay rate λ 1 U < 0. It suggests a uniform asymptotic decay of components other than the larger k components preserved by λ 0 U .

Special cases and numerical evaluation
Using the approximate rates of convergence provided by linearization, we can now study the time scale of the T cell pool's adjustment to a new export rate. While some T cell clones will expand and attain a large size, most are small. Thus, in both the cases of thymic atrophy and recovery, we take as a proxy for the rate at which the T cell diversity converges to equilibrium the eigenvalue that dictates the rate of convergence of c 1 , typically given by the quantity p(N*) − μ(N*) (this tends to also be the dominant eigenvalue). Recalling that p(N*) < μ(N*), high proliferation rates (p(N*) ~ μ(N*)) lead to small values of |p(N*) -μ(N*)|, thus slower adjustment of diversity to the changing γ. That is, in a proliferation-dominant scenario, a drop in thymic export leads to repopulation via clonal expansion. In this section, we study several specific models arising from canonical choices of p and μ, and compute the changing convergence rate as gamma varies.

The logistic mode: Regulated proliferation, constant death
We In Figure 5, the quantity λ 1 S is plotted against γ for several different combinations of p 0 ,μ 0 , showing the unboundedness of the convergence rate as γ increases. Within this physiological range of γ values, the dependence of λ 1 S on γ presents as linear on a log-log plot, indicating a power law relationship. Indeed, the power law is described in detail in the caption of Figure 5.

Constant proliferation, regulated death
Let us now assume that p(N) = p 0 > 0 and μ(N) = μ 0 + μ 1 N 2 K 2 + N 2 , with μ 0 , μ 1 > 0, in the determination of N(t) via Eq 1. We assume that p 0 > μ 0 and p 0 − μ 0 + μ 1 < 0, so that the action of the proliferation-death mechanism results in net cellular birth at low cell counts and net cellular death at high cell counts. The steady states of Eq 1 are given by the roots of the following cubic P (N) = p 0 − μ 0 + μ 1 N 3 + γN 2 + p 0 − μ 0 K 2 N + γK 2 . (30) First note that P(0) = γK 2 > 0, and the highest order coefficient satisfies (p 0 − (μ 0 + μ 1 )) < 0 by assumption, so that P(N) → −∞ as N → ∞, and P(N) has at least one positive, real root. From Descartes' rules of signs, the polynomial has at most one positive real root, so we may conclude that it has precisely one positive real root. This root corresponds to the only physically relevant stable fixed point of  From the fact that N* → ∞ as γ → ∞, we have that |p(N*) − μ(N*)| → p 0 − (μ 0 + μ 1 ) as γ → ∞. This limiting behavior is reflected in the eventual plateau seen in Figure  6, which plots the quantity λ 1 S in this case. Before the plateau occurs, γ and λ 1 S are again related by a power law. The transition from power law to plateau occurs at a "threshold" value, γ*, of γ, at which the rate of T cell adjustment becomes sensitive to a changing thymic export rate. If, for some i, γ i − 1 , γ i ≥ γ*, λ 1 S is unaffected by the transition from thymic export rate γ i−1 to thymic export rate γ i -that is, the T cell pool adjusts to the new thymic export rate γ i as quickly as it had adjusted to the previous thymic export rate γ i−1 . If, however, γ i − γ* γ i − 1 − γ* < 0, then a dramatic shift in the adjustment rate will occur. Thus, parameter choices that result in a low threshold value γ* might correspond to physiological conditions under which an instance of acute thymic atrophy actually does not affect T cell adjustment rates. Likewise, a high threshold value of γ* indicates potential sensitivity of adjustment rates to the changing level of thymic export, with adjustment rates obeying a power law dependence on γ.

Discussion and conclusions
In this paper, we formulated a model of how the naive T cell pool adjusts to changes in the rate of thymic export of new T cells during a cycle of stress-induced atrophy and recovery, and how it may be reconstituted following an instance of severe lymphopenia induced by a state of immune disease, or treatments such as chemotherapy. Our underlying model is a birth-death-immigration process studied under a mean-field approximation for the mean clone abundance distribution (or clone count) c k (t). A recent investigation into the fully stochastic model indicates that the true c k differs from the c k derived using the mean-field assumption (Eqs 2-5) only for very large k ≈ N [52]. Thus, our analyses may be inaccurate only if a single large clone dominates the whole population.
Another modeling choice we made is that TCRs are generated one naive T cell at a time. Successive emigrations from the thymus are uncorrelated with the TCRs that are produced. However, emigration can be clustered, with cell proliferation generating Δ ~ 2 − 4 copies of naive T cells carrying the same TCR during each emigration event. In this case, we simply modify the immigration terms in Eqs 2-3. For Eq 2, the immigration term proportional to γ/Ω is removed, while the immigration term γ/Ω c k − 1 − c k is replaced by γ/Ω c k − Δ − c k in Eq 3. By setting c l < 0 = 0, the solution to Eqs 2 and 3 can be numerically evaluated but a closed-form analytic solution is not possible. Solutions with clustered immigration (Δ > 1) show no qualitative difference from Δ = 1, with minor quantitative differences arising only for very small k.
In section 3, we found that our mean-field ODE model admitted one stable equilibrium solution when γ > 0. From an analysis of the eigenvalues and eigenvectors of the system linearized around this stable equilibrium, we found that for small k, perturbations in c k about a steady-state solution are weighted more strongly in the slowest mode (slowest eigenvalue) λ 1 S = p N* − μ N* < 0. As also shown in Figure 1, the variation in c k for larger k contains higher weights of the faster modes corresponding to more negative (faster) eigenvalues λ l S = l p N* − μ N* < 0. Similarly, in section 4, we analyzed the eigenvalue and eigenvector decomposition of the solution for γ = 0, for which two equilibrium points, N* = 0 and N* > 0, arise. For p(N*) − μ(N*) < 0, we find the same decomposition of c k (t) as in the γ > 0 case in section 3. Thus, in terms of the relaxation of c k (t) towards a finite steady-state, our eigenvalue/eigenvector analysis suggests that the counts of large clones might evolve faster towards the new steady-state. For the unstable equilibrium state N* = 0, the eigenvalue/eigenvector decomposition of c k is shown in Figure 4. In this unstable case, the slowest eigenvalue λ 0 U ≈ 0 has a corresponding eigenvector z 0 i with elements that decay with i. This result predicts that large-population (M − i) clones relax slowly (recall that the labeling is inverted: i = M − 1 corresponds to c 1 ).
In addition to decomposing the linearized solutions in terms of eigenvalues and eigenvectors, we explicitly plot trajectories of c k (t) following small, abrupt changes in γ.
By plotting the log of the deviation δc k = c k (t) − c k (∞) /c k (∞) in all cases (Figure 7a), we can see that counts of small clones evolve faster after perturbation. Although this seems to contradict the eigen decomposition of the γ > 0 case, the coefficients and eigenvector elements corresponding to larger l are negative (see Figure 1b,d), and convert fast decaying modes into short bursts of growth and conspire to cancel the fast dynamics indicated by the more negative eigenvalues. In fact, we see that the clone counts of large clones actually evolve more slowly than the rate associated with the largest eigenvalue. This behavior can also be heuristically understood by noticing that under a given change in γ, the immigration term γ c k − 1 − c k /Ω is larger for smaller k because c k ~ 1/k. Thus, for a given change in γ, the perturbation is larger in the equations with smaller k and induces changes in c k (t) that appear larger.
In section 5, we infer the rate of convergence of the counts of the smallest (but most common) clones to equilibrium by computing the dominant eigenvalue as a function of γ > 0 for two choices of regulated functions p(N),μ(N). In section 5.1, we test the logistic model, which assumes a constant death rate but a population-dependent proliferation rate.
From the explicit form of λ 1 S in Eq 29, p N* − μ N* ∞ as γ ∞ for fixed values of the other parameters, which produces the power-law relationship between γ and λ 1 S depicted in Figure 5. In section 5.2, we assumed instead a constant rate of cellular proliferation with an N-dependent death rate. In this case, which differs from the logistic formulation in that regulation is incorporated into μ(N) via a Hill-type function, γ and λ 1 S are related by a power law for low γ, before reaching a plateau at higher γ.
Since regulation through death is typically associated with the actual mechanism of naive T cell survival, we expect this mechanism to be more realistic than a population-dependent proliferation rate p(N). Thus, jumps in the thymic export rate that either cross the threshold value of γ, or occur between two values of γ both in the power-law region, can be expected to produce changes in both the equilibrium values of c k and also the convergence rates. If a jump in γ occurs between two values of γ that are both in the plateau region, the equilibrium values shift, but the convergence rates stay the same. The presence of the power law region indicates the robustness of the T cell diversity during a time of severe thymic atrophy. That is, the slower convergence of the T cell diversity to equilibrium at low γ values protects the pool from quick shifts to the lower diversity associated with lower γ values.     shows that clone counts of small clones appear to evolve faster than counts of larger clones. (b) Similarly, the number of small clones also evolve faster away from an unstable empty equilibrium state N* = 0.