How a well-adapting immune system remembers

A Lymphocyte Dynamics for Approximating Optimal Sequential Inference.

For concreteness, we consider a drift-diffusion model of environmental change (Eq. 11). The drift-diffusion model, while clearly a much-simplified model of real evolution, captures two key features of changing pathogenic environments: the coexistence of diverse pathogens and the temporal turnover of dominant pathogen strains. The aim of this model of pathogen evolution is not to provide a realistic description of short-term pathogen dynamics within a population such as during an epidemic, but rather to capture overall features of the long-term dynamics of many pathogens over a host’s lifetime. The drift-diffusion model is mathematically equivalent to a classical neutral stochastic evolution of pathogens (25) driven by genetic drift happening on a characteristic time scale τ and immigration from an external pool with immigration parameters θ=(θ1,…,θK) (Eq. 11). Generally, pathogens are under selective pressure to evade host immunity, and strains are replaced faster than under the sole action of genetic drift. Matching the time scale of pathogen change to those observed experimentally, the model then provides a simple, effective description of the pathogen dynamics.

Can a plausible dynamics of lymphocyte receptor clones approximate the optimal repertoire dynamics? In particular, is there an approximate autonomous dynamics for the repertoire composition, which does not require access to the full latent high-dimensional belief distribution B(Q,t)? Consider the simple case of a logarithmic cost function, c(P̃a)=−lnP̃a, and uniquely specific receptors, fa,r=δa,r, in which the optimal match between the receptor and pathogen distributions (Eq. 1) is the identity, Pa*=Q^a (ref. 15 and Materials and Methods).

For the pathogen dynamics of Eq. 10, we show, using a decomposition of B(Q,t) into a mixture of Dirichlet distributions (SI Appendix, section 1B), that there exists an approximate but autonomous dynamics for the number of lymphocyte receptors N. It takes the following simple form (parallel to Eqs. 2–4):Na(0)=χθa for all a,[5]Na(t+)=Na(t−)+χ when a is encountered,[6]τdNadt=−12(χ−1|N|−1)(Na−χθa) for all a,[7]where χ is a scale factor controlling the total population size and where P=N/|N|. Here and in the following, we use |x|=∑ixi to denote the l1-norm of a vector x. We compared these approximate dynamics to an exact solution computed by spectrally expanding the generator of the stochastic dynamics (SI Appendix, section 1C) and found that they approximate the optimal dynamics closely (Fig. 2).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

A plausible repertoire dynamics can implement approximate Bayesian inference. Frequency Qa of a pathogen over time (blue line) in a dynamically changing environment along with the exact optimal repertoire dynamics (orange line) and the approximate dynamics based on Eqs. 5–7 (Pa, green line). The approximate dynamics closely follows the exact Bayesian dynamics implemented as described in SI Appendix, section 1C (Pa*, orange line). Parameters are as follows: K=500, θ=0.02, τ=20 year, λ=10/year, c(P̃a)=−lnP̃a, and fa,r=δa,r.

These dynamics have a plausible biological implementation. Each pathogen encounter leads to a fixed increment χ of the number of specific lymphocytes (Eq. 6), implying a regulation mechanism that controls the number of cell divisions upon clonal expansion as a function of precursor frequency. Through thymic output, the repertoire starts at and is renewed by the naive repertoire—encoded in χθa (Eqs. 5 and 7). By definition, this naive repertoire is optimally tuned to the baseline frequencies of pathogens θa/|θ|. The total number of lymphocytes |N| regulates the time scale of cell decay and renewal (Eq. 7), which can be done through some global homeostatic mechanisms.

To gain intuition about what the dynamics of Eqs. 5–7 mean, let us define the rescaled variables na=χ−1Na. Eq. 6 increases na by 1 whenever pathogen a is sampled. n can thus be interpreted as counting the number of times the different pathogens were sampled and their initial value na(0)=θa as pseudocounts corresponding to the prior belief in the pathogen frequencies. In a changing environment, experience gained from previous encounters eventually loses its predictive power, and thus Eq. 7 discounts these old encounters, giving more relative weight to recent ones. If there are no new pathogen encounters, na eventually relaxes back to the prior θa.

The assumption that Pa*=Q^a can be relaxed, provided that the nonlinear relation G between P and Q (Eq. 1) remains approximately invertible, leading to autonomous dynamics of the lymphocyte clone sizes Na(t) similar to Eqs. 5–7 (SI Appendix, section 2C).

Learnability of Pathogen Distribution Implies a Sparse Pathogenic Landscape.

The immune system must be prepared to protect us not just from one pathogen but a whole distribution of them. Even restricting recognition to short peptides and accounting for cross-reactivity (26), estimates based on precursor frequencies for common viruses give an effective antigen environment of size K∼105–107 (27). How can the immune system learn anything useful about such a high-dimensional distribution from a limited number of pathogenic encounters? Naively, one might expect that the number of samples needed to learn the distribution of pathogens must be larger than of the number of pathogens, i.e., λt∼K, where t is the time over which learning takes place. Although little is known about the receptor–antigen encounter rate λ, this estimate suggests that the pathogenic environment is not easily learnable, and therefore memory has limited utility.

This apparent paradox can be resolved by the fact that the pathogenic environment may be sparse, meaning that only a small fraction of the possible pathogens are present at any given time. In our model of the pathogen dynamics, this sparsity is controlled by the parameter θ. In the scenario that we are considering, typical pathogen landscapes Q are drawn from the steady-state distribution ρs(Q) of the immigration-drift dynamics, which is a Dirichlet distribution parametrized by θ (SI Appendix, section 1, Eq. S1). When θa is small, the distribution is peaked at Qa=0, meaning that pathogen a is absent a majority of the time. For instance, for uniform θa≡θ≪1, the effective number of pathogens present at any given time is Kθ (SI Appendix, section 3C). Since the system only needs to learn about the pathogens that are present, the condition for efficient learning should naively be λt∼Kθ, which is much easier to achieve realistically for small θ.

Our theory can be used to quantify the benefit of memory as a function of the different immunological parameters. We compute the optimized cost function c(t)=∑aQa(t)c(P̃a⋆(t)) and study how it decreases as a function of age, t, relative to the cost at birth c0=c(t=0), as the organism learns from pathogen encounters. This relative cost also depends on the encounter rate λ, the size of the pathogenic space K, the sparsity of pathogenic space θ, and time scale of change in the environment τ. Fig. 3A shows that the benefit of memory increases with pathogen sparsity—when θ is small, even a few encounters suffice to seed enough memory to reduce the cost of future infections. The cost saturates with age to a value c∞, either because memory approaches optimality or because memory eventually gets discarded and renewed as the environment changes. Fast-changing environments lead to an earlier and higher saturation of the cost with age (Fig. 3B) since learning and prediction are limited by decorrelation of the environment. The pathogen dynamics is sped up when there are strong selection pressures to evade immunity. Faster dynamics decrease learning efficiency and in turn reduce selective pressures. The effective time scale should in practice be set by a coevolutionary balance between both effects.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Advantage of immunological memory depends on sufficient sampling. The mean expected cost of an infection in a changing environment is a function of the age of the organism t, the time scale τ on which the environment changes, and the sparsity 1/θ of the environment. (A) Relative cost as a function of age for environments with different sparsity (for fixed λτ=105, K = 2,000). (B) Relative cost as a function of age for environments changing more or less rapidly (for fixed θ=0.02, K = 2,000). (C) Collapse of the data by plotting the relative cost gap against the number of samples per effective (eff.) dimension. Simulations with logarithmically spaced parameters in the ranges λt from 1 to 5,000, λτ from 1 to 1,000,000, K from 100 to 2,000, and θ from 0.01 to 0.2 (color-coded as in A). Cost is c(P̃a)=1/P̃a, and receptors are assumed to be uniquely specific, fa,r=δa,r.

Analytical arguments show that in the limit of few samples the relative cost c/c0 achievable in a static environment scales as λt/Kθ (SI Appendix, section 3). In general, we find that the cost is a function of λte/Kθ, where the effective time te is defined via λte=|n(t)|−|n(0)|≈2λτtanh(λt/2λτ) with n(t) being the vector of the encounter counts discussed above (SI Appendix, section 1D for derivation from Eqs. 6 and 7). Plotted in terms of this variable, the relative cost gap as function of time, (c−c∞)/(c0−c∞), collapses onto a single curve for all parameter choices (Fig. 3C). Fig. 3C shows that the cost drops by a factor of ∼2, when λte/Kθ∼1. Thus, there is a substantial benefit to memory already when the effective number of encounters is comparable to the effective number of pathogens. At young ages (small t) or with slowly changing environments (large τ), te≈t, and so this condition is simply λt∼Kθ, i.e., the total number of encounters should be comparable to effective number of pathogens that are present.

Optimal Attrition Time Scale.

Our theory suggests that there is an optimal time scale for forgetting about old infections which is related to the time scale over which the environment varies. Eq. 7 shows that memory should optimally be discounted on an effective time scale τmem=2τ/(|n|−1). Comparing this to the slowest time scale of environmental variation, τc=2τ/|θ| (SI Appendix, Eq. 37), where |θ|=∑aθa, we haveτmem=τc|θ||n|−1.[8]The time scale on which old memories should be forgotten scales with the environmental correlation time scale. The two time scales are equivalent when the immune system has little information about the pathogenic environment (|n|∼|θ|). Given the long time scales over which many relevant pathogens change, immune memory should generally be long-lived (with the time scale of decay being of the order of years or decades). Indeed, despite the relatively short life span of memory cells (28), constant balanced turnover keeps elevated levels of protection for decades after an infection, even in the absence of persistent antigens (29⇓–31).

Interestingly, our theory predicts that memory should be discounted more quickly when the immune system has gathered more information (larger |n|). Using the mean-field equations for |n(t)|−|n(0)| from the previous section, we can derive how the memory time scales at steady state at high sampling rate. Using that for large times |n(t)|≫|n(0)| holds in the high sampling rate limit, one can simplify the mean-field result to |n|∼2λτ in steady state (t→∞). Combined with Eq. 8, τmem=τc|θ|/λ follows, which shows that a larger sampling rate leads to a faster discounting of past evidence. This is reminiscent of results in optimal cellular signaling where there are similar trade-offs between noise averaging and responsiveness to changes in the input signal (32).

Memory Production in Sparse Environments Should Be Large and Decrease with Prior Exposure.

The theory can be used to make quantitative and testable predictions about the change in the level of protection that should follow a pathogen encounter. Consider an infection cost function that depends as a power law on the coverage, c(P̃a)=1/P̃aα, with a cost exponent α that sets how much attention the immune system should pay to recognizing rare threats. (Below, we will use the shorthand α=0 to indicate logarithmic cost.) Cost functions of this form can be motivated by considering the time to recognition of an exponentially growing antigen population by the immune system (15) or, alternatively, by considering the time delay of the expansion of the precursor cells to some fixed number of effector cells (SI Appendix, section 4).

In the simplest model for repertoire updates, recognition of pathogens leads to proliferation proportionally to the number of specific precursor cells, followed by a homeostatic decrease of the memory pool (18, 33). Thus, the fold change P̃a(t+)/P̃a(t−)=const where t−, t+ are times just before and after the encounter. By contrast, our Bayesian theory predicts that the fold change upon encountering pathogen a should beP̃a⋆(t+)P̃a⋆(t−)=(1+κ/P̃a⋆(t−)(1+α))1/(1+α),[9]where κ depends on prior expectations about the antigenic environment and previous pathogen encounters (Materials and Methods). Setting α=0 gives the result for a logarithmic cost function.

To understand this prediction, first consider the effect of a primary infection on a naive repertoire, θa≡θ, P̃⋆(0)=1/K, and |n(0)|=Kθ, where the receptors are uniquely specific (fa,r=δa,r). In this case, κ=1/K1+αθ (SI Appendix, section 2A) and Eq. 9 predicts a fold change of (1+1/θ)1/(1+α). We have argued previously that their learnability implies that pathogenic environments are sparse, i.e., θ≪1. Therefore, we predict that primary antigenic encounters should lead to a large memory production. Experimentally, memory production typically leads to the proliferation of antigen-specific cells by a factor of 100- to 1,000-fold (14), in qualitative agreement with this prediction. Turning the argument around, such a large increase in protection upon an encounter is only adaptive in highly sparse environments. Quantitatively, it implies a sparsity parameter θ∼10−6–10−4 (here, taking α=1 for definiteness) (Fig. 4B). Combined with the estimate K∼105–107 (27), this suggests that the effective number of pathogens at any given time ranges from Kθ=0.1 to 1,000.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Changes in protection levels upon infections for cost functions c(P̃a)=1/P̃aα (α=0 indicates logarithmic cost). (A) The prediction of the Bayesian model (Eq. 9 for different α) closely fits experimental data on antibody titers, while a constant fold-change model (which would correspond to a straight horizontal line) does not. Data on prevaccination and postvaccination antibody titers against stem and head hemagglutinin epitopes following a booster vaccination with inactivated H5N1 in humans from ref. 34. (B) Optimal fold change of coverage for ni(t−)=θ for all i. The fold change increases with the sparsity of the environment controlled by parameter θ. (C) Absolute (Abs.) change in coverage for primary infection and reinfections (for α=1), normalized such that the change for the primary infection is 1 and neglecting attrition.

To test Eq. 9 on immunological data, we fit the Bayesian update model to experiments reporting fold changes in antigen titers upon booster vaccinations against influenza from ref. 34 (Fig. 4A) using least squares. Titers correspond to the concentration of antibodies that are specific to the antigen a and can thus be viewed as an experimental estimate of P̃a. The optimal Bayesian strategy explains the data, accounting for the larger boosting at small prevaccination titers and showing no increase for large titers, while the proportional model predicts constant boosting for all titers. Similar experimental results have been reported for antibody titers before and after a shingles vaccination (35). Mechanistic models have been proposed to explain how the population dynamics of expanding lymphocytes might give rise to nonproportional boosting for both B and T cells (22, 33, 36, 37).

Interestingly, for T cells, Quiel et al. (38) have shown that fold expansion to peak cell numbers in an adoptive transfer experiment depends on the initial number of T cells as a power law with exponent ∼−1/2. That scaling, which is for the peak expansion, predicts more expansion at high precursor number than Eq. 9, which is for memory production. This implies a nonlinear relationship between peak T-cell level and memory production, which further suppresses memory production at high precursor numbers. This prediction might be checked in experiments measuring memory production after infection clearance, as well as the expansion peak.

Long-Term Dynamics of a Well-Adapting Repertoire.

Our model makes predictions for the dynamics of growth and attrition of memory over time, with consequences for immunity and for the diversity of the immune repertoire. We quantify the dynamics in terms of a memory fraction defined as a sum of the coverage fractions P̃ai over all previously encountered pathogens {ai}. The memory fraction measures the size of memory relative to the size of the whole immune repertoire. Early in life, every infection is new, and even modest increases in the memory fraction lead to large drops in infection susceptibility (measured by the expected cost of new infections in Fig. 5A). At the same time, the memory fraction increases rapidly (Fig. 5B), but the growth of memory slows as subsequent infections lead to less memory production following the optimal fold-change rule in Eq. 9 and as attrition begins to play a role. The fraction of the repertoire devoted to memory in midlife is largely determined by how the cost of infections scales with coverage. The observed memory fraction of ∼50% at midlife suggests a cost exponent of α≈0.5 (Fig. 5B). The diversity of the memory repertoire increases with time at a rate that slows with age (quantified in Fig. 5C by richness, which measures the number of unique specificities and the Shannon entropy of the repertoire frequency distribution).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Relative cost (A and D), memory-cell fraction (B and E), and memory diversity (C) as a function of age. A–C show long-term dynamics for a repertoire following optimal Bayesian update dynamics for three different α. Memory diversity is plotted as richness, i.e., the number of unique memory specificities, as well as the exponential of Shannon entropy S of the memory compartment defined for a probability distribution pi with ∑ipi=1 as S=−∑pi⁡lnpi. D and E compare the optimal Bayesian dynamics with a constant fold change update by a factor of 30 and the same proportional update but with a cap to 105 of the total fold expansion for any clone. (F) Comparison of the cost as a function of age for the optimal dynamics (solid lines) and a dynamics without attrition (dotted lines) for two environmental correlation time scales τc=2τ/Kθ (SI Appendix, section 1C). Parameters are as follows: encounter rate λ=40/year, antigen space dimensionality K=105, and antigen sparsity θ=2.5⋅10−4. In A–E, we used τ→∞, and in D–F, α=0.5. To reduce fluctuations, the statistics were averaged over multiple runs of the dynamics. Data in B are from refs. 39 and 40.

To gain insight into these dynamics of our model, we average the stochastic equations over the statistics of pathogen encounters. We show in SI Appendix, section 2B that this mean-field approximation yields a differential equation for the population fraction of different clones with two opposing contributions which balance alignment of the immune repertoire with the current pathogenic environment (i.e., memory production) against alignment with the long-term mean environment (i.e., attrition). Interestingly, the mean-field equation broadly coincides with dynamics that were proposed in ref. 15 to self-organize an optimal immune repertoire. The essential difference here is that the time scale of learning slows down with increasing experience following the rules of optimal sequential update in Eq. 9.

We then asked which features of the proposed repertoire dynamics are most relevant to ensure its effectiveness. How important is the negative correlation between fold expansion and prior immune levels, and how important is attrition? Furthermore, if the immune system follows Bayesian dynamics, it must have integrated on an evolutionary time scale a prior about composition and evolution of the pathogen environment through the parameters θ and τ—however, the prior may be inaccurate. How robust is the benefit of memory to imperfections of the host’s prior assumptions about pathogen evolution? To answer these questions, we compare the long-term immune-repertoire dynamics using the optimal Bayesian scheme to other simplified schemes. We find that a constant fold expansion dynamics quickly leads to very suboptimal repertoire compositions (Fig. 5D, pink line), since the exponential amplification of cells specific to recurrent threats quickly leads to a very large fraction of the repertoire consisting of memory of those pathogens (Fig. 5E, pink line). This suboptimality persists, even if we assume that some global regulation caps the constant fold expansion such that no individual receptor clone can take over all of the repertoire (Fig. 5 D and E, gray line). Thus, negative feedback in T-cell expansion to individual antigens is very important to maintain a properly balanced diverse repertoire. In contrast, within dynamics with a negative correlation, the precise levels of updating do not need to be finely tuned to the environmental statistics: Varying the assumed sparsity of the pathogen distribution, which controls fold expansion upon primary infection in the optimal dynamics, leads to a relatively modest deterioration of the convergence speed of the learning process (SI Appendix, Fig. S2A) and does not matter asymptotically (SI Appendix, Fig. S2B). Attrition does not matter at young age, but can play an important role for long-term adaptation to relatively rapidly changing pathogen distributions (Fig. 5F). However, the attrition time scale need not be finely tuned to get close to optimal dynamics (SI Appendix, Fig. S3).

Adapting a Cross-Reactive Repertoire.

Above, we described adaptation of immune repertoires in terms of changes in the effective coverage P̃a=∑rfa,rPr, where the cross-reactivity matrix F=(fa,r) reflects the ability of each receptor to recognize many antigens and also the propensity of each antigen to bind to many receptors (26). Because of cross-reactivity, each pathogen encounter should result in the expansion of not just one but potentially many receptor clones. Here, we ask how the optimal immune response is distributed among clones with different affinities.

Following Perelson and Oster (41), we will represent the interaction of receptors and antigens by embedding both in a multidimensional metric recognition “shape space,” where receptors are points surrounded by recognition balls. Antigens that fall within a ball’s radius will be recognized by the corresponding receptor. In this presentation, a and r are the coordinates of antigens and receptors, respectively, and their recognition propensity depends on their distance, fa,r=f(|a−r|).

Earlier sections have already discussed the optimal dynamics of the coverage P̃a, which is a convolution of the cross-reactivity matrix with the receptor clone distribution Pr. Thus, the optimal dynamics of the clone distribution can be derived by deconvolving the cross-reactivity subject to the constraint that Pr cannot be negative. Carrying out this analysis in SI Appendix, section 2C reveals a general qualitative phenomenon—competitive exclusion between clones expressed in the repertoire and their close neighbors within the cross-reactivity radius (SI Appendix, Fig. S4, blue line). This exclusion is not an assumption of the model, but rather stems from the optimal Bayesian theory. Given a receptor clone that covers one region of antigenic shape space, the global likelihood of detecting infections increases by placing other clones to cover other regions. This can be shown analytically when cross-reactivity is limited, memory updates are small in magnitude, and the pathogen distribution is assumed to be uncorrelated (SI Appendix, section 2C).

In general, the frequencies of pathogens might be correlated in antigenic space, for example, because mutations from a dominant strain give rise to new neighboring strains. An optimally adapting immune system should incorporate such correlations as a prior probability favoring smoothness of the pathogen distribution. Such priors work their way through the optimal belief update scheme that we have described and weaken the competitive exclusion between clones with overlapping cross-reactivity (SI Appendix, Fig. S4, orange line).

In general, when cross-reactivity is wide or the required clone fraction update is large, numerical analysis shows that achieving optimally predictive immunity after a pathogen encounter requires a global reorganization of the entire repertoire (SI Appendix, Fig. S5, blue line). There is no plausible mechanism for such a large-scale reorganization since it would involve up- and down-regulation, even of unspecific clones. However, in SI Appendix, section 2C we show that the optimal update can be well-approximated by changes just to the populations of specific clones with pathogen-binding propensities fa,r that exceed a threshold. The optimizing dynamics with this constraint exhibit strong competitive exclusion, where only the highest-affinity clones proliferate, while nearby clones with lower affinity are depleted from the repertoire (SI Appendix, Fig. S5, orange line). The local update rule provides protection that comes within 1% of the cost achievable by the best global update. Thus, reorganization of pathogen-specific receptor clone populations following an infection, as seen in vertebrates, can suffice to achieve near-optimal predictive adaptation of the immune repertoire.