How Naive T-Cell Clone Counts Are Shaped By Heterogeneous Thymic Output and Homeostatic Proliferation

Abstract

The specificity of T cells is that each T cell has only one T cell receptor (TCR). A T cell clone represents a collection of T cells with the same TCR sequence. Thus, the number of different T cell clones in an organism reflects the number of different T cell receptors (TCRs) that arise from recombination of the V(D)J gene segments during T cell development in the thymus. TCR diversity and more specifically, the clone abundance distribution, are important factors in immune functions. Specific recombination patterns occur more frequently than others while subsequent interactions between TCRs and self-antigens are known to trigger proliferation and sustain naive T cell survival. These processes are TCR-dependent, leading to clone-dependent thymic export and naive T cell proliferation rates. We describe the heterogeneous steady-state population of naive T cells (those that have not yet been antigenically triggered) by using a mean-field model of a regulated birth-death-immigration process. After accounting for random sampling, we investigate how TCR-dependent heterogeneities in immigration and proliferation rates affect the shape of clone abundance distributions (the number of different clones that are represented by a specific number of cells, or "clone counts"). By using reasonable physiological parameter values and fitting predicted clone counts to experimentally sampled clone abundances, we show that realistic levels of heterogeneity in immigration rates cause very little change to predicted clone-counts, but that modest heterogeneity in proliferation rates can generate the observed clone abundances. Our analysis provides constraints among physiological parameters that are necessary to yield predictions that qualitatively match the data. Assumptions of the model and potentially other important mechanistic factors are discussed.

Keywords: T-cell receptor; clone-count distributions; heterogeneity; immigration-proliferation model; mathematical modeling; naive T cells; repertoire diversity.

Figure 1

Normalized naive T cell clone count data from one patient in Oakes et al. (12) plotted on a log-log scale. Values of the normalized clone counts along the vertical axis are the average of three samples among CD4 and CD8 cell subgroups. Clones are defined by different nucleotide sequences associated with different alpha or beta chains of the TCR.

Figure 2

Schematic of a multiclone birth-death-immigration process. Clones are defined by distinct TCR sequences i. Each clone carries its own thymic output and peripheral proliferation rates, α_i and r_i , respectively. We assume all clones have the same population-dependent death rate μ(), where is the total number of cells in the organism that influence the death rate. Since Q ≫ 1, we impose a continuous distribution over the rates α and r. Theoretically, there may be Q ≳ 10¹⁵ (6) or more (30, 31) possible viable V(D)J recombinations. The actual, effective number of different selected TCRs sequences is expected to be much less since extremely low probability sequences may never be formed during the organism’s lifetime. A strict lower bound on Q is the actual number of distinct clones Ĉ in an entire organism [Ĉ ∼ 10⁶ – 10⁸ for humans (, , –34)].

Figure 3

The effects of sampling on two different neutral-model relative clone counts $c_k^s/C^s$ plotted using the dimensionless proliferation rate r = 1/2 in Eqs. 11 and 12 or Eq. 10 and S9 from Section 2 of the Supplementary Material . In (A), we used α = 10^-5, λ =0.01. The effect of sampling is illustrated for η = 1 (no subsampling), 10^-3, 10^-4, 10^-5, and 10^-6. All clone counts are qualitatively similar, with subsampling increasing the exponential decay in c_k . In (B), we use a physiologically unrealistic set of parameters, α = λ = 10, which leads to a qualitatively different unsampled clone count pattern that exhibits a peak. However, under small subsampling fractions η, the clone count loses its peak as it shifts to a rapidly decreasing patterns $c_k^s$ that are not significantly different from sampled clone counts predicted using the parameters α and λ in (A). This indicates inferring parameters using clone counts is ill-conditioned (rather insensitive to parameters) if η is too small.

Figure 4

An exploration of the effects of proliferation rate heterogeneity on the mean clone counts c_k with Q = 10¹³. (A) Various box distributions π_r (r|w) for w = 0, 0.2, 0.4, 0.6, 0.8, and 1. (B) Using Eq. 14 and the dimensionless values $\overline{\alpha }={10}^{-3},\lambda =8\times {10}^{-3}$ such that $\lambda /\overline{\alpha }=8$ , we plot, using the same color spectrum as (A), the corresponding clone counts C_k and show that wider distributions typically generate longer-tailed c_k . However, if λ is set even larger such that $\lambda /\overline{\alpha }=80$ even modest values of w can generate a very long-tailed c_k , as shown in (C). The color spectrum in (C) is for visualization only and not associated with that in (A, B). In the limit of very large $\lambda /\overline{\alpha }$ , the effects of heterogeneous proliferation saturate at very small w beyond which it has negligible effect in further extending the tail.

Figure 5

Ordered integer-valued frequencies j, plotted on a log-log scale, of the C_* distinct (A) alpha and (B) beta chains drawn using OLGA. The index 1 ≤ i ≤ C_* < N_* labels the distinct sequences drawn while b_j is defined as the number these sequences that exhibit the specific frequency j [b ₁ and b ₂ are explicitly indicated in (B)]. The highest frequency clone appears J times such that b_j>J = 0. Since C_* is comparable to N_* , the drawn sequences are dominated by the low probability ones that appear only once. The insets display the frequencies on a linear scale and indicate the long-tailed behavior of the frequencies. The shape of the frequency spectra is self-similar once N_* ≳ 10⁷, allowing us to use this sampling procedure to reliably estimate ${\pi }_{\alpha }\left(\alpha |\overline{\alpha }\right)$ .

Figure 6

The error $H\left(\overline{\alpha },\lambda ,w,\eta \right)$ plotted as a function of α (on a log₁₀ scale) and λ. Darker colors represent smaller values of error as shown by the scale bar on the right. The data used are the clone counts of beta chain sequences of naive CD4 cells from one patient, averaged over three samples. Panels (A–C) use the simple neutral model (Eqs. S9 and S10) and sampling fractions η = 10^-4, 10^-5, and 10^-6, respectively. Since $\overline{\alpha }$ is on a log scale, the error is minimal along a line ${\lambda }_{\min }\propto \overline{\alpha }$ ; the error does not change appreciably along this path and only slightly decreases as λ and $\overline{\alpha }$ become smaller. For the neutral model (w = 0), the error is very sensitive to the sampling fraction η. Here, a fixed, physiologically reasonable value of $\overline{\alpha }$ results in a minimizing λ _min that is unreasonably large, in excess of one and that does not agree well with our expectations of $\lambda =N^{\ast }/Q\ll 1$ . Panels (D–F) show results for the distributed proliferation rate model at full width (w = 1). In this case, the errors are insensitive to the specific choice of η and the minimizing λ _min values are much smaller, consistent with our estimates of N ^∗ and repertoire size. For w = 1, the values of the errors H are also smaller along the ${\lambda }_{min}-\overline{\alpha }$ minimum valley.

Figure 7

Log-log plots of λ _min values as functions of α for η = 10^-4, 10^-5, 10^-5, and 10^-6 for (A) the neutral model, w = 0, and (B) the full-width distributed proliferation rate model, w = 1. These curves trace the values of λ _min along the minimum valley in and show the relative insensitivity of the distributed proliferation rate model to the subsampling fraction η. In both (A, B), the minimum line slopes are near one, with (B) showing a slightly greater slope, indicating λ _min is approximately proportional to $\overline{\alpha }$ over the entire range of w. The color intensity along the lines in (A, B) indicates variation in the total error along the minimum valley; their uniformity shows that the errors are nearly constant along each line. (C) Log-linear plot of λ _min as a function of proliferation rate heterogeneity w for $\overline{\alpha }=2\times 10^{-5},10^{-4}$ . The lower darker curves in each pair correspond to η = 10^-4 while the lighter curves correspond to η = 10^-6. The curves show that even a small heterogeneity w quickly reduces λ _min to below one; however, if λ is forced to be even smaller, the required heterogeneity w increases.

Figure 8

The error $H\left(\overline{\alpha },\lambda ,w,\eta \right)$ using CD4 alpha data from Oakes et al. (12) plotted as a function of w for various $\lambda /\overline{\alpha }$ . We fixed λ = 10^-3 and varied, from left to right, $\overline{\alpha }=2\times {10}^{-5}$ (red), 6 × 10^-5 (green), 10^-4 (blue) and 1.4 × 10^-4 (black). From (A–C), η = 10^-4, 10^-5, and 10^-6. Smaller values of $\lambda /\overline{\alpha }$ result in larger best-fit values of w.

Figure 9

Plots of the representative optimal solutions of clone counts $f_k^s$ from Eq. 19 (using η = 10^-4 and λ = 10^-3 unless otherwise indicated) plotted along side the shown data from Oakes et al. (12). The model predictions and CD4 beta chain data are shown in both (A) log-log and (B) linear scales (there are no zero-values clone counts in this dataset). In (A), the best fit model for the neutral model (w = 0 and ${\pi }_{\alpha }\left(\alpha |\overline{\alpha }\right)=\delta \left(\alpha -\overline{\alpha }\right)$ ) using $\overline{\alpha }=10^{-4}$ is given by λ ≈ 3 shown by the solid black curve. The dashed curves represents best-fit curves using the values associated with the error minima in, where $\overline{\alpha }=2\times 10^{-5}$ , w ≈ 0.09 (red), 6 × 10^-5, w ≈ 0.3 (green), 10^-4, w ≈ 0.53 (blue) and 1.4 × 10^-4, w ≈ 0.76 (black). Note that the neutral model fits well for only the first 2-3 k-points, while the heterogeneous model (w > 0) fits better at larger k.