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-15]. Acute atrophy can occur under a plethora of conditions associated with a state of disease or stress [16-18], including viral, bacterial and fungal infection [19-21], malnutrition [22], cancer treatment [23], bone marrow transplant [24], psychological stress and pregnancy [25-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-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-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-51]. At steady-state, when a healthy, homeostatic cell count N* is achieved, . Note that when γ = 0, simple decreasing functions p(N) (and/or 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(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 the marginalized probability that any single clone i has population k. The master equation for P(n, t) is difficult to solve with regulation terms. Here, we provide a heuristic derivation of the equations obeyed by c(t) by using a mean-field approximation that the total population N is uncorrelated with any n, (although we know ). Under this mean-field approximation, the evolution of P(n = k, t), and hence c(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 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 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, represents the rate at which clones move from c to c 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 fully accounts for changes to c due to thymic export. The term p(N)kc denotes the rate at which cells in size-k clones proliferate, which in turn corresponds to the rate at which clones move from c to c due to peripheral proliferation. Analogously, p(N)(k − 1)c denotes the rate at which clones enter c from c due to proliferation, so that p(N)[(k − 1)c − kc] accounts for changes to c due to proliferation. The death term in Eq 3, given by μ(N) [(k + 1)c − kc], is defined analogously. Eqs 2 and 4 represent “boundary conditions” for Eq 3. In Eq 2, the term gives the number of clonotypes unrepresented in the periphery, so that 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 c1 due to thymic export, proliferation, and death, and the addition of clones into c1 due to death of cells in c2. Finally, Eq 4 retains terms accounting for the introduction of clones into c via thymic export to and proliferation within clones in c, as well as loss of clones from c 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 (Eq 5) differ by and thus 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 = 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 → 0+, appears to violate the 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 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 = [t0, t]. We assume that this individual’s thymic export rate undergoes S abrupt changes, at times t1, t2,..., t, where t0 < t1 < t2 < … < t < t. Letting I = [t, t], so that we assume that γ = γ ≥ 0, on I. If the initial condition represents the size of each c compartment at the start of the process, we then let represent the solution of the ODE in Eqs 1–4, on I, with γ = γ and initial condition Thus, the solution represents the time evolution of the c compartments after a transition to a thymic activity level γ. 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(t) evolve from their steady-state values after a small perturbation in the system (through changes in γ). The dynamics of c are not equivalent–but qualitatively related to–those of n(t), the number of cells in clone i. For example, when large values of n(t) increase, decreases while increases. Thus, increases in large n convect c forward, especially for larger k. As will be explicitly shown, since c is typically monotonically decreasing in k, the dynamics of small(large) clone populations are correlated with the dynamics of c 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 γ, 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 of Eqs 1–4. If we take M → ∞ and consider instead the infinite dimensional system, the c compartments can be obtained through a generating function in a conjugate variable q, defined as

Note that which allows us to extract c with k ≥ 0. In addition, the total population can also be recovered from the generating function via the expression

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/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 where q0 = q(0) and

Along each trajectory q(t), the generating function obeys and can be expressed as

By allowing all possible initial values q0, we can express the full solution as

The solutions c(t) can then be extracted from Q(q, t) by taking a series expansion of Q(q, t) and identifying coefficients with c(t) according to Eq 6. These exact solutions c(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 under a finite dimensional formulation as in Eqs 2–4 is sufficiently close to that of c 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 For a given N*(γ), the have the form

In the discussion below, we will write for simplicity, unless desiring to emphasize the γ-dependence. To identify the stability of this equilibrium solution–and the rates of convergence of solutions to equilibria under the linearized model later on–we consider the linearization of the system around this generic equilibrium solution, represented by the with component s given by

For clarity, an example of the matrix L with M = 4 is

We now apply a simplifying assumption to the matrix L to analytically compute its eigenvalues more easily. In general, using previous estimates for γ [55] and Ω [56,57], (when rates are measured in units of year−1), thus, we assume and define a new matrix which is L but with only in elements in which it appears explicitly. is defined by Eq 14 with in the first five terms set to zero, but with N* determined under the appropriate general immigration rate γ ≥ 0.

Numerical computation confirms that the eigenvalues of the new matrix are essentially identical to those of the original matrix validating our assumption that the term may be neglected in L. Note that this assumption does not cause us to omit the constant γ from the linearization matrix L entirely, as the steady state values depend on γ. Denote by the eigenvalues of , and note that that the entry is an eigenvalue. We denote this eigenvalue by Since N* is the stable equilibrium solution of Eq 1, we assume that As shown in [58], the eigenvalues of all have strictly negative real part, as do those of L. The eigenvalues of the minor of L can be approximated by with corresponding eigenvectors shown in the following proposition:

Proposition 3.1.

If the eigenvalues of the matrix are well approximated by the terms in the sense that there exist vectors such that

Proof. We begin by assuming that the terms are themselves eigenvalues of , and search for their corresponding eigenvectors, Choosing we then choose for inductively so as to force the i-th component of the residual vector, which we denote by to equal zero for We then verify that so that for where ||·|| is any p-norm. (Trivially, We first note that the components of the approximate eigenvector corresponding to eigenvalue are defined by the recurrence relation,

The solution of this recurrence relation is then where we let whenever such a term appears in the above sum. (This is verified in Appendix 6.) As previously mentioned, We now compute obtaining

Each term in the sum above has the form is a polynomial of degree k in the variable M, and a = p(N*)/μ(N*). Recalling that p(N*)/μ(N*) < 1, then so that at large values of This demonstrates that may be regarded as an approximation to the true eigenvalue assuming that the eigenvalues of are stable under small perturbations. □

The solutions 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 We choose so that at homeostatic population levels, the death rate exceeds the proliferation rate, preventing exponential growth of the population.

Figure 1.

Eigenvalues and eigenvectors of Here, and in all subsequent evaluations, we parameterize our model using values in a range based qualitatively on human data, where the unit of rates are expressed in 1/year [55,59,60]. However, the model can have arbitrary units in general cases. (a) Numerically computed eigenvalue spectrum of the matrix with . Dots identify the locations of the eigenvalues (b) First 100 components of the eigenvectors with indices k = 1, 100, 200, the eigenvalues corresponding to which are marked on the spectral curve in (a). (c) Comparison of true eigenvalues and approximate eigenvalues . Approximation is strong for . (d) (Top) comparison of showing that the approximation is strong. (Bottom) comparison of showing that the accuracy of the approximation breaks down, but the qualitative behavior of is captured in

Figure 1a presents a sample plot of the eigenvalue spectrum in the case M = 500; on the spectral curve, several eigenvalues are marked, for which the first 100 components of the corresponding eigenvectors are plotted in Figure 1b. As discussed previously, all eigenvalues are negative. The eigenvector plots in Figure 1b indicate that as k increases, the oscillatory “mass” of the corresponding eigenvectors occurs at increasingly large values of j. Figure 1c presents a comparison of the true eigenvalue spectrum and the approximate eigenvalue spectrum of the matrix The eigenvalue approximation is very strong for and this remains true as M varies. The quantities over-approximate the . Figure 1d depicts a comparison of for two values of k, one below the crossover ~ M/10 (k = 50, top) and one above M/10 (k = 200, bottom). As expected, the eigenvector approximation is accurate precisely when the corresponding eigenvalue approximation is accurate. Even for the approximate eigenvectors present an appearance similar to that of the for smaller k. The diminished accuracy of the eigenvalue/eigenvector pairs at higher k (for fixed M) is attributable to the slower convergence of to 0 as M → ∞ for larger k, which is immediately apparent from the form in Eq 17. For a system of a fixed dimension M, increasing μ(N*) relative to p(N*) causes the approximate eigenvalues to become increasingly accurate at larger k, clearly due to the quicker convergence of the residual quantity when p(N*)/μ(N*) ≪ 1, as indicated by Eq 16.

Increases to μ(N*) relative to p(N*) also cause an intensified dampening of the oscillations in the eigenvector 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 compartments.

Figure 2.

Eigenvalues and eigenvectors of varying μ. (a) Numerically computed eigenvalues, when μ(N*) = 0.16 and μ(N*) = 0.6. In both cases, p(N*) = 0.15, M = 100. (b) Numerically computed eigenvectors y81.

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 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 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, γ > 0, the solution The convergence of the total population occurs at rate Based on the eigenvalues and eigenvectors of the approximate linearization, of Eqs 1–4 around this equilibrium, we can formally construct approximations to the time-dependent solutions for c(t).

In the fully nonlinear system, however, the linearized eigenvalues provide only a priori rates of convergence of solution trajectories initialized near equilibrium. The accuracy of the eigenvalues in providing convergence rates of solutions depends on the initial conditions. If the initial conditions, then the solutions begin near the stable equilibrium, and the eigenvalues provide accurate rates of convergence. If the initial conditions are far from equilibrium, the eigenvalues may not provide accurate rates of convergence of the entire solution trajectory. When trajectories are far from equilibrium at time t, further information about the speed of convergence can be discerned from the relationship between 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 As the naive T cell population evolves from for which the T cell pool will experience a brief period of higher cellular death. As N(t) approaches 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(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(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 and Note that for depleted initial conditions Eq 18 leads to at all times. Indeed, the T cell pool is expected to remain empty since there is no thymic export.

Figure 3a depicts a numerical computation of the solution, c, of the infinite-dimensional formulation, obtained from the generating function. We include values of c for k = 1, 2,…, 50 at times t = 30, 60, 90. As a function of k, the c present as linear on a logarithmic scale, as expected. To compare the infinite dimensional system with the truncated system, we also compute solutions, y, of the truncated Eqs 2–4 (not pictured). Figure 3b depicts the relative error, As we see, the error is several orders of magnitude smaller than c, y 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.

Figure 3.

Computation of c from method of characteristics, comparison with truncated system. (a) Plots of at times t = 30, 60, 90. Solutions c were computed numerically from the analytic method described in 4.1, based on the infinite-dimensional system. As a function of k, c presents as linear on a logarithmic scale. (b) Relative error, at times t = 30, 60, 90, where c denotes the solutions depicted in (a), and y denotes the numerically computed solutions of the truncated system in Eqs 2–4. From (b), the disparity between the solutions of the infinite dimensional systems (c) and the finite dimensional truncated systems (y) is negligible, validating our decision to use them interchangeably. Coefficient functions: Parameter values: Initial condition: .

Equilibrium solutions and linearization

We now investigate the solution c(t) near the fixed points that arise when we take γ = 0 in Eqs 1–4. In this case, the system has two possible equilibrium solutions. Denoting generic equilibrium solutions by the unstable solution is for all k ≥ 1 and N* = 0, and the asymptotically stable solution is for all 1 ≤ k ≤ M. However, in the stable state we define where satisfies To verify the stability of these solutions, we consider the linearization of the system around this equilibrium, which is represented by the matrix we call The components u of L are given explicitly by:

As before, is an eigenvalue with eigenvector (0,...,0,1), and the remaining eigenvalues are those of the minor of L. All eigenvalues of the (M + 1) × (M + 1) minor have negative real part, and the stability of an equilibrium solution depends on the sign of If N* = 0, then as described previously, and the equilibrium is unstable. On the other hand, if then N* represents a positive, homeostatic cell count, and as p(N), μ(N) are assumed to be non-increasing and non-decreasing, respectively. Therefore, we have that is a stable equilibrium solution.

If γ = 0, the solution will evolve away from the unstable equilibrium and towards the equilibrium 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 but are accounted for in N, causing N(t) → N* despite the fact that c(t) → 0 for all k. As before, we wish to explore further the rates at which individual functions c 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, evaluated at the two equilibria.

First, consider the eigenvalues of L 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 linearized model, we derive approximations to the dominant eigenvalues of LU. Under the assumption that p(0) > μ(0), we denote the true eigenvalues of for k = 0, 1,…M, with corresponding eigenvectors (note that we have reversed the index ordering here). Assign to the eigenvalue the label and to its eigenvector (0,...,0,1) the label z. What remains is to find approximations to the other M eigenvalues of L, which are precisely the eigenvalues of the (M + 1) × (M + 1) minor. For k = 0,1,...,M − 1, denote the approximation to the eigenvalue and the approximation to the eigenvector We begin by establishing that the eigenvalue of the (M + 1) × (M + 1) minor with the smallest magnitude, is well approximated by

Proposition 4.1.

The eigenvalue of L, is well approximated by in the sense that there exists a vector

Proof. (Note: The components of are written above in “descending” order for notational convenience, as this reflects the order in which they will be chosen recursively below. The i-th component from the left of , denoted explicitly by still corresponds to the function c.) We begin by considering the matrix and searching for an appropriate eigenvector, We define and once again choose the components inductively via a three-term recurrence relation so as to force the i-th component of which we denote as before by to satisfy for i = 2, 3, …, M + 1. While we show that to, so that may be regarded formally as an “approximate” eigenvector corresponding to the approximate eigenvalue

Defining 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

By construction of the recurrence relation, The first component, satisfies as M → ∞. Thus, when and we may conclude that are suitable approximations to respectively. □

Recalling that all eigenvalues of L have negative real parts, the true eigenvalue has a negative real part of very small magnitude. We now identify which entries in the eigenvector z0, as approximated by are particularly large in magnitude in comparison with the others. Recalling that the i-th component of z0 is given by decay nearly exponentially in i, so that c for large k are preserved by this slow eigenvalue, in agreement with the described “build up” of clones at the boundary k = M when γ = 0.

The eigenvalue of second smallest magnitude is well separated from and it encodes information about how the number of small clones evolves. Similar analysis of an eigenvalue and eigenvector approximating indicates that c 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 except those with k ~ M, which had been preserved by the slow eigenvalue empty at nearly the same rate,

Proposition 4.2.

In the case p(0) > μ(0), the matrix L has an eigenvalue, which is well approximated by in the sense that there exists a vector such that

Proof. First define Then for i = 1, 2, ... , M − 2, let be given by the solutions to the following recurrence relation:

It is worth nothing that if we were to instead choose then the recurrence relation in Eq 23 would have a constant solution, Although for our purposes, solutions of the recurrence relation do converge rapidly to constants, as will be discussed later.

By construction, Additionally, To show that we use Eq 23 to derive an explicit bound on the quantity Consider

Taking i = M − 2 in the above relation, we find as M → ∞. Thus, for M ≫ 1, and we find that is “almost” an eigenvector of L corresponding to the approximate eigenvalue □

We now wish to identify which of the c empty at the rate determined by the second approximate eigenvalue, As it turns out, the eigenvector corresponding to this eigenvalue is “nearly” constant, and thus all c 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 exhibits “Cauchy-like” behavior, mimicking “convergence” to a limiting value. We make this more precise in the following proposition:

Proposition 4.3.

Let be an approximate eigenvector of L corresponding to approximate eigenvalue where is generated by the recurrence relation in exhibit “Cauchy-like” behavior: for any ε > 0 and 0 < c < 1, we may choose M ∈ N such that for all cM ≤ m ≤ M and cM ≤ n ≤ M.

Proof. Recalling the bound on obtained in Eq 24, we find

Recalling that we may choose M large enough that allowing us to conclude that for all M ≥ m, n ≥ cM. □

By taking we find that “most” components of are within an ε-distance of each other, so that the eigenvector is nearly constant. (The constancy breaks down at the components representing large-k compartments.) With this, we are able to classify the rates at which all of the c are lost from the pool. As the eigenvector corresponding to the second smallest magnitude eigenvalue, is nearly constant, all c, except those with k ~ M, empty at nearly the same rate, on a time scale

The remaining eigenvalues are not treated analytically, but numerical computation indicates that a similar general approximation to the k-th eigenvalue, may be made, taking the form The true eigenvalues, are depicted in Figure 4a (comparison of the true and approximate spectra is omitted, as the result is similar to that depicted in Figure 1c. That is, The oscillatory behavior observed in the approximate eigenvectors in the case p(N*) < μ(N*) of section 3 is absent here; although the subsequent approximate eigenvectors do not share the Cauchy-like behavior of the components corresponding to c for small k do not vary much in magnitude, and are thus interpreted as being nearly constant themselves (Figure 4b). Thus, we conclude that for small k, the functions c all have very similar dynamics, diverging at the rate while for large k, the c converge very slowly, at a rate governed by the dominant, near-zero eigenvalue

Figure 4.

Eigenvalues and eigenvectors of (a) Numerically computed eigenvalue spectrum of the matrix L, with p(0) = 0.17, μ(0) = 0.12, and M = 500. Dots identify the locations of the eigenvalues The first several eigenvalues of the full system, including the positive eigenvalue (open circle when p(0) > μ(0)) and the eigenvalue are highlighted in the zoomed-in axis. (b) First 30 components of the eigenvectors with indices k = 1, 6, 21, the eigenvalues corresponding to which are marked on the spectral curve in (a). Note that the eigenvectors were defined in a “reverse-order”, so that corresponds to compartment to compartment c6, and to compartment c21. Generally, corresponds to compartment c.

We now consider the stable equilibrium solution, In this case, the linearization around this equilibrium may be expressed as matrix with component given by

For reference, an example of the matrix with M = 4, is given below:

Denote by the eigenvalues of the matrix evaluated at the stable equilibrium solution, and their corresponding eigenvectors. As before, let Then, the remaining eigenvalues are those of the which are independent of the parameters of the system except M. (Of course, the eigenvalues of .) These eigenvalues and eigenvectors are not treated analytically. 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 studied analytically in section 3.

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, and studied the linearization of the system around this equilibrium. Under the linearized model, if γ = 0 for some i, the eigenvalue/eigenvectors pairs suggest that solutions diverge away from this equilibrium, with for small k evolving at a rate for k − M evolving at the very slow rate given by the small-magnitude eigenvalue, The total population N(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, is asymptotically stable, after diverging from Under the linearized model, the eigenvalue/eigenvector pairs suggest that slowly for small k, and much more quickly for large k.

As before, the validity of the eigenvalues in providing accurate convergence rates of to and from equilibria depends on the initial condition in the full nonlinear model. If the human is in a state of immune health for so that γ > 0 and the initial conditions satisfy we expect that The higher rate of death than proliferation at t may cause a transient period of quick collapse, with N(t) decreasing to convergence occurs at the rates dictated by the linearized eigenvalues. If γ > 0 but 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 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 for larger k is sustained by a near-zero eigenvalue 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 It suggests a uniform asymptotic decay of components other than the larger k components preserved by

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 c1, 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 begin with the canonical logistic growth model, taking p(N) = p0(1 − N/K), μ(N) = μo, where p0, μ0 > 0 are basal rates of cellular proliferation and death, respectively, and K > 0 is an inherent carrying capacity. Under this model, Eq 1 has a positive steady state, N*, given by

In this case, so that the assumption always applies. Moreover, so it is clear that is the dominant eigenvalue. Then,

In Figure 5, the quantity is plotted against γ for several different combinations of p0,μ0, showing the unboundedness of the convergence rate as γ increases. Within this physiological range of γ values, the dependence of 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.

Figure 5.

Dominant eigenvalue, plotted against γ, case 1. In (a), p0 = 0.18, μ0 = 0.17. In (b), p0 = 0.018, μ0 = 0.017. In (c), μ0 = 0.0018, μ0 = 0.0017. There is an approximate power law relationship between within this range of parameter values. In (a), for example, the best fit line to the curve K = 107 is given by The curve K = 109 is fit by with R = 0.9752, and the curve K = 1011 is fit by with R2 = 0.9801.

Constant proliferation, regulated death

Let us now assume that with μ0, μ1 > 0, in the determination of N(t) via Eq 1. We assume that 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

First note that P(0) = γK2 > 0, and the highest order coefficient satisfies (p0 − (μ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 By regarding this root, N*, as the intersection of the line and the rational expression we see that N* → ∞ as γ → ∞.

We also verify that the eigenvalues satisfy is, in fact, the dominant eigenvalue. We first check that Recalling that we see after some simple algebraic manipulation that the condition is equivalent to

But to, along with the fact that P(N) has only one real positive root indicates that N* does, in fact, satisfy Eq 31, so that It is easily verified that and consequently that

From the fact that N* → ∞ as γ → ∞, we have that |p(N*) − μ(N*)| → p0 − (μ0 + μ1) as γ → ∞. This limiting behavior is reflected in the eventual plateau seen in Figure 6, which plots the quantity in this case. Before the plateau occurs, γ and 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 is unaffected by the transition from thymic export rate γ to thymic export rate γ-that is, the T cell pool adjusts to the new thymic export rate γ as quickly as it had adjusted to the previous thymic export rate γ. If, however, 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 γ.

Figure 6.

Dominant eigenvalue, plotted against γ case 2. In (a), p0 = 0.18, μ0 = 0.17, and μ1 = 0.004. In (b), p0 = 0.018, μ0 = 0.017, and μ1 = 0.004. In (c), p0 = 0.0018, μ0 = 0.0017, and μ1 = 0.0004. The relationship between follows a power law for low values of γ, before reaching a plateau for high values of γ. In (a), the best fit line to the curve K = 107 over the power law region is given by with R2 = 0.9989. The curve K = 109 over the power law region is fit by log γ − 8.708, with R2 = 0.9995, and the curve K = 1011 over the power law region is fit by log γ − 10.72, with R2 = 0.9998.

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(t). A recent investigation into the fully stochastic model indicates that the true c differs from the c derived using the mean-field assumption (Eqs 2–5) only for very large [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 is replaced by in Eq 3. By setting 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 about a steady-state solution are weighted more strongly in the slowest mode (slowest eigenvalue) As also shown in Figure 1, the variation in c for larger k contains higher weights of the faster modes corresponding to more negative (faster) eigenvalues 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(t) as in the γ > 0 case in section 3. Thus, in terms of the relaxation of c(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 is shown in Figure 4. In this unstable case, the slowest eigenvalue has a corresponding eigenvector 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 c1).

In addition to decomposing the linearized solutions in terms of eigenvalues and eigenvectors, we explicitly plot trajectories of c(t) following small, abrupt changes in γ.

By plotting the log of the deviation 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 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 is larger for smaller k because c ~ 1/k. Thus, for a given change in γ, the perturbation is larger in the equations with smaller k and induces changes in c (t) that appear larger.

Figure 7.

Time-dependence of c(t) near fixed points. (a) The explicit time-evolution of following a small abrupt change γ = 1.8 × 1010 → 0.99 × 1.8 × 1010 at t = 0. The evolution from one nonzero steady state to another nonzero steady state 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.

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 in Eq 29, for fixed values of the other parameters, which produces the power-law relationship between γ and 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 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 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.