How heterogeneous thymic output and homeostatic proliferation shape naive T cell receptor clone abundance distributions

TL;DR

This study models how heterogeneity in thymic output and T cell proliferation influences naive T cell receptor clone abundance distributions, revealing that proliferation heterogeneity is key to matching observed data.

Contribution

It introduces a mean-field birth-death-immigration model incorporating TCR-dependent heterogeneity, highlighting proliferation variability as crucial for clone abundance patterns.

Findings

Heterogeneity in proliferation rates significantly impacts clone abundance distributions.

Immigration rate heterogeneity has minimal effect on predicted clone counts.

Model constraints align with experimental clone abundance data.

Abstract

The set of T cells that express the same T cell receptor (TCR) sequence represents a T cell clone. The number of different naive T cell clones in an organism reflects the number of different T cell receptors (TCRs) arising 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 is an important factor in immune function. 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. Using a mean-field approximation to the solution of a regulated birth-death-immigration model and a modification arising from sampling, we systematically quantify 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 comparing predicted clone counts derived from our heterogeneous birth-death-immigration model with experimentally sampled clone abundances, we show that although heterogeneity in immigration rates causes very little change to predicted clone-counts, significant heterogeneity in proliferation rates is necessary to generate the observed abundances with reasonable physiological parameter values. 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.

Equations

$$C=\sum _{k\ge 1}{\hat{c}}_k=Q-{\hat{c}}_0\text{and}N=\sum _{k\ge 1}k{\hat{c}}_k,$$

$$\frac{\text{d}{n}_{\alpha ,r}(t)}{\text{d}t}=\alpha +r{n}_{\alpha ,r}(t)-\mu (N(t))n_{\alpha ,r}(t).$$

$$N(t)=Q{\int }_0^{\infty }\text{d}\alpha {\int }_0^R\text{d}r{n}_{\alpha ,r}(t)\pi (\alpha ,r).$$

$$n_{\alpha ,r}^{\ast }=\frac{\alpha }{\mu (N^{\ast })-r}$$

$$N^{\ast }=Q{\int }_0^R\text{d}r{\int }_0^{\infty }\text{d}\alpha \frac{\alpha \pi (\alpha ,r)}{\mu (N^{\ast })-r},$$

$$\frac{N^{\ast }}{\bar{\alpha }Q}\equiv \frac{\lambda }{\bar{\alpha }}={\int }_0^1\text{d}r\frac{{\pi }_r(r)}{{\mu }^{\ast }-r},$$

$$\frac{\text{d}{c}_k(\alpha ,r)}{\text{d}t}=\alpha \left[c_{k-1}(\alpha ,r)-c_k(\alpha ,r)\right]+r\left[(k-1)c_{k-1}(\alpha ,r)-k{c}_k(\alpha ,r)\right]+\mu (N)\left[(k+1)c_{k+1}(\alpha ,r)-k{c}_k(\alpha ,r)\right],$$

$$c_k(\alpha ,r|{\mu }^{\ast })=Q{\left(1-\frac{r}{{\mu }^{\ast }}\right)}^{\alpha /r}{\left(\frac{r}{{\mu }^{\ast }}\right)}^k\frac{1}{k!}\prod _{\ell =0}^{k-1}\left(\frac{\alpha }{r}+\ell \right).$$

$$\mathbb{P}\left[k|\ell \right]\approx \left(\genfrac{}{}{0ex}{}{\ell }{k}\right){\eta }^k{\left(1-\eta \right)}^{\ell -k},k\le \ell .$$

$$c_k^{\mathrm{s}}(\alpha ,r|{\mu }^{\ast },\eta )\approx \sum _{\ell \ge k}{c}_{\ell }(\alpha ,r|{\mu }^{\ast })\mathbb{P}\left[k|\ell \right]=\sum _{\ell \ge k}{c}_{\ell }(\alpha ,r|{\mu }^{\ast })(\genfrac{}{}{0ex}{}{\ell }{k}){\eta }^k(1-\eta {)}^{\ell -k},$$

$$c_k^{\mathrm{s}}(\alpha ,r|{\mu }^{\ast },\eta )=\frac{Q}{k!}{\left(\frac{\eta r/{\mu }^{\ast }}{1-(1-\eta )(r/{\mu }^{\ast })}\right)}^k{\left(\frac{1-r/{\mu }^{\ast }}{1-(1-\eta )(r/{\mu }^{\ast })}\right)}^{\frac{\alpha }{r}}\prod _{j=0}^{k-1}\left(\frac{\alpha }{r}+j\right).$$

$$C^{\mathrm{s}}(\alpha ,r|{\mu }^{\ast },\eta )=\sum _{k=1}^{\infty }{c}_k^{\mathrm{s}}(\alpha ,r|{\mu }^{\ast },\eta )=Q\left[1-{\left(\frac{1-r/{\mu }^{\ast }}{1-(1-\eta )r/{\mu }^{\ast }}\right)}^{\alpha /r}\right].$$

$${\pi }_r(r)=\{\begin{array}{ll}1/w& \text{if }|r-1/2|<w/2\\ 0& \text{otherwise},\end{array}$$

$${\mu }^{\ast }=\frac{\left(\frac{1}{2}+\frac{w}{2}\right)e^{\lambda w/\bar{\alpha }}-\left(\frac{1}{2}-\frac{w}{2}\right)}{e^{\lambda w/\bar{\alpha }}-1},$$

$${\alpha }_j=\frac{\bar{\alpha }{f}_j}{(N_{\star }/C_{\star })},j=1,2,\dots J,$$

$$c_k(\bar{\alpha },\lambda ,w,\eta )$$

$$=\frac{Q}{k!}\sum _{j=1}^J\frac{b_j}{C_{\star }}{\int }_{(1-w)/2}^{(1+w)/2}\frac{\text{d}r}{w}{\left(\frac{\eta r/{\mu }^{\ast }}{1-(1-\eta )r/{\mu }^{\ast }}\right)}^k{\left(\frac{1-r/{\mu }^{\ast }}{1-(1-\eta )r/{\mu }^{\ast }}\right)}^{\frac{{\alpha }_j}{r}}\prod _{j=0}^{k-1}\left(\frac{{\alpha }_j}{r}+j\right),$$

$$f_k^{\mathrm{s}}(\bar{\alpha },\lambda ,w,\eta )\equiv \frac{k{c}_k^{\mathrm{s}}(\bar{\alpha },\lambda ,w,\eta )}{N^{\mathrm{s}}}=\frac{k{c}_k^{\mathrm{s}}(\bar{\alpha },\lambda ,w,\eta )}{{\sum }_{\ell =1}^{\infty }\ell c_{\ell }^{\mathrm{s}}(\bar{\alpha },\lambda ,w,\eta )}=\frac{k{c}_k^{\mathrm{s}}(\bar{\alpha },\lambda ,w,\eta )}{Q\eta \lambda },$$

$$H(\bar{\alpha },\lambda ,w,\eta )=\sum _{k=1}^{\infty }|f_k^{\mathrm{s}}(\mathrm{data})-f_k^{\mathrm{s}}(\bar{\alpha },\lambda ,w,\eta ){|}^2,$$

$$c_k=Q{\left(1-\frac{\bar{r}}{\mu (N^{\ast })}\right)}^{\bar{\alpha }/\bar{r}}{\left(\frac{\bar{r}}{\mu (N^{\ast })}\right)}^k\frac{1}{k!}\prod _{k^{\prime }=0}^{k-1}\left(\frac{\bar{\alpha }}{\bar{r}}+k^{\prime }\right).$$

$$\frac{c_k}{C}=\frac{c_k}{{\sum }_{\ell \ge 1}{c}_{\ell }}$$

$$\mathbb{P}\left[X=k\right]={\left(1-\frac{\bar{r}}{\mu (N^{\ast })}\right)}^{\bar{\alpha }/\bar{r}}{\left(\frac{\bar{r}}{\mu (N^{\ast })}\right)}^k\frac{1}{k!}\prod _{\ell =0}^{k-1}\left(\frac{\bar{\alpha }}{\bar{r}}+\ell \right).$$

$$\mathbb{P}\left[X=k|X>0\right]=\frac{\mathbb{P}\left[X=k\right]}{{\sum }_{\ell \ge 1}\mathbb{P}\left[X=\ell \right]}=\frac{c_k/Q}{{\sum }_{\ell \ge 1}{c}_{\ell }/Q}=\frac{c_k}{C}.$$

$$\mathbb{E}\left[e^{\xi X}|X>0\right]=\frac{\mathbb{E}\left[e^{\xi X}\right]-\mathbb{E}\left[e^{\xi X}|X=0\right]\mathbb{P}\left[X=0\right]}{\mathbb{P}\left[X>0\right]}.$$

$$\mathbb{E}\left[e^{\xi X}|X>0\right]=\frac{{\left(\frac{1-\bar{r}/\mu (N^{\ast })}{1-e^{\xi }\bar{r}/\mu (N^{\ast })}\right)}^{\bar{\alpha }/\bar{r}}-{\left(1-\frac{\bar{r}}{\mu (N^{\ast })}\right)}^{\bar{\alpha }/\bar{r}}}{1-{\left(1-\frac{\bar{r}}{\mu (N^{\ast })}\right)}^{\bar{\alpha }/\bar{r}}}.$$

$$x^{\bar{\alpha }/\bar{r}}$$

$$\mathbb{E}\left[e^{\xi X}|X>0\right]$$

$$=\frac{\log \left(1-e^{\xi }\frac{\bar{r}}{\mu (N^{\ast })}\right)}{\log \left(1-\frac{\bar{r}}{\mu (N^{\ast })}\right)}+o\left(1\right),$$

$$\underset{\bar{\alpha }/\bar{r}\to 0}{\lim }{c}_k=\frac{C}{\log \left(\frac{1}{1-\bar{r}/\mu (N^{\ast })}\right)}\frac{1}{k}{\left(\frac{\bar{r}}{\mu (N^{\ast })}\right)}^k.$$

$$\frac{\bar{r}}{{\mu }^{\ast }}\to \frac{\lambda }{\lambda +2\bar{\alpha }},\frac{\bar{\alpha }}{\bar{r}}\to 2\bar{\alpha },$$

$$c_k^{\mathrm{s}}(\bar{\alpha },\lambda ,\eta )\equiv \frac{2\bar{\alpha }Q}{k!}{\left(\frac{\eta \lambda }{\eta \lambda +2\bar{\alpha }}\right)}^k{\left(\frac{2\bar{\alpha }}{\eta \lambda +2\bar{\alpha }}\right)}^{2\bar{\alpha }}\prod _{j=1}^{k-1}(2\bar{\alpha }+j).$$

Full text

How heterogeneous thymic output and homeostatic proliferation shape naive T cell receptor clone abundance distributions

Renaud Dessalles1, Yunbei Pan1,2, Mingtao Xia3, Davide Maestrini1 Maria R. D’Orsogna1,2, Tom Chou1,3,

1 Dept. of Computational Medicine, UCLA, Los Angeles, CA 90095-1766

2 Dept. of Mathematics, CalState-Northridge, Los Angeles, CA 91330

3 Dept. of Mathematics, UCLA, Los Angeles, CA 90095-1555 $\ast$ E-mail: tomchou@ucla.edu

Abstract

The set of T cells that express the same T cell receptor (TCR) sequence represents a T cell clone. The number of different naive T cell clones in an organism reflects the number of different T cell receptors (TCRs) arising 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 is an important factor in immune function. 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. Using a mean-field approximation to the solution of a regulated birth-death-immigration model and a modification arising from sampling, we systematically quantify 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 comparing predicted clone counts derived from our heterogeneous birth-death-immigration model with experimentally sampled clone abundances, we show that although heterogeneity in immigration rates causes very little change to predicted clone-counts, significant heterogeneity in proliferation rates is necessary to generate the observed abundances with reasonable physiological parameter values. 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.

Author Summary

It is not known how naive T cell receptor (TCR)-dependent thymic output, proliferation, and death rates control the overall diversity of T cells in the body. Using a reasonable range of physiological parameters within a birth-death-immigration (BDI) model, we quantitatively show how heterogeneity influences the expected sampled clone abundances. We find that heterogeneity in clone-specific immigration rate only weakly affects sampled clone counts while modest heterogeneity in proliferation rates can dramatically affect clone counts. Fitting our model to recently published data, we impose quantitative constraints on biologically meaningful rate parameters in the model. These findings are consistent with the known proliferation-driven maintenance of T cell population in humans.

Introduction

Naive T cells play a fundamental role in the immune system’s response to pathogens, tumors, and other infectious agents. These cells are produced in the thymus, circulate through the blood, and migrate to the lymph nodes where they may be presented with different antigen proteins from various pathogens. Naive T cells mature in the thymus where the so-called V, D, and J segments of genes that code T cell receptors undergo rearrangement. Most T cell receptors (TCRs) are comprised of an alpha chain and a beta chain that are formed after VJ segment and VDJ segment recombination, respectively. The number of possible TCR gene sequences is extremely large, but while recombination is a nearly random process, not all TCRs are formed with the same probability. Before export to the periphery, T cells undergo a selection process, during which T cells with TCRs that react to self proteins are eliminated.

The unique receptors expressed on the cell surface of circulating TCRs enable them to recognize specific antigens; well-known examples include the naive forms of helper T cells (CD4+) and cytotoxic T cells (CD8+). The set of naive T cells that express the same TCR are said to belong to the same T cell clone. Upon encountering the antigens that activate their TCRs, naive T cells turn into effector cells that assist in eliminating infected cells. Effector cells die after pathogen clearance, but some develop into memory T cells. Because of the large number of unknown pathogens, TCR clonal diversity is a key factor for mounting an effective immune response. Recent studies also reveal that human TCR clonal diversity is implicated in healthy aging, neonatal immunity, vaccination response and T cell reconstitution following haematopoietic stem cell transplantation [1, 2]. Despite the central role of the naive T cell pool in host defense, and broadly speaking in health and disease, TCR diversity is difficult to quantify. For example, the human body hosts a large repertoire of T cell clones, however the actual distribution of clone sizes is not precisely known [3]. Only recently have experimental and theoretical efforts been devoted to understanding the mechanistic origins of TCR diversity [4, 5, 6, 7, 8, 9]. The goal of this work is to formulate a realistic mathematical model that includes heterogeneous naive T cell generation and reproduction rates and that we will use to describe recent experimental results.

A well-established way to describe the T cell repertoire is by determining the clone abundance distribution or “clone count” $\hat{c}_{k}$ (for $k\geq 1$) that measures the number of distinct clones represented by exactly $k$ T cells: $\hat{c}_{k}\equiv\sum_{i=1}^{\infty}\mathds{1}(n_{i},k)$, where $n_{i}$ is the discrete number of T cells carrying TCR $i$ and the indicator function $\mathds{1}(n,k)=1$ if $n=k$ and [math] otherwise. This distribution captures the entire pattern of the clonal populations. Several summary indices for T cell diversity such as Shannon’s entropy, Simpson’s index, or the whole population richness can be deduced from the distribution $\hat{c}_{k}$ [10]. Note that $\hat{c}_{k}$ counts only the number of clones of a specific population $k$ and does not carry any TCR sequence or identity information.

Complete clone counts $\hat{c}_{k}$ and the total number of circulating naive T cells are difficult to measure in humans. Nonetheless, high-throughput DNA sequencing on samples of peripheral blood containing T cells [11, 12, 13, 14] have provided some insight into TCR diversity. A commonly observed feature is that clone counts $\hat{c}_{k}$ often exhibit a power-law distribution in the clone abundance $k$, albeit for approximately one or at most two decades (see Fig. 1).

Several models have been developed to explain certain features of observed clone counts [3, 15, 16, 17, 4], including the apparent power-law behavior. One proposal is that T cells in different clones have TCRs that have different affinities for self-ligands that are necessary for peripheral proliferation [4, 5, 6], leading to clone specific replication rates. An alternative hypothesis [7] is that specific TCR sequences are more likely to arise in the V(D)J recombination process in the thymus [18] leading to a higher probability that these TCRs are produced. De Greef et al. [7] estimated the probability of production of a given TCR sequence by using the Inference and Generation of Repertoires (IGoR) simulation tool that quantitatively characterizes the statistics of receptor generation from both cDNA and gDNA data [18]. However, the data indicate that power laws, which represent a simple summary of sampled clone counts, arise only across one or two decades of clone sizes, as shown in Fig. 1. Moreover, the above studies have not systematically incorporated and compared heterogeneity in both immigration and replication rates, sampling, or fitting to measured TCR clone abundance distributions.

In this paper, we analyze the effects of heterogeneity and sampling within a mean-field model based on a stochastic multiclone birth-death-immigration (BDI) process. The model allows for TCR-sequence dependent replication and immigration rates. Our model is derived from an underlying continuous-time Markovian birth-death-immigration (BDI) process [10] described by: (i) immigration representing the arrival of new clones from the thymus; (ii) birth during homeostatic proliferation of naive T cells that yield newborn naive T cells with the same TCR as their parent; and (iii) death representing cell apoptosis. We also include a regulation, or “carrying capacity,” mechanism through a total population-dependent death rate which may represent the global competition for cytokines, such as Interleukine-7 [19, 20, 21, 22, 23], needed for naive T cell survival and homeostasis [24, 25]. This interaction will ensure a homeostatic finite naive T cell population and will be assumed to be clone-independent since these cytokine signals are TCR-independent [21].

We derive analytical results of our heterogeneous BDI model that are applicable on the scale of the entire organism. These results are then extended to predict the expected clone counts in a small subsample of the population. Finally, this result is averaged over distributions of TCR immigration (thymic output) and homeostatic proliferation rates. Heterogeneity in TCR production rates are extracted from new computational tools: Inference and Generation of Repertoires (IGoR) [18] and Optimized Likelihood estimate of immunoGlobulin Amino-acid sequences (OLGA) [26]. Since there are no equivalent tools that measures proliferation rates, we will assume simple functional forms for the distribution of homeostatic proliferation rates.

Our calculations provide insight into how parameters describing the shape of the distribution of immigration and proliferation rates affect the shape of the expected clone counts. To compare with experimental measurements, we also quantify the random sampling process that describes actual measurements derived from blood draws. Our model allows us to estimate the values of relevant physiological parameters and reveals that the shape of the sampled expected clone counts are much less sensitive to TCR-dependent thymic output than to TCR-dependent proliferation. However, we show how the width of a simple uniform proliferation rate distribution can inform the observed $\hat{c}_{k}$.

Materials and Methods

To understand the observed clone counts, we focus on the clone abundance distribution $\hat{c}_{k}$ associated only with naive T cells, the first type of cells produced by the thymus that have not yet been activated by any antigen. Antigen-mediated activation initiates a one-way irreversible cascade of differentiation into effector and memory T cells that we can subsume into a death rate. Thus, we limit our analysis to birth, death, and immigration within the naive T cell compartment. Here, we first present the mathematical framework of the BDI process to provide an initial qualitative understanding for clone counts.

Non-neutral Birth-Death-Immigration model

The multiclone BDI process is depicted in Fig. 2.

We define $Q$ to be the theoretical number of possible functional naive T cell receptor clones that can be generated by V(D)J recombination in the thymus which is estimated to be $Q\sim 10^{13}-10^{18}$ [6, 26]. As we will later show, results of our model will not be depend the explicit value of $Q$ as long as $Q\gg 1$. Due to naive T cell death or removal from the sampling-accessible pool, not all possible clone types will be presented in the organism, so we denote the number of clones actually present in the body (or “richness”) by $C\ll Q$, where estimates of $C$ range from $\sim 10^{6}-10^{8}$ in mice and humans [29, 30, 1, 6, 32].

Although naive T cells are difficult to distinguish from the entire T cell population, the total number of naive T cells (across all clones present) in human has been estimated to be about $N\sim 10^{11}$. However, circulating naive T cells number perhaps $10^{9}$ but can exchange, at different time scales, with those that reside in peripheral tissue, which may carry their own proliferation and death rates. The effective pool that is ultimately sampled is thus difficult to estimate, but measurements show the theoretical number of different clones $\gg$ the total number of naive T cells $\gg$ the total number of different T cells clones actually in the body ($Q\gg N\gg C$). Regardless of the precise values of the discrete quantities $Q,C,N$, they are related to the discrete clone counts $\hat{c}_{k}$ through

[TABLE]

where $\hat{c}_{0}$ is the number of possible clones that are not expressed in the organism.

Each distinct clone $i$ (with $1\leq i\leq Q$) is characterized by an immigration rate $\alpha_{i}$ and a per cell replication rate $r_{i}$. The immigration rate $\alpha_{i}$ is clone-specific because it depends on the preferential V(D)J recombination process; the replication rate $r_{i}$ is also clone-specific due to the different interactions with self-peptides that trigger proliferation. Since both the numbers of theoretically possible ($Q\gg 1$) and observed ($C\gg 1$) clones are extremely large, we can define a continuous density $\pi(\alpha,r)$ from which immigration and proliferation rates $\alpha$ and $r$ are drawn. This means that the probability that a given clone has an immigration rate between $\alpha$ and $\alpha+\textrm{d}\alpha$ and replication rate between $r$ and $r+\textrm{d}r$ is $\pi(\alpha,r)\textrm{d}\alpha\textrm{d}r$. Since $Q$ is finite and countable, there are theoretical maximum values $A$ and $R$ for the immigration and proliferation rates, respectively, such that $\pi(\alpha,r)=0$ for $\alpha\geq A$ and $r\geq R$. In the BDI process, the upper bound $R$ on the proliferation rate prevents unbounded numbers of naive T cells and is necessary for a self-consistent solution. Conversely, results are rather insensitive to the upper bound $A$ on the immigration rate provided $\pi(\alpha,r)$ decays sufficiently fast such that the mean value of $\alpha$ is finite for all $r\leq R$. Thus, for simplicity, we henceforth take the $A\to\infty$ limit. The heterogeneity in the immigration and replication rates allows us to go beyond typical “neutral” BDI models, where both rates are fixed to a specific value for all clones.

Finally, we assume the per cell death rate $\mu(N)$ is clone-independent but a function of the total population $N$. This dependence represents the competition among all naive T cells for a common resource (such as cytokines), which effectively imposes a carrying capacity on the population [22, 28, 33]. The specific form of the regulation will not qualitatively affect our findings since we will ultimately be interested in only its value $\mu(N^{*})\equiv\mu^{*}$ at the steady state population $N^{*}$.

Mean-Field Approximation of the BDI Process

The exact steady-state probabilities of configurations of the discrete abundances $\hat{c}_{k}$ for a fully stochastic neutral BDI model with regulated death rate $\mu(N)$ were recently derived [10]. To incorporate TCR-dependent immigration and replication rates in a non-neutral model, we must consider distinct values of $\alpha_{i}$ and $r_{i}$ for each clone $i$. In this case, an analytic solution for the probability distribution over $\hat{c}_{k}$, even at steady state, cannot be expressed in an explicit form. However, since the effective number of naive T cells ($N\sim 10^{9}-10^{11}$ [32]) is large, we can exploit a mean-field approximation to the non-neutral BDI model and derive expressions for the mean values of the discrete clone counts $\hat{c}_{k}$. We will show later that under realistic parameter regimes, the mean-field approximation is quantitatively accurate. Breakdown of the mean field approximation has been carefully analyzed in other studies [34].

Deterministic approximation for total population and effective death rate:

To implement the mean-field approximation in the presence of a general regulated death rate $\mu(N)$, we start by writing the deterministic, “mass-action” ODE for the mean density $n_{\alpha,r}(t)$ of cells having immigration rate near $\alpha$ and proliferation rate near $r$ in a BDI process

[TABLE]

Since $Q\gg 1$, the total number of clones that carry an immigration rate between $\alpha$ and $\alpha+\textrm{d}\alpha$ and a replication rate between $r$ and $r+\textrm{d}r$ is approximately $Q\pi(\alpha,r)\textrm{d}\alpha\textrm{d}r$. The total mean number $N(t)$ of naive T cells can then be estimated as a weighted integral over all $n_{\alpha,r}(t)$

[TABLE]

At steady-state, the solution to Eq. 2 can be simply expressed as

[TABLE]

in which $N^{*}$ is the steady-state value of $N(t)$. Thus, upon averaging Eq. 4 over $\alpha$ and $r$, we find

[TABLE]

a self-consistent equation for $N^{*}$ if the functional form of $\mu(N^{*})$ is known. Without the finite upper bound $R$ of the density $\pi(\alpha,r)$, the integral in Eq. 5 diverges. The self-consistent determination of $N^{*}$ depends implicitly on the parameters that define the distribution $\pi(\alpha,r)$.

Note that for factorisable $\pi(\alpha,r)=\pi_{\alpha}(\alpha)\pi_{r}(r)$, the self-consistent death rate $\mu^{*}$ depends only on the combination $N^{*}/(\bar{\alpha}Q)=\int_{0}^{R}\!\textrm{d}r\,\pi_{r}(r)/(\mu^{*}-r)$, where $\bar{\alpha}\equiv\int_{0}^{\infty}\alpha\pi_{\alpha}(\alpha)$ and $\mu^{*}\equiv\mu(N^{*})$ is the effective death rate at steady-state. We now nondimensionalize all rates with respect to the maximum proliferation rate $R=2\bar{r}<\mu^{*}$. In these units, $\bar{r}=1/2$ and $0<w\leq 1$. Henceforth, we will use the notation $\bar{\alpha}/R\to\bar{\alpha}$ and $\mu^{*}/R\to\mu^{*}$ so that the quantities $\bar{\alpha}$ and $\mu^{*}$ are meant to be dimensionless, unless otherwise stated. The steady-state self-consistent condition becomes

[TABLE]

where $\lambda\equiv N^{*}/Q$ is the dimensionless ratio of the total naive T cell population to the total effective number of clones $Q$. We expect $\lambda\ll 1$. For a given $\pi_{r}(r)$, one can then self-consistently determine $\mu^{*}$ in terms of $\lambda/\bar{\alpha}$ and the parameters defining $\pi_{r}(r)$.

Mean-field model of clone counts:

We now use the results for the mean steady-state population $N^{*}$ to find the clone counts averaged over all realizations of the underlying stochastic process. The mean-field equations for the dynamics of these mean clone counts in the neutral model were derived in [35, 34] and are reviewed in Appendix 1. In the neutral model, we assume that all effective clones $Q$ are associated with the same rates $\alpha$ and $r$ so that the mean field evolution equation for $c_{k}(\alpha,r)$ is [35, 34]

[TABLE]

along with the constraint $\sum_{k=0}^{\infty}c_{k}(\alpha,r)=Q$. This equation assumes that both $c_{k}(\alpha,r)$ and $N$ are uncorrelated, allowing us to write the last term as a product of functions of the mean population $N=\sum_{k}kc_{k}$ and $c_{k+1},c_{k}$.

At steady state we replace $\mu(N)$ with $\mu(N^{*})\equiv\mu^{*}$ self-consistently found by solving Eq. 6. The solution follows a negative binomial distribution with parameters $\alpha/r$ and $r/\mu^{*}$ [10, 34]

[TABLE]

The dependence on $\lambda/\bar{\alpha}$ arises from the determination of $\mu^{*}(\lambda/\bar{\alpha})$ through Eq. 6. Although $c_{k}(\alpha,r|\mu^{*})$ has not yet been averaged over $\alpha,r$, it also implicitly depends on parameters that define $\pi_{r}(r)$ through the solution for $\mu^{*}$ from Eq. 6.

Sampling

Unless an animal is sacked and its entire naive T cell population is sequenced, TCR clone distributions are typically measured from sequencing TCRs in a small blood sample. In such samples, low population clones may be missed. In order to compare our predictions with measured clone abundance distributions, we must revise our predictions to allow for random cell sampling.

We define $\eta$ as the fraction of naive T cells in an organism that is drawn in a sample and assume that all naive T cells in the organism have the same probability $\eta$ of being sampled. This is true only if naive T cells carrying different TCRs are not preferentially partitioned into different tissues and are uniformly distributed within an animal. Let us assume that a specific TCR is represented by $\ell$ cells in an organism. If $N^{*}\eta\gg\ell$, the probability that $k$ cells from the same clone are sampled approximately follows a binomial distribution with parameters $\ell$ and $\eta$ [36, 37, 38, 39, 40]

[TABLE]

The associated mean sampled clone count $c^{\rm s}_{k}$ depends on the clone count $c_{\ell}(\alpha,r|\mu^{*})$ predicted in the body via the formula

[TABLE]

where $c_{\ell}(\alpha,r|\mu^{*})$ is determined by Eq. 8. As mentioned, the clone counts $c_{\ell}(\alpha,r|\mu^{*})$ and $c_{k}^{\rm s}(\alpha,r|\mu^{*},\eta)$, before averaging over $\pi(\alpha,r)$, depend implicitly on $\lambda/\bar{\alpha}$ and other parameters that define $\pi_{r}(r)$ through the determination of $\mu^{*}$ via Eq. 6. The sum in Eq. 10 can be explicitly performed to yield

[TABLE]

The total expected number of clones predicted in the sample comprised of cells within $[\alpha,\alpha+\textrm{d}\alpha]$ and $[r,r+\textrm{d}r]$ can be found from direct summation:

[TABLE]

The effects of sampling are shown in Fig. 3. Subsampling greatly affects the observed clone counts, increasing the exponential decay of $c_{k}$ [38]. This shifts $c_{k}$ at large $k$ to much smaller values of $k$, while reducing the values of $c_{k}$ for small $k$. In Fig. 3(a) and (b) we use two very different parameter values $\alpha=10^{-5},r=1/2,\lambda=0.01$ and $\alpha=\lambda=10,r=1/2$ to generate two qualitatively very different $c_{k}$. If the subsampling $\eta\ll 1$ is sufficiently small, the resulting $c_{k}^{\rm s}$ corresponding to qualitatively different $c_{k}$ can become qualitatively similar. This implies that sampling may make the estimation of whole-body clone counts from sampled data somewhat ill-conditioned, i.e., different whole-body clone counts, upon sampling, may yield similar sampled clone counts.

Although sampling has a dominating effect on the inference of $c_{k}$, immigration and proliferation rate distributions may also strongly affect $c_{k}$. By averaging the clone counts $c_{k}(\alpha,r|\mu^{*})$ (Eq. 8), the sampled clone counts $c_{k}^{\rm s}(\alpha,r|\lambda,\eta)$ (Eq. 11), and the total clone count $C^{\rm s}(\alpha,r|\mu^{*},\eta)$ (Eq. 12) over the distribution $\pi(\alpha,r)$, we can make predictions on expected clone counts in the body and in the sample.

Heterogeneity and determination of $\pi(\alpha,r)$

The fundamental result given in Eq. 11 applies only to the clone count density (expected clone counts with immigration rate within $[\alpha,\alpha+\textrm{d}\alpha]$ and proliferation rates within $[r,r+\textrm{d}r]$). Now, we propose forms for $\pi(\alpha,r)$ in order to evaluate $c_{k}^{\rm s}(\mu^{*},\eta)=\int_{0}^{\infty}\!\!\textrm{d}\alpha\int_{0}^{1}\!\!\textrm{d}r\,\pi(\alpha,r)c_{k}^{\rm s}(\alpha,r|\mu^{*},\eta)$. For simplicity, we assume that $\alpha$ and $r$ are uncorrelated and that the distribution factorises: $\pi(\alpha,r)=\pi_{\alpha}(\alpha)\pi_{r}(r)$.

Proliferation rate heterogeneity:

First, we consider a distribution of TCR sequence-dependent proliferation rates. Since TCR-antigen affinity depends on the receptor amino-acid sequence, the rate of T cell activation and subsequent proliferation is clone-specific [41, 28]. Thus, the specific interactions between TCRs and low-affinity MHC/self-peptide complexes maps to a distribution of proliferation rates among all the $Q$ possible clones. Since there are no data (known to us) that can be used to infer this mapping or the specific shape of $\pi_{r}(r)$, we assume, for simplicity, a simple uniform “box” distribution centered about a mean value $\bar{r}=1/2$:

[TABLE]

where $0\leq w\leq 1$ represents the relative width of the uniform box distribution. The dimensionless self-consistency condition (Eq. 6) thus yields

[TABLE]

where $\bar{\alpha}$ is the mean per-clone immigration rate normalized by the maximum proliferation rate $R$.

Before we use these formulae to understand the experimentally measured clone counts, we explore the effects of proliferation rate heterogeneity on the clone counts. In Fig. 4, we plot the predicted, whole-organism ($\eta=1$) clone counts for different values of $w$. Since the function $c_{k}$ has an exponentially decaying term $(r/\mu^{*})^{k}$, conditions that lead to larger $\mu^{*}$ lead to the slowest decaying $c_{k}$.

A fixed normalized value of $r=1/2$ leads to an exponential decay in $c_{k}$ of approximately $2^{-k}$. However, if $r$ is integrated over its full width of values $(0,1)$, the contribution from $r\approx 1$ generates a longer-tailed distribution $c_{k}$. Note that for large $\lambda/\bar{\alpha}\gtrsim 1$, smaller $w$ gives rise to smaller $\mu^{*}$, and even a modest $w$ is enough to make the $c_{k}$ decay slowly.

Immigration rate heterogeneity:

Next, we propose a way of evaluating formulating a distribution $\pi_{\alpha}(\alpha)$ that represents the TCR sequence-dependent output from the thymus. A specific form for $\pi_{\alpha}(\alpha)$ can be obtained from previous studies that predict V(D)J recombination frequencies associated with each TCR sequence. The statistical model for differential V(D)J recombination in humans is implemented in the Optimized Likelihood estimate of immunoGlobulin Amino-acid sequences (OLGA) software [26], which is an updated version of the Inference and Generation of Repertoires (IGoR) software [18]. Below, we estimate $\pi_{\alpha}(\alpha)$ by sampling a large number of TCRs from OLGA that draws sequences according to their generation probability. Our working assumption is that thymic selection is uncorrelated with V(D)J recombination so the relative probabilities of forming different TCRs provide an accurate representation of the ratios of the TCRs exported into the periphery.

Both IGoR and OLGA can be used to generate the probabilities corresponding to each drawn sequence but this requires significant computational time and memory. However, since the sequence draws are proportional to the underlying probabilities, we simply drew $N_{\star}=10^{9}$ sequences and exploited the probability-biased selection of sequences by counting the frequencies of each sequence. Out of $N_{\star}$ sequence draws from IGoR or OLGA, there will be $C_{\star}$ distinct sequences with $b_{j}$ of them occurring $j$ times.

For the alpha and beta chains, we found the frequencies $f_{j}$, $j=1,2,3,\ldots,J$, where $J$ is the total number of different integer-valued frequencies observed. The frequencies of the drawn clones are plotted in Figs. 5(a) and (b) in decreasing order. For the alpha chain, we found $C_{\star}=372,806,648\approx 3.72\times 10^{8},J=52,294$ while for the beta chain we found $C_{\star}=875,920,705\approx 8.76\times 10^{8},J=6430$. The probability that a sequence has integer-valued frequency $f_{j}$ is $b_{j}/C_{\star}$. From these numbers, we can effectively average over $\pi_{\alpha}(\alpha)$ by taking a sum over discrete values of $\alpha$: $\int\pi_{\alpha}(\alpha)c_{k}^{\rm s}(\alpha)\textrm{d}\alpha\to\sum_{j=1}^{J}{b_{j}\over C_{\star}}c_{k}^{\rm s}(\alpha_{j})$ with

[TABLE]

where here, $\bar{\alpha}$ is the dimensionless per-clone immigration rate, averaged over the set of $C_{\star}$ unique sequences.

Now that we have defined the possible distributions for $\pi_{\alpha}(\alpha)$ and $\pi_{r}(r)$, we can compute the mean, sampled, immigration- and proliferation-averaged clone counts and compare them with measurements. The full formula for the clone counts is thus

[TABLE]

where $\alpha_{j}(\bar{\alpha})$ is given by Eq. 15 and $\mu^{*}$ is given by Eq. 14. Eq. 16 is thus our “full model” from which we make predictions of clones count-related quantities and compare them with data. Using this equation, we can mathematically study the importance of the heterogeneities in $\alpha$ and $r$ by comparing predictions from simpler forms of $\pi_{\alpha}(\alpha),\pi_{r}(r)$. In Appendix 2 in the SI, we provide explicit dimensionless formulae for different combinations of $\pi_{\alpha}(\alpha)$ and $\pi_{r}(r)$. For example, we can consider fixed $\alpha$ or $r$ by using $\pi_{\alpha}(\alpha)=\delta(\alpha-\bar{\alpha})$ or $\pi_{r}(r)=\delta(r-{1\over 2})$, respectively.

The effective distribution $b_{j}$ over different values $\alpha_{j}$ indicates that most drawn clones only appear once since $b_{1}\gg b_{j\geq 2}$. Thus, we hypothesize that for sufficiently small $\bar{\alpha}$, our formulae for $c_{k}$ and all subsequent quantities can be simplified by taking the limit $\bar{\alpha}/r\ll 1$. Indeed, as we show in Appendix 3, such a simpler expression remains highly accurate provided the dimensionless $\bar{\alpha}<10^{-2}$ and allows efficient computation. Moreover, this implies that, mathematically, the result arising from weighting $c_{k}(\alpha,r)$ over $\pi_{\alpha}(\alpha)$ can be approximated by a single effective value $\bar{\alpha}$, reinforcing our overall conclusion that proliferation rate heterogeneity is the dominant factor. Nonetheless, we will use the full summation over $\alpha_{j}$ (Eq. 16) in our subsequent analyses.

Results and Analysis

Before performing a quantitative comparison with measured clone counts from Oakes et al. [12], we discuss the qualitative features of our model and typical physiological parameter ranges. While even the basic model parameters are difficult to measure and vary widely across different organisms and individuals, our nondimensionalized model unifies the mechanisms and concepts common to the maintenance of diversity in the T cell repertoire across different organisms.

When considering the data, we observe that even after significant subsampling, there are appreciable clone counts at reasonably large clone sizes $k$, whereas the unsampled clone counts decay exponentially in $k$ with rate $\log(\mu^{*}/r)$. Even though $r$ may take on a range of values, as determined by $\pi_{r}(r)$, the slowest decay of $c_{k}$ in $k$ arises from the largest possible values of $r$. Thus, a larger proliferation rate heterogeneity $w$ will generally yield a longer-tailed $c_{k}$, as illustrated in Fig. 4.

Since the data we will analyze are derived from human samples, we will use the following arguments as a guide to roughly estimate the relevant range of parameters:

• •

The average total number of statistically relevant circulating naive T cells is not completely known but estimates are $N^{*}\sim 10^{11}$ [32]. However, the circulating population in the peripheral blood is approximately two orders of magnitude smaller. These circulating naive T cells nonetheless exchange with those in the much larger population in the lymph and other tissues. The timescale of this exchange (relative to the age of the organism being sampled or the intersample times) will determine the effective $N^{*}$ relevant for sampling clone counts $c_{k}^{\rm s}$. Thus, we will use an order-of-magnitude estimate on the lower range of measurements and estimate $N^{*}\sim 10^{10}-10^{11}$.

• •

The total possible number $Q$ of TCRs of either alpha or beta chains may be in the range $10^{13}-10^{18}$ [42], but the effective number is probably much smaller as many clones may never be generated in a lifetime. Thus, an effective value of $Q$ may reside at the lower range which leads to $\lambda\equiv N^{*}/Q\sim 10^{-4}-10^{-2}$.

• •

The average (dimensional) immigration rate per clone $\bar{\alpha}$ can be deduced from the total thymic output of all clones $\bar{\alpha}Q$, which has been estimated across a wide range of values $\bar{\alpha}Q\sim 10^{7}-10^{8}$/day [43, 44, 45, 46, 47]. If we use an effective repertoire size of $Q\sim 10^{13}-10^{14}$, the average per clone immigration rate becomes $\bar{\alpha}\sim 10^{-7}-10^{-5}/\text{day}$.

• •

The mean proliferation rate $\bar{r}$ is difficult to measure but has been estimated to be on the order of $\bar{r}\sim 10^{-4}-10^{-3}/\text{day}$ [47]. If we nondimensionalize using $2\bar{r}$, the dimensionless $\bar{\alpha}\sim 10^{-4}-10^{-1}$.

• •

The sampling fraction $\eta$, although in principle determined experimentally, is also hard to quantify due to the uncertainty in $N^{*}$. Blood sampling volume fractions from humans are typically $\eta\sim 10^{-3}$; however, in recent experiments [12] the number of enumerated sequences $\sim 10^{5}$, which, given rough estimates of effective $N^{*}\sim 10^{10}-10^{11}$, yield $\eta\sim 10^{-6}-10^{-4}$. Due to this uncertainty in $\eta$, we will explore different fixed values of $\eta$ around $10^{-5}$.

Using the above guide for reasonable parameter ranges, we now consider fitting our results in Eqs. 16, S9, S10, S11, S12, S13, and S14 to some of the available data [12]. Before doing so, note that although the log-log plots shown in Figs. 1(a,b) provide a simple visual for $\log c_{k}^{\rm s}$ or $\log[c_{k}^{\rm s}/C^{\rm s}]$, fitting must be performed on the linear scale.

The measured data includes data at values of $k$ for which no clones were detected so that $c_{k}^{\rm s}=0$. On the log scale these data points include $\log c_{k}^{\rm s}\to-\infty$ which nonetheless should be included in the fitting as they represent realizations of the system. Numerical fitting on the log-log scale would be misleading once a value of $c_{k}^{\rm s}=0$ is encountered. Thus, will fit our mean-field model on the linear scale to the fraction $f_{k}^{\rm s}$ of the total number of sampled cells that are in clones of size $k$:

[TABLE]

where the denominator $Q\eta\lambda$ comes directly from the definition $\sum_{\ell=1}^{\infty}\ell c_{\ell}^{\rm s}(\bar{\alpha},\lambda|\eta)\equiv N^{\rm s}$, the sampling relation $N^{\rm s}=\eta N^{*}$, and Eq. 6. Rather than using $N^{\rm s}$ directly from the number of reads in an experimental sample, equivalently, we use the model expression $N^{\rm s}=Q\eta\lambda$ to arrive at the last equality in Eq. 17. This form ensures strict normalization $\sum_{k=1}^{\infty}f_{k}^{\rm s}=1$ and is independent of the unknown repertoire size $Q$ ($c_{k}^{\rm s}$ is proportional to $Q$). Note that for $f_{k}^{\rm s}\equiv kc_{k}^{\rm s}/N^{\rm s}$, the implicit factor of $Q$ in $c_{k}^{\rm s}$ (Eq. 11) can cells the explicit $Q$ in the denominator of Eq. 17. Thus, $f_{k}^{\rm s}$ as well as $c_{k}^{\rm s}/C^{\rm s}$ depend on the effective value of $Q$ only through the determination of $\mu^{*}$ through $\lambda\equiv N^{*}/Q$ in Eq. 6.

Our mathematical framework provides only mean sampled clone counts while each sample of the data represents one realization. Large sample-to-sample variations in the clone counts would render the fitting less informative, but these large variations were not seen in the triplicate samples in Oakes et al. [12]. Mechanistically, we expect that for large $k$ the number of cells contributing to $f_{k}^{\rm s}$ is also large so demographic stochasticity is relatively small and results in small uncertainties in the value of $k$, and not in the magnitude of $f_{k}^{\rm s}$. Large clones are also likely include memory T cells that have been produced after antigen stimulation of specific clones. Memory T cells are difficult to accurately distinguish from naive T cells [12] but we will see that large $k$ components of $f_{k}^{\rm s}$ negligibly influence the fitting.

For small $k$, the numbers of clones are large so we expect sampling and intrinsic stochasticity to also be controlled. Thus, we averaged the clone counts from three different blood samples from each patient in Oakes et al. [12] and only compared them to the expected value of the sampled clone counts in our model.

By comparing our model $f_{k}^{\rm s}(\bar{\alpha},\lambda,w,\eta)$ with the data $f_{k}^{\rm s}({\rm data})$ via the total least-squares error

[TABLE]

we explore how $H(\bar{\alpha},\lambda,w,\eta)$ depends on the parameters $\bar{\alpha}$, $\lambda$, $w$, and sampling fraction $\eta$.

Since the variation was not large and there were already too many uncertainties (particularly the sampling fraction $\eta$, the shape of $\pi_{r}(r)$, and the accuracy of a steady-state approximation), we did not pursue a much more technical uncertainty analysis. Our goal was to simply consider the data using a reasonable initial model that captures the known general features of naive T cell maintenance (immigration, birth, death) and find constraints among the parameters $\lambda,\bar{\alpha}$, and $w$.

In Figs. 6(a-c) the data $f_{k}^{\rm s}(\rm data)$ were derived from the average of three measurements from beta chain CD4 sequences sampled from one patient [12]. Using this data, we compute and plot the error $H(\bar{\alpha},\lambda,w=0,\eta=10^{-4},10^{-5},10^{-6}))$ as a function of $\lambda$ for various values of $\bar{\alpha}$ using the neutral model ($w=0$, Eq. S9 in Appendix 2). For reasonable values of $\bar{\alpha}$ and sampling fraction $\eta$, we see that the predicted $\lambda\equiv N^{*}/Q$ is typically ${\cal O}(1)$ or larger. In Figs. 6(d-f) we use the full-width ($w=1$) distribution of $\pi_{r}(r)$ to show the analogous error for the same data using the same sampling fractions $\eta=10^{-4},10^{-5},10^{-6}$. Note that the optimal (low error) values of $\lambda$ are significantly smaller that those in found using $w=0$ (Figs. 6(a-c)) and that the results are rather insensitive to sampling fraction $\eta$. These smaller values of $\lambda=N^{*}/Q$ are more consistent with known physiological understanding. Thus, the distributed proliferation rate model provides a much more self-consistent fit to the data than the fixed proliferation rate neutral model.

Note that the value of the errors along the minimum valley are nearly constant, only slightly decreasing as $\lambda,\bar{\alpha}\to 0$. This means that the model and associated data constrains $\lambda$ to $\bar{\alpha}$, but cannot independently determine them. The constraint arises as a linear relationship as shown in Fig. 7 that shows the relationship between $\lambda$ and $\bar{\alpha}$ along the minimizer of $H(\bar{\alpha},\lambda,\eta)$.

Figs. 7 show the optimal $\lambda$ as a function of $\bar{\alpha}$ for (a) the neutral model and (b) the full-width $\pi_{r}(r)$ model. Although the variation in error is negligible across $\bar{\alpha}$ for each model $w$ (color variation along the lines), the full-width distributed model ($w=1$) generates a larger error than the neutral model, but it provides much more reasonable values of $\lambda,\bar{\alpha}$ and is much less sensitive to $\eta$. To investigate intermediate values of $0<w<1$ to determine if a low error can be achieved at small values of $\lambda\ll 1$, we plot in Fig. 7(c) the optimal $\lambda$ as a function of proliferation rate heterogeneity $w$. Note that even a small heterogeneity significantly reduces the predicted $\lambda$. However, if $\lambda\approx 10^{-3}$ or $10^{4}$ say, the required $w$ can become quite large.

To explore the dependence of the error on the proliferation rate heterogeneity, we fix $\bar{\alpha},\lambda$, and $\eta$, and plot the error $H(\bar{\alpha},\lambda,w,\eta)$ as a function of $w$.

Fig. 8 shows how the minimizing value of $w$ is sensitive to $\lambda/\bar{\alpha}$, suggesting larger proliferation heterogeneity as $\lambda/\bar{\alpha}$ is decreased. The minimum of the error, however, is rather insensitive to $\lambda/\bar{\alpha}$. For sufficiently small $\lambda/\bar{\alpha}$, the minimum error at $w=1$ is higher than those with minimum $w<1$. Since the minima in the interior $0<w<1$ are nearly invariant to $\lambda/\bar{\alpha}$, we indeed show that near-optimal solutions with $\lambda\ll 1$ can be found when the proliferation rate heterogeneity $w$ is appreciable.

Using the parameters associated with the minima in Fig. 8, we plot our predicted $f_{k}^{\rm s}$ against the data $f_{k}^{\rm s}({\rm data})$ in Fig. 9.

Although the neutral model appears to fit better, especially for small $k$, the necessary values of $\lambda$ and $\bar{\alpha}$ are too large and small, respectively. Conversely, when proliferation rate heterogeneity is allowed, the best-fits have small error and are found using $\lambda\ll 1$. Note that most of the information in the data lies in how $f_{k}^{\rm s}({\rm data})$ decreases over the first few values of $k$. The goodness of fit of our model to the data depends mostly on the predicted initial decreases in $f_{k}^{\rm s}(\bar{\alpha},\lambda,w,\eta)$.

Discussion

Here, we review and justify a number of critical biological assumptions and mathematical approximations used in our analysis. The effects of relaxing our approximations are also discussed.

Factorisation of $\pi(\alpha,r)$.

For mathematical tractability, we have assumed $\pi(\alpha,r)=\pi_{\alpha}(\alpha)\pi_{r}(r)$. Given the typical physiological values of $\bar{\alpha}$, the clone count formulae derived from our model can be accurately approximated by a single value of $\bar{\alpha}$. Thus, we expect that the immigration rate distribution can be approximated by $\pi_{\alpha}(\alpha)=\delta(\alpha-\bar{\alpha})$. This allows further approximation of our formulae as shown in Appendix 3. In Appendix 4, we explicitly show that factorisation is an accurate approximation.

We have also assumed that selection is uncorrelated with the generation probabilities of the TCR nucleotide sequences encoded in IGoR/OLGA. The assumption is that the recombination statistics are uncorrelated with the statistics of thymic selection, a process that is based on TCR amino acid sequences. However, we note that it has been suggested that selection pressure may induce a correlation between TCRs generated and selected [48]. The corresponding statistics of the frequencies of selected TCRs would be modified from those of the generated TCRs shown in Figs. 5. Nonetheless, we assume that the resulting distribution can still be approximated by a single-$\alpha$ model which will not qualitatively alter our conclusions.

Mean-field approximation.

Our mean-field approximation for the mean clone count $c_{k}$ is embodied in Eq. 7, where correlations between fluctuations in the total population $N=\sum_{k}kc_{k}$ in the regulation term $\mu(N)$ and the explicit $c_{k}$ terms are neglected. This approximation has been shown to be accurate for $k\lesssim N^{*}$ when $\bar{\alpha}Q^{2}>\mu(N^{*})$ [34]. The mean-field results overestimate the clone counts for $k\gtrsim N^{*}$. Moreover, when the total steady-state T cell immigration rate is extremely small, the effects of competitive exclusion dominate and a single large clone arises due to competitive exclusion [49, 50, 34]. Nonetheless, an accurate approximation for the steady-state clone abundance $c_{k}$ can be obtained using a variation of the two-species Moran model as shown in [34].

For the naive T cell system, because $Q$ is so large, the mean immigration rate $\bar{\alpha}$ is such that competitive exclusion is not a dominant feature. Moreover, since $N^{*}\gtrsim 10^{10}$, clones counts at comparable sizes are not observed and predicted to be negligible in all models.

The values of $f_{k}^{\rm s}({\rm data})$ become exponentially smaller for large $k$, our inference is most sensitive to the values of $f_{k}^{\rm s}({\rm data})$ for small to modest $k$. The information in the data is primarily manifested by how the $f_{k}^{\rm s}({\rm data})$ decays in $k$, we before the mean-field approximation deviates from the exact solution. Thus, the parameters associated with the human adaptive immune system satisfy the conditions for the mean-field approximation to be accurate, justifying its use in the BDI model.

Steady state assumption.

In this study, we only considered the steady state of our birth-death-immigration model in Eq. 8 because this limit allowed relatively easy derivations of analytical results. This was also the strategy for previous modeling work [35, 4, 6, 7, 34]. However, the per-clone immigration and proliferation times may be on the order of months or years, a time scale over which thymic output diminishes as an individual ages [47, 51, 52, 33]. Indeed, clone abundance distributions have been shown to show specific patterns as a function of age [53, 54, 55].

Although $N(t)$, with fixed $\bar{\alpha}$ and $\bar{r}$ relaxes to steady-state quickly, on a timescale of months, the different subpopulations of specific sizes described by their number $c_{k}$ relax to quasi-steady steady-state across a spectrum of time scales depending on the clone sizes $k$ [56, 34]. The timescales of relaxation of the largest clones can be estimated through the eigenvalues of linearized the Eqs. 7 to be $\sim 10$ years [33].

In addition to time-dependent changes in $\bar{\alpha}Q$, more subtle time-inhomogeneities such as changes in proliferation and death rates have been demonstrated [51, 52]. Thus, our steady-state assumption should be relaxed by incorporation of time-dependent perturbations to the model parameters $\mu^{*}$ and/or $\pi(\alpha,r)$. Longitudinal measurements of clone abundances or experiments involving time-dependent perturbations would provide significant insight into the overall dynamics of clone abundances. Finally, notice that the error slowly decreases as $\lambda,\bar{\alpha}\to 0$, especially for small $w$. The requirement that $\lambda\ll 1$ requires a best fit $\bar{\alpha}$ that is extremely small. At such small $\bar{\alpha}$, the time for the system to reach steady state scales as $\bar{\alpha}$ and the validity of the steady-state approximation is nullified.

Thus, given the multiple known underlying time-inhomogeneous processes that are known to exist, the steady-state may not be strictly reached and our analysis should be considered as an approximation.

Clustered Immigration.

Our mean field model assumed that each immigration event introduced a single naive T cell in the immune system. However, T cells can divide before leaving the thymus and reach a homeostatic state in the periphery. This process can be described by the simultaneous immigration of more than one naive T cell with the same TCR. Clustered immigration of $q$ cells can be implemented in the core model for $c_{k}$ (Eq. 7) via an immigration term of the form $\alpha_{q}(c_{k-q}(\alpha_{q},r)-c_{k}(\alpha_{q},r))$, where $c_{k-q}=0$ for $k-q<0$ (see Appendix 5). For $q>1$, an informative analytic expression for $c_{k}$ is not available. In Fig. S2 in the Appendix, we show the predicted clone abundance $c_{k}$ for a neutral model in which $q=5$. When compared to the case where there is only one cell per immigration, the clone abundance $c_{k}$ will have a larger slope for $k\lesssim q$, making it kink more downward near $k\approx q$. Thus, from Figs. S2 and  9(a), we can see that paired immigration ($q=2$) would increase $f_{k}^{\rm s}$ for $k=2$, providing an improved fitting to data over single copy immigration ($q=1$).

Thus, in addition to appreciable sensitivity of the predicted clone counts to $\pi_{r}(r)$, we also expect clustered immigration defined through the immigration rates $\alpha_{q}$, $q>1$ to control the goodness of fit to data. Indeed, Fig. S2 suggests that the distribution of immigration cluster sizes $q$, in addition to the proliferation rate heterogeneity $w$, is an important determinant of measured clone counts and that $\alpha_{q}$ may be constrained by data. We leave this for future investigation.

General conclusions.

We developed a heterogeneous multispecies birth-death-immigration model and analyzed it in the context of T cell clonal heterogeneity; the clone abundance distribution is derived in the mean-field limit. Unlike previous studies [4], our modeling approach incorporated sampling statistics and provided simple formulae, allowing us to predict clone abundances under different rate distributions for arbitrarily large systems ($N^{*}\sim 10^{10}-10^{11}$), without the need for simulation. The properties of the BDI model and the overall shape of the sampled clone count data renders the first few $k$-values of $c_{k}^{\rm s}$ or $f_{k}^{\rm s}$ the most important for determining the constraints among the model parameters. In other words, only the initial rate of the decrease in $f_{k}^{\rm s}({\rm data})$ for small $k$ governs the quality of fitting to the model, and one should not expect to be able to explicitly infer more than one or two free parameters.

Our heterogeneous BDI model produced mean sampled clone count distributions that we could directly compare with measured clone counts. The unsampled clone counts $c_{k}$ of the neutral model (homogeneous $\alpha$ and $r$) follow a negative binomial distribution which is further modified upon sampling and distribution over the heterogeneous immigration and proliferation rates. Although we determined $\pi_{\alpha}(\alpha)$ through a tool that used recombination statistics inferred from cDNA and gDNA sequences [18, 26], we found that the behavior of the model is rather insensitive to distributions $\pi_{\alpha}(\alpha)$ with mean values $\bar{\alpha}$ much smaller than the largest proliferation rates $r$. The model results are dominated by many low immigration-rate clones and a model that replaces $\alpha$ with its mean value $\bar{\alpha}$ is sufficient.

Conversely, we find that the shape of the clone count profiles $c_{k}$ are quite sensitive to the proliferation rate heterogeneity $w$. A small amount of heterogeneity quickly reduces the best-fit values of $\lambda$ to reasonable values. For estimated values $\eta\sim 10^{-6}-10^{-4}$, $\bar{\alpha}\sim 10^{-4}$, and small values of $\lambda=N^{*}/Q\lesssim 10^{-3}$, requires a best-fit width $w\approx 1$. Heterogeneity is needed to generate clones of sufficiently large size that persist after sampling. Although the number of TCR clones with large proliferation rates $r$ may be small, such clones proliferate more rapidly contributing to higher clone counts at larger sizes. In particular, we found that the shape of expected clone abundance is sensitive to the behavior of the proliferation rate distribution near the maximum dimensional proliferation rate $R$, $\pi_{r}(r\approx R)$.

The predicted clone counts are also modestly sensitive to the distribution of immigration cluster sizes $q$ (representing transient proliferation just before thymic output). When $q>1$ cells of a clone are simultaneously exported by the thymus, the predicted mean clone counts decay much more slowly for small $k\lesssim q$ (see Fig. S2). This modification will allow for better fitting since clustered immigration increases the predicted clone counts for larger $k$, $c_{2}^{\rm s},c_{3}^{\rm s}$, etc., and eventually $f_{2}^{\rm s},f_{3}^{\rm s}$, etc. Thus, we expect that a model containing multiple clustered immigration rates $\alpha_{q\geq 1}$ will lower the error and provide better fitting, particularly at larger $w$. Additional analysis using a distribution of immigration cluster sizes may allow this type of clone count data to reveal more information about the physiological mechanism of naive T cell maintenance.

Even assuming modest heterogeneity, our work lead to the conclusion that proliferation heterogeneity is the more important mechanism driving the emergence of the sampled clone count distributions $c_{k}^{\rm s}$ (and $f_{k}^{\rm s}$) [12]. These results are consistent with the finding that naive T cells in humans are maintained by proliferation rather than thymic output [9]. Since we have only investigated the effects of a uniform distribution for $\pi_{r}(r)$, further studies using more complex shapes of $\pi(\alpha,r)$ can be easily explored numerically using our modeling framework. Different parameter values and rate distributions appropriate for mice, in which naive T cells are maintained by thymic output, should also be explored within our modeling framework.

Acknowledgements

This work was supported by grants from the NIH through grant R01HL146552 and the NSF through grants DMS-1814364 (TC) and DMS-1814090 (MD). The authors also thank the Collaboratory in Institute for Quantitative and Computational Biosciences at UCLA for support to RD.

How heterogeneous thymic output and

homeostatic proliferation shape naive T cell receptor clone abundance distributions

S1: Mathematical Appendices

1 Neutral model

Here, we review the neutral model to provide insight into the properties of our heterogeneous BDI model. When there is no heterogeneity in either proliferation or immigration rates, $\pi(\alpha,r)=\delta(\alpha-\bar{\alpha})\delta(r-\bar{r})$. Upon inserting this expression for $\pi(\alpha,r)$ in Eq. 8, we find that the clone abundance $c_{k}$ follows a negative binomial distribution [10]:

[TABLE]

We can also express $c_{k}/C$, the clone abundance distribution normalized by the mean richness $C$ in the body as

[TABLE]

which is a negative binomial distribution of parameters $\bar{\alpha}/\bar{r}$ and $\bar{r}/\mu(N^{*})$. Using $\bar{\alpha}\approx 1.6\times 10^{-8}/\text{day}$, $\bar{r}\sim 5\times 10^{-4}/\text{day}$, and $\mu(N^{*})\approx 6.4\times 10^{-4}$, we find $\bar{\alpha}/\bar{r}\ll\bar{r}/\mu(N^{*})$. In this regime, $c_{k}/C$, for $k\geq 1$, can be approximated by a log-series distribution with parameter $\bar{r}/\mu(N^{*})$.

To mathematically show that ${c_{k}}/{C}$ converges to a log-series distribution when $\bar{\alpha}/\bar{r}\to 0$, consider a random variable $X$ that follows a negative binomial distribution of parameters $\bar{\alpha}/\bar{r}$ and $\bar{r}/\mu(N^{*})$

[TABLE]

Note that the probability mass function of $X$ is given by $c_{k}/Q$ as can be seen from Eq. S1, the clone abundance distribution for all possible $Q$ clones, which includes $c_{0}$, the number of all clones that are not represented in the organism. To find the clone abundance distribution $c_{k}/C$, for all the $C$ clones present in the organism, we must exclude the case $k=0$ by marginalizing the distribution of $X$ over all $X>0$:

[TABLE]

What remains is to show that the distribution of converges to a log-series distribution of parameter $\bar{r}/\mu(N^{*})$ when $\bar{\alpha}/\bar{r}\to 0$. Consider the moment generating function of $X|X>0$ given by

[TABLE]

Since the moment generating function of a negative binomial distribution $\mathbb{E}\left[e^{\xi X}\right]$ is known, and since $\mathbb{P}\left[X>0\right]=1-\mathbb{P}\left[X=0\right]$ (see Eq. S3), we can write

[TABLE]

For any $x>0$, the limit of $\bar{\alpha}/\bar{r}\to 0$, yields

[TABLE]

If we apply this result to Eq. S6 for $\mathbb{E}\left[e^{\xi X}|X>0\right]$, we find

[TABLE]

which we recognize as the moment generating function of a log series distribution of parameter $\bar{r}/\mu(N^{*})$. Thus, we finally have

[TABLE]

2 Explicit forms using different $\pi_{\alpha}(\alpha)$, $\pi_{r}(r)$

In the following, we propose four simplifying expressions for the heterogeneity-averaged clone counts $c_{k}^{\rm s}(\lambda|\eta)$ derived from Eq. 16.

Clone-independent Neutral model: $\pi(\alpha,r)=\delta(\alpha-\bar{\alpha})\delta(r-\bar{r})$

First, consider the simplest case where all naive T cells carry the same immigration and proliferation rates $\bar{\alpha}$ and $\bar{r}$, respectively, and define $\pi(\alpha,r)=\delta(\alpha-\bar{\alpha})\delta(r-\bar{r})$. This case corresponds to $w\to 0$ and $r\to\bar{r}=1/2$. The self-consistent condition for $\mu^{*}$ and $\bar{\alpha}/\bar{r}$ become

[TABLE]

and the clone count can be explicitly simplified to

[TABLE]

The total sampled clone count is then

[TABLE]

Fixed $\alpha$, distributed $\pi_{r}(r)$:

$\pi(\alpha,r)=\delta(\alpha-\bar{\alpha})\pi_{r}(r)$

Next, consider a common immigration rate $\bar{\alpha}$ for all T cell clones and a box distribution $\pi_{r}$ of full width $w=R=1$. The integral over the dimensionless proliferation rate $r$ is now over $(0,1)$ and $\mu^{*}=(1-e^{-\lambda/\bar{\alpha}})^{-1}$. The averaged clone counts are now explicitly

[TABLE]

The total sampled clone count can also be explicitly expressed the integral over $C^{\rm s}(\bar{\alpha},r,\lambda|\eta)$ from Eq. 12:

[TABLE]

Clone-specific immigration, fixed $r$:

$\pi(\alpha,r)=\pi_{\alpha}(\alpha)\delta(r-\bar{r})$

Finally, using the same rate dimensionalization as before (Eqs. S8), we find explicitly

[TABLE]

where here, the $\bar{\alpha}$ dependence also implicitly arises in $\alpha_{j}$ through Eq. 15. As in Eq. S9, the factor of $2\alpha_{j}$ arises from our choice of nondimensionalization using $2\bar{r}$, twice the naive T cell proliferation rate. Similarly, the total sampled clone count can be explicitly expressed as

[TABLE]

Note that here, $\bar{\alpha}$ is interpreted as the mean immigration rate over the $C_{\star}$ different clones.

3 Small $\bar{\alpha}$ approximation

We show that if the support of $\alpha$ is sufficiently small, terms of order $\alpha/r$ can be neglected in the equation for $c_{k}^{\rm s}$. While $\alpha$ is summed or integrated over, for reasonable distributions $\pi_{\alpha}(\alpha)$, the lowest few rates contribute the most and the average over $\pi_{\alpha}(\alpha)$ can be replaced by a single effective value $\bar{\alpha}$. Even though $r$ is integrated over $(0,1)$, and the region near $0^{+}$ would make $\alpha/r$ large, the contribution from the $c_{k}^{\rm s}(\alpha,r,\lambda,\eta)$ is also small near $r=0$. We have numerically checked that for all cases of $\bar{\alpha}\ll 1/2$, we can approximate $c_{k}^{\rm s}$ by

[TABLE]

This form is derived by approximating the product term as $\sim(k-1)!$ and the exponential term $(\,\cdot\,)^{\alpha/r}\sim 1$. Thus, $f_{k}^{\rm s}$ can also be approximated by

[TABLE]

where $\lambda\equiv N^{*}/Q$ and $\mu^{*}$ is given by

[TABLE]

Note that the $\bar{\alpha}$ prefactor in Eq. S16 derived from $\int_{0}^{\infty}\alpha\pi_{\alpha}(\alpha)\textrm{d}\alpha$ or $\sum_{j}(b_{j}/C_{\star})\alpha_{j}$. Since the rest of the integrand is independent of $\alpha$, the irrelevance of the shape of $\pi_{\alpha}(\alpha)$ is apparent. Only the mean value $\bar{\alpha}$ arises in this approximate calculation of $f_{k}^{\rm s}$.

We have explicitly shown that for small $\bar{\alpha}\ll 1$, this approximation is quantitatively accurate. This simpler form speeds up our numerical analysis and fitting to data. The parameters to infer thus $\theta=\{\bar{\alpha},\lambda=N^{*}/Q,w\}$. Even though $\eta$ is uncertain, it is in principle controlled by the experiment.

Although we evaluate the error between the model and date as a function of the parameters $\theta$, we have not implemented priors for these parameters. We simply investigated the structure of the error $H(\theta)$ and considered physiologically realistic regions in which $\bar{\alpha},\lambda\ll 1$.

4 Correlated immigration $\alpha$ and proliferation $r$.

Hitherto, we have considered independent immigration and proliferation, and assumed a factorisable rate distribution $\pi(\alpha,r)=\pi_{\alpha}(\alpha)\pi_{r}(r)$. However, immigration and proliferation rates may be correlated for certain clones. For example, a frequent realization of V(D)J recombination may also result in a TCR that is more likely to be activated for proliferation. In this case, $\alpha$ would be positively correlated with $r$. In Fig. S1 we consider the effect of correlated $\pi(\alpha,r)$. For $\bar{r}/2\leq r\leq 2\bar{r}$, we considered normalized, positively/negatively correlated box distributions as shown in Fig S1(a):

[TABLE]

Within our mean field model, these correlated distributions $\pi(\alpha,r)$ result in very similar expected clone abundance distributions $c_{k}$ (Fig S1(b)). This insensitivity to correlations between immigration and proliferation can be qualitatively understood by considering the “line integral” over dominant paths of $\pi(\alpha,r)$ in the $uv/(1-v)=\alpha/(\mu(N^{*})-r)$ vs. $u=\alpha/r$ diagram. As shown in Fig. S1(c), both line integrals remain in the log-series distribution regime, indicating that the clone abundance distributions are qualitatively similar to that predicted by a model with proliferation heterogeneity alone.

5 $c_{k}$ for clustered immigration

We explore how clustered emigration from the thymus affects the mean clone count $c_{k}$. Suppose that $q$ cells of the same clone (TCR nucleotide or amino acid sequence) are simultaneously exported by the thymus. The equation for the mean clone count $c_{k}$ is thus

[TABLE]

This equation does not admit a simple analytic solution so we numerically solved the equation assuming $\alpha_{q}=\alpha_{5}\mathds{1}(q,5)$ and $Q=10^{11}$. The shape of $c_{k}$ for single cell immigration ($q=1$) is compared to that of $q=5$.

For $q>1$, $c_{k}$, and ultimately $c_{k}^{\rm s}$ and $f_{k}^{\rm s}$ is flatter up to $k\approx q$, making the clone counts kink more downwards near $q$. Thus, as can be seen from Fig. 9(a,b), we can reasonably conclude that some level of paired immigration would provide even better fits to the data at appropriately small values of $\lambda$, especially for the first few $k$-points.