Quantitative Immunology for Physicists

TL;DR

This paper reviews recent quantitative approaches in immunology, highlighting how physicists model the complex, multiscale adaptive immune system to predict immune responses and understand its dynamics.

Contribution

It provides an overview of recent theoretical and experimental contributions by physicists to modeling the adaptive immune system, emphasizing quantitative methods.

Findings

Physicists have developed models predicting immune cell behavior.

Quantitative methods help understand immune system dynamics.

The review summarizes theoretical tools used in the field.

Abstract

The adaptive immune system is a dynamical, self-organized multiscale system that protects vertebrates from both pathogens and internal irregularities, such as tumours. For these reason it fascinates physicists, yet the multitude of different cells, molecules and sub-systems is often also petrifying. Despite this complexity, as experiments on different scales of the adaptive immune system become more quantitative, many physicists have made both theoretical and experimental contributions that help predict the behaviour of ensembles of cells and molecules that participate in an immune response. Here we review some recent contributions with an emphasis on quantitative questions and methodologies. We also provide a more general methods section that presents some of the wide array of theoretical tools used in the field.

Figures

Tables

Table 1. Table 1: Summary of different neutral models discussed in the review, with their clone size distributions. ν < μ 𝜈 𝜇 \nu<\mu refer to the division and death rates respectively, and ρ = 1 / ( 1 − α ) 𝜌 1 1 𝛼 \rho=1/(1-\alpha) , where α 𝛼 \alpha the probability of mutation upon division in the Yule process.

Neutral model	$N_C$	Comment	Eq.
Lymphocyte dynamics	$\sim {\left(\mu /\nu \right)}^C/C$	Source	Eq. 63
Population genetics	$\sim 1/C$	Fixed pop. size	Eq. 141
Yule process	$\sim 1/C^{1+\rho }$	Expanding pop.	Eq. 207

Equations

$${\rm Ligand}+{\rm Receptor}\begin{array}[]{c}k_{\rm on}\\ \rightleftharpoons\\ k_{\rm off}\end{array}{\rm Complex}\longrightarrow{\rm Activation}$$

$$k_{\mathrm{diffusion}}=4\pi \left(R_{\mathrm{Ligand}}+R_{\mathrm{Receptor}}\right)\left(D_{\mathrm{Ligand}}+D_{\mathrm{Receptor}}\right),$$

$$D_{\mathrm{Ligand}}=D_{\mathrm{solution}}=\frac{k_{B}T}{6\pi \eta R_{\mathrm{Ligand}}}\approx 300\mu m^2{s}^{-1}$$

$$k_{\mathrm{diffusion}}=4\pi R_{\mathrm{Ligand}}{D}_{\mathrm{Ligand}}=\frac{2k_{B}T}{3\eta }\approx 2\cdot 10^9{M}^{-1}{s}^{-1}.$$

$$k_{\mathrm{diffusion}}=4\pi \left(R_{\mathrm{Ligand}}+R_{\mathrm{Receptor}}\right)\left(D_{\mathrm{solution}}\right).$$

$$k_{\mathrm{on}}=k_{\mathrm{diffusion}}{e}^{-\Delta G_{\mathrm{association}}/\left(k_{B}T\right)},$$

$$[R{]}_{\mathrm{total}}=\frac{\#R_{\mathrm{cell}}{N}_{\mathrm{cell}}}{V{\mathcal{N}}_{\mathrm{Avogadro}}},$$

$$[\mathrm{IL}-7R{]}_{\mathrm{total}}\approx 3\cdot 10^{-10}M=300pM,$$

$$\Delta G(a,r)=\sum _{k=1}^N{J}_k(a_k,r_k),$$

$$S(a)=\sum _{a^{\prime }\in D}{e}^{k{s}_{a,a^{\prime }}},$$

$$I(a)=\frac{S(a)}{S_0+S(a)}.$$

$$C_i=C_0{\left(\frac{k_{\mathrm{P}}}{k_{\mathrm{P}}+k_{\mathrm{off}}}\right)}^i.$$

$$C_{\mathrm{total}}\sim \frac{(k_{\mathrm{on}}/(V{\mathcal{N}}_A))TL}{k_{\mathrm{off}}+(k_{\mathrm{on}}/(V{\mathcal{N}}_A))T}.$$

$$\frac{C_{N_{\mathrm{p}}}(\mathrm{agonist})}{C_{N_{\mathrm{p}}}(\mathrm{non}-\mathrm{agonist})}={\left(\frac{{\tau }_{\mathrm{off}}(\mathrm{agonist})}{{\tau }_{\mathrm{off}}(\mathrm{non}-\mathrm{agonist})}\right)}^{N_p}.$$

$$\begin{array}{rl}\mathcal{F}& =\sum _i\frac{{\lambda }_i}{2}\int dxdy{C}_i(x,y){\left[z(x,y)-z_i^{(0)}\right]}^2\\ & +\frac{1}{2}\int dxdy\left[\gamma (\nabla z(x,y){)}^2+\kappa ({\nabla }^{2}z(x,y){)}^2\right],\end{array}$$

$$\frac{\partial C_i}{\partial t}$$

$$\frac{\partial z}{\partial t}=-M\frac{\delta \mathcal{F}}{\delta z}+\zeta ,$$

$$\frac{d}{dt}\#Recepto{\mathrm{r}}^{\ast }=k_{\mathrm{on}}[\mathrm{Cytokine}]\#Receptor-k_{\mathrm{off}}\#Recepto{\mathrm{r}}^{\ast },$$

$$\mathrm{JA}{\mathrm{K}}^{\ast }={\mathrm{JAK}}_{\mathrm{total}}\frac{[\mathrm{Cytokine}]}{K_D+[\mathrm{Cytokine}]},$$

$$\frac{d}{dt}[\mathrm{pSTAT}]=k_{\mathrm{cat}}[\mathrm{JA}{\mathrm{K}}^{\ast }]\frac{[\mathrm{STAT}]}{[\mathrm{STAT}]+K_m}-k_{\mathrm{dephos}}[\mathrm{pSTAT}],$$

$${\rm L+R_{1}+R_{2}\begin{array}[]{c}k_{1}^{\rm on}\\ \rightleftharpoons\\ k_{1}^{\rm off}\end{array}LR_{1}+R_{2}\begin{array}[]{c}k_{2}^{\rm on}\\ \rightleftharpoons\\ k_{2}^{\rm off}\end{array}LR_{1}R_{2}.}$$

$$K_1[L{R}_1]$$

$$K_2[L{R}_1{R}_2]$$

$$[L{R}_1{R}_2]=\frac{[R_2{]}_{\mathrm{total}}}{1+\frac{K_2}{[R_1{]}_{\mathrm{total}}}\left(1+\frac{K_1}{[L]}\right)}.$$

$$[L{R}_1{R}_2]=[R_2{]}_{\mathrm{total}}\frac{[L]}{[L]+\frac{K_1{K}_2}{[R_1{]}_{\mathrm{total}}}}.$$

$$\frac{d}{dt}[\mathrm{Cytokine}]={\kappa }_{\mathrm{prod}}-{\kappa }_{\mathrm{consum}}.$$

$$\frac{\partial c(\vec{r},t)}{\partial t}=D{\nabla }^{2}c(\vec{r},t)+{\kappa }_{\mathrm{prod}}(\vec{r},t)-{\kappa }_{\mathrm{consum}}(\vec{r},t),$$

$${\kappa }_{\mathrm{consum}}=k_{\mathrm{on}}{N}_{\mathrm{R}}{n}_{\mathrm{consum}}c(\vec{r},t)=k_{c}c(\vec{r},t),$$

$$\frac{\partial c(\vec{r},t)}{\partial t}=D\left[{\nabla }^{2}c(\vec{r},t)-{\xi }^{-2}c(\vec{r},t)\right],$$

$$\tau (r)=\frac{1}{2k_c}\left(1+\frac{r}{\xi }\right).$$

$$c(r)=c(R_{\mathrm{cell}})\frac{R_{\mathrm{cell}}}{r}{e}^{(R_{\mathrm{cell}}-r)/\xi },$$

Full text

1. I Introduction

2. II Physical chemistry of ligand-receptor interaction: specificity, sensitivity, kinetics.

3. II.1 Diffusion-limited reaction rates

4. II.2 Extrapolating collision rates in solution to association rates in the physiological context

5. II.3 Rates and numbers in the physiological context

6. II.3.1 Association rates

7. II.3.2 Dissociation rates

8. II.3.3 Numbers of receptors per cell

9. II.3.4 Typical binding rates and some biology

10. II.4 Receptor-antigen specificity

11. II.4.1 Cross-reactivity

12. II.4.2 Models of receptor-antigen binding

13. II.4.3 Data-driven receptor-antigen binding models

14. II.4.4 Modeling immunogeneticity

15. III Antigen discrimination

16. III.1 T cells

17. III.1.1 Kinetic proofreading for ligand discrimination

18. III.1.2 Adaptive kinetic proofreading

19. III.1.3 Coupling mechanics and biochemistry: the significance of forces for ligand discrimination

20. III.2 B cells

21. III.3 Coarse-graining of molecular details and model reduction

22. IV Cell-to-cell communication through cytokines

23. IV.1 Cytokine signaling and the JAK-STAT pathway

24. IV.1.1 Cytokine binding and signaling at equilibrium

25. IV.1.2 Tunability of cytokine responses.

26. IV.1.3 Regulation by cytokine consumption

27. IV.1.4 Other regulations

28. IV.2 Communication across space and time, and cytokine niches

29. IV.2.1 Cytokine dynamics as diffusion-degradation

30. IV.2.2 Screening by cytokine-consuming cells

31. IV.2.3 Probability of autocrine capture

32. IV.2.4 Size of cytokine niches

33. V Cell fate

34. V.1 Gene regulation and cell differentiation

35. V.1.1 Th1/Th2 differentiation

36. V.1.2 Other differentiation switches

37. V.1.3 Experimental test of bistability in cell differentiation

38. V.2 Hematopoiesis

39. V.2.1 Timescales

40. V.3 Inferring the hematopoiesis differentiation tree

41. V.4 Cell fate during the immune response

42. V.4.1 Choice and timing of cell fate under stimulation

43. V.4.2 An aside on cell types

44. V.4.3 Inferring cell fate timelines during the immune response

45. V.4.4 Quorum sensing

46. VI Repertoires

47. VI.1 Size of immune repertoires

48. VI.2 Inference of the stochastic repertoire generation process

49. VI.3 Thymic selection and central tolerance

50. VI.4 Diversity and the clone size distribution

51. VI.5 Repertoire sharing

52. VI.6 Optimal receptor distribution

53. VI.7 Repertoire response to an immune challenge

54. VII Lymphocyte population dynamics

55. VII.1 Neutral dynamics

56. VII.2 A note about “neutral processes”

57. VII.3 Population dynamics model with external signals

58. VII.4 In host HIV dynamics.

59. VIII Affinity maturation

60. VIII.1 Modeling the hypermutation process

61. VIII.2 Cycles of selection

62. VIII.3 Evolution of broadly neutralizing antibodies

63. VIII.4 Population genetics approaches to affinity maturation

64. VIII.4.1 Evolutionary analysis of repertoire dynamics

65. VIII.4.2 Models of co-evolution of phenotypic traits

66. IX Population dynamics of pathogens and hosts

67. IX.1 Viral fitness models

68. IX.2 Co-evolution between host and pathogen populations

69. X Discussion

70. XI Glossary

71. XII Methods

72. XII.1 Physical kinetics

73. XII.1.1 Diffusion-limited reaction rate

74. XII.1.2 The rates of dissociation between two biomolecules

75. XII.1.3 The formation of ligand-receptor pairs: equilibrium and kinetics

76. XII.2 Gene regulation

77. XII.2.1 Basic model

78. XII.2.2 Auto-amplification with a single transcription factor

79. XII.2.3 Auto-amplification with multiple transcription factors

80. XII.3 Population dynamics, genetics

81. XII.3.1 Deterministic mutation-selection balance

82. XII.3.2 Genetic drift

83. XII.3.3 Wright-Fisher model

84. XII.3.4 Probability of extinction

85. XII.3.5 Moran model, continuous limit, and time varying selection

86. XII.3.6 Branching processes

87. XII.3.7 Coalescence process

88. XII.3.8 Site frequency spectra and tree balancing

89. XII.3.9 Clonal interference

90. XII.3.10 Quantitative traits

91. XII.3.11 Lineage reconstruction

92. XII.3.12 Population growth rates

93. XII.4 Ecological models

94. XII.4.1 Generalized Lotka–Volterra models

95. XII.4.2 Susceptible-Infected-Recovered (SIR) Models

96. XII.4.3 Solution of stochastic population dynamics with a source

97. XII.4.4 Solution to foward jump process with opposing drift and source

98. XII.4.5 The Yule process

99. XII.5 Inference

100. XII.5.1 Probabilistic inference, maximum likelihood and Bayesian statistics

101. XII.5.2 Model selection

102. XII.5.3 Expectation-Maximization

103. XII.5.4 Hidden Markov models

104. XII.5.5 Information theory

105. XII.5.6 Maximum entropy models

106. XII.5.7 Machine learning and Neural networks

107. XII.6 High-throughput repertoire sequencing

108. XIII Acknowledgments

I Introduction

The role of the immune system is to detect potential pathogens, confirm they really are undesirable pathogens, and destroy them. The goal is in principle well defined. However recognizing molecular friends from foes is not easy, and organisms have evolved many complementary ways of dealing with this problem. Immunologists separate the molecularly non-specific response of the “innate” immune system, which includes everything from scratching to the recognition of protein motifs characteristic of bacteria, and the molecularly specific “adaptive” response, by which specialized cells recognize evolving features of never encountered before pathogens. From another angle we can consider different ways of destroying a pathogen: either swallowing pathogens whole, which is done by cells of the innate immune system called neutrophils and macrophages; or killing our own cells that have been infected by a pathogen or are tumourous — as performed by representatives of the adaptive immune system called killer T-cells. Alternatively, the adaptive immune system produces specialized molecules called antibodies that smother the invader: they attach to pathogens to prevent them from entering cells and multiplying; they bind to bacterial toxins, thereby disarming them; or they bind directly to bacterial cells, flagging them for consumption by macrophages.

That short overview gives us an idea of the many strategies that both pathogens and the host organism have at their disposal (Fig. 1). Pathogenic cells (a bacteria, virus, or a tumour cell) are programmed to divide. The immune system is there to prevent this. To a large extend, its main challenge is to recognize the unknown. Pathogens are constantly evolving to escape recognition by the immune system. Although the host organism does have a certain list of “warning sign features,” most of which are taken care of by the innate immune system, it is up against a large set of constantly moving targets — as illustrated in our everyday experience by the evolving influenza virus, which requires a new vaccination every year. For this reason, the strategies developed by the immune system are mostly statistical, and require multiple interactions between different types of cells, a lot of checks and balances that leads to a multiplicity of time and length scales (Fig. 2). Despite this complexity, the immune system works remarkably reliably. How do these interactions on many scales dynamically come together in a self-organized way to build a complex sensory system against a high dimensional moving target? This review breaks this high level question into smaller problems and presents some results and concepts contributed by physicists. It also presents a summary of the current methods, experimental and theoretical, used in the field. We do not shy away from the biology of the immune system, but to help the physicist navigate the complexity of immunology, biological details will be introduced as we go along. For an introduction of the immune system for the non-specialist, we refer the reader to the short but excellent book by Lauren Sompayrac Sompayrac1999 . This review does not aim to be exhaustive, but rather focuses on the important physical concepts of immune function, reducing biological complexity to a minimum whenever appropriate.

From an evolutionary perspective, all organisms have some form of protection. Bacteria protect themselves from viruses through specific CRISPR (Clustered Regularly Interspaced Short Palindromic Repeats)-Cas, or unspecific restriction modification systems. The innate immune system is shared by many animals, and is largely similar between us and flies. The adaptive immune system evolved in jawed vertebrates and also has changed little between fish and mammals. Plants also have a well-studied innate system, and have recently been shown to have elements of adaptive immunity. Immunity is therefore a basic element of life. However the details of its implementation and spatial organisation scale with the organism.

Burnet’s clonal selection theory Burnet1957 , built on previous observations and ideas Ehrlich1900 ; Jerne1955 ; Hodgkin2007 , provides a theory of the adaptive immune system in the same sense that Darwin’s theory provides a theory or framework for evolution. Burnet’s theory states that molecules of the pathogens stimulate specific B and T lymphocytes (cells of the adaptive immune system) among a pre-existing ensemble of possible cells, which leads to proliferation of this specifically selected clone. Molecules thus recognized are called antigens. It explains the diversity and specificity of the adaptive immune system, also highlighting its adaptive nature. The framework is often summarized in four assumptions: (i) each T and B lymphocyte cell has one type of receptor; (ii) receptor-antigen binding is required for cell or receptor proliferation; (iii) offspring of the stimulated cell have the same receptor as their parents; (iv) cells that have receptors that recognize the organism’s self molecules (self-antigens) are removed early in their development. The theory was validated by showing that B-cells always produce one receptor NossalLederberg and later by experiments showing immunological tolerance to factors introduced in the embryonic period or immediately after Medawar1953 ; Medawar1956 . Burnet’s theory provides an incredibly successful framework to understand the adaptive immune system, yet it is not a quantifiable theory that can be tested against concrete measurements. While what remains to be “filled in” may seem like details, these details still provide a lot of puzzles and uncertainties, also at the conceptual scale. In this review, we present certain examples of quantified ideas that have recently emerged.

We also present some (albeit not all) of the quantitative puzzles, starting from the smallest molecular (sections II and III) and cellular (sections IV and V) scales and moving to the organismal (sections  VI, VII and  VIII) and population wide (section IX) scales. How do cells discriminate between self and non-self in such an exquisite way integrating processes taking part on many timescales (section III on antigen discrimination)? How do cells communicate and coordinate to orchestrate the immune response as a collective phenomenon (section  IV on cytokines)? How do cells adopt specialized phenotypes through cell differentiation to span the range of functions they must fulfil (section V on cell fate)? How does the system cover the space of possible specificities allowing for complete protection against the unknown (section  VI on repertoires)? How does the immune system as a whole adapt to the changing environment (section VIII) and how does it influence virus evolution (section IX)?

This review covers a lot of topics, many of them connected. We try to mention these connections, however we have tried to make the sections stand alone and the reader does not have to (and probably should not attempt to) read them all at once, or in the presented order. At the end we include a glossary of the biological terms (section XI) to help navigate the immunological terminology. Detailed appendices (Sec. XII) present general methods that are used in physical immunology as well as other fields.

II Physical chemistry of ligand-receptor interaction: specificity, sensitivity, kinetics.

Our exploration of the immune system starts at the molecular scale, through the binding of ligands and receptors expressed at the surface of immune cells. Ligands carry information about the pathogenic threat and can be of two types: antigens, i.e. bits of proteins that are recognized by the immune system; and cytokines, which are small molecules secreted by immune cells to communicate with each other about their current state and experience. Binding events between these ligands and their cognate receptors provide the raw signals that cells must interpret to adapt their individual and collective behavior, eventually mounting an immune response in the case of recognition.

All immune decisions start with the classical physical chemistry of ligand–receptor interactions:

[TABLE]

This extremely simple reaction leads to strong limits on how fast signals can spread in the immune system (and in biological systems in general), which in turn allows us to discriminate possible regulatory scenarios. This constraint stems from physical-chemical limits on the parameters of these reactions, typically in the picomolar to millimolar range for the equilibrium dissociation constant $K_{D}=k_{\rm off}/k_{\rm on}$ (we recall that $1nM=10^{-9}{\rm mol}/L=10^{-24}\mathcal{N}_{\rm Avogadro}\mu m^{-3}\approx 0.6$ molecules per $\mu m^{3}$, with $\mathcal{N}_{\rm Avogadro}\approx 6\cdot 10^{23}{\rm mol}^{-1}$). In this section, we discuss the physical considerations leading to estimates of the biochemical rates driving in these reactions. Such a discussion of the physical chemistry of ligand-receptor interactions is necessary to understand key quantitative aspects of immune activation, as well as discrimination between self and non-self antigens by B & T cells (as we will discuss in Sec. III).

II.1 Diffusion-limited reaction rates

The basic laws of physical kinetics Landau_Physkin can be used to obtain an estimate of the diffusion-limited rate of molecular association (see Sec. XII.1.1 for a derivation): $k_{\rm on}\leq k_{\rm diffusion}$, with

[TABLE]

where $R_{\rm Ligand}$ and $R_{\rm Receptor}$ are the radii of the ligand and of receptor binding pockets (both assumed spherical), while $D_{\rm Ligand}$ and $D_{\rm Receptor}$ are their diffusion coefficients. Since the receptor in Eq. 1 is usually embedded in the lipid bilayer of an immune cell, we can assume it is relatively immobile compared to its ligand $L$, which diffuses more rapidly in the extracellular medium or in the cytoplasm, $D_{\rm Receptor}=D_{\rm membrane}\ll D_{\rm Ligand}=D_{\rm solution}$.

The ligand’s diffusion coefficient can be estimated from Stokes-Einstein’s formula:

[TABLE]

for a small ligand of size $R_{\rm Ligand}=1nm$, diffusing in extracellular medium (whose viscosity is given by that of water, $\eta=0.7\,mPa\cdot s$), at body temperature $T=310K$ with Boltzmann’s constant $k_{B}=1.38\cdot 10^{-23}$J/K). Then, from Eq. 2, assuming a small target ($R_{\rm Receptor}\ll R_{\rm Ligand}$), we obtain:

[TABLE]

Eq. 4 gives a general upper bound for any ligand-receptor association rate, whether it is immunological or not. This limit is conceptually equivalent to the speed of light (although clearly not as fundamental).

Yet, there are examples when ligand-receptor associations between two small proteins seem to “beat” this diffusion limit: $k_{\rm on}>k_{\rm diffusion}$ Hager2009 . Such an apparent paradox can be resolved when considering the limiting step for this association. If the ligand interacts weakly with a large object, (e.g. the entire plasma membrane of a cell, DNA coils), inducing directed diffusion in a constrained space before reaching its specific targets Halford2004 ; Gorman2008 , this pre-equilibration step implies that the cross-subsection of the object our ligand needs to hit can be much larger. The diffusion molecules then can “hop” around these non-specific binding sites to accelerate their search for the specific binding sites. In these cases, the effective $R_{\rm Receptor}$ of collision is the macroscopic scale associated with the large object of weak/non-specific interactions ($5\mu m$ for cells, $>0.1\mu m$ for DNA coils), hence, an apparent increase in the rate of collision according to Eq. 2. A careful calculation (see e.g. Slutsky2004 ) yields the appropriate bound.

The diffusion coefficient will be slowed down for larger macromolecular complexes, and for molecules diffusing inside the cells (there, the apparent viscosity can be increased by 10-fold). Additionally, when considering biomolecules embedded within cell membranes, one must take into consideration the dramatically increased viscosity that leads to slowed-down diffusion for receptor proteins on the surface of immune cells ($D_{\rm solution}\approx 200\mu m^{2}/s$ vs $D_{\rm membrane}\approx 5\mu m^{2}/s$ for a typical protein of radius 2nm). For this reason, when considering a soluble protein interacting with proteins embedded within the plasma membrane of immune cells (*e.g. *extracellular cytokines being captured by a cytokine receptor, intracellular enzymes –such as kinases or phosphatases, that catalyze the phosphorylation/dephosphorylation of proteins– interacting with an activated receptor, etc.), the diffusion of the membrane proteins is so small ($D_{\rm membrane}\ll D_{\rm solution}$), than it can be neglected and

[TABLE]

II.2 Extrapolating collision rates in solution to association rates in the physiological context

While the rate of collision calculated above is an upper bound to the association rate, a very important limitation must be taken into account: not every molecular collision will lead to their association, and we must estimate the probability of a successful association event. This probability is given by Arrhenius’s law:

[TABLE]

where $\Delta G_{\rm association}$ is the free energy barrier (entropic and enthalpic) that molecules must overcome to associate.

The estimation of $\Delta G_{\rm association}$ is tricky as it requires a deep structural understanding at the atomic level of the entropic, energetic and conformational changes associated with bond formation between two large biomolecules. For the entropic contribution, one rule of thumb is to estimate the numbers of degrees of freedom constrained by the association between biomolecules: simply aligning two biomolecules (in rotation or in translation) incurs an entropic cost of at least $6k_{B}T$, reducing the probability of association of each collision by $e^{6}\sim 400$. Hence, the basal association rate (before taking into account more subtle molecular constraints) is $\sim 5\cdot 10^{6}M^{-1}s^{-1}$ rather than the $\sim 2\cdot 10^{9}M^{-1}s^{-1}$ estimated in Eq. 4. Additional corrections should be made for each pair of biomolecules (as we will discuss in Sec. II.3.4).

II.3 Rates and numbers in the physiological context

II.3.1 Association rates

The formula in Eq. 2 is particularly useful in the context of quantitative immunology, when one must estimate the kinetics of molecular interactions in different context (solution, intracellular cytoplasm, surface plasma membrane etc.). Most kinetic parameters for biomolecular interactions are measured in solution, as experimentalists routinely purify ligand–receptor pairs and test their interactions in their soluble form (e.g. by calorimetry or by surface plasmon resonance – in the latter case, one molecule must be immobilized). Such measurements can be used to estimate the hard-to-predict activation barrier $\Delta G_{\rm association}$ for molecular association by inverting Eq. 6, using the measured in vitro $k_{\rm on}$, and the estimated $k_{\rm diffusion}$. Eq. 2 can be used to rescale the viscosity $\eta$ and cross-section $R_{\rm Ligand}+R_{\rm Receptor}$ of molecular interactions in such a way that one can translate soluble measurements into physiologically relevant parameters.

II.3.2 Dissociation rates

As discussed below, immunological interactions span the range of very short lived interactions, $k_{\rm off}>10\,s^{-1}$ or $\tau_{\rm off}=(k_{\rm off})^{-1}<0.1\,s$ for antibody binding to its target in the initial phase of an immune response, to extremely long-lived interactions, $k_{\rm off}<10^{-4}s^{-1}$ or $\tau_{\rm off}>3$ hours for cytokines interacting with their receptors. These estimates set a huge range of time scales the immune system must deal with, even before taking into consideration delays in cellular responses and how these affect the ligand environment on time scales ranging from hours to days. These considerations form the crux of the matter for quantitative immunology at the cellular scale: while the physical chemistry of ligand-receptor interactions is straightforward, the immune system builds a response of devilish complexity from such elementary interactions.

II.3.3 Numbers of receptors per cell

In some situations, such as in T cell antigen recognition, where both the T cell receptor and the antigen exist in a membrane-bound form. All dynamics of interactions must be estimated taking into account the surface concentrations of molecules, with proper adjustment for slower diffusion in the association rates.

One often measures the number of receptors per cell $\#R_{\rm cell}$, e.g. by assessing the number of cytokine receptors on the surface of cells using quantitative flow cytometry, while the ligand is provided in soluble forms, as is the case with cytokines. In this case, one must estimate the solution-level concentrations of available receptors $[R]_{\rm total}$ in the reaction volume $V$:

[TABLE]

where $N_{\rm cell}$ is the number of cells in volume $V$. For example, the effective “concentration” of receptors binding the IL-7 cytokine within a lymph node can be estimated, knowing that each of $10^{7}$ T cells within this 50 $\mu$L volume express $10^{3}$ receptors:

[TABLE]

while $K_{D}=10^{-11}M=10pM$ for IL-7 binding to its receptor. Hence $[{\rm IL-7R}]_{\rm total}\gg K_{D}$ and any secreted IL-7 will rapidly be captured by the cells within a lymph node. The high density of receptors and cells within lymphoid tissues makes for an interesting regime of competition for soluble ligands, as discussed below (Sec. IV). Such normalization by the Avogadro number, while straightforward, is of critical significance in many immunological configurations, as one must bridge the molecular scale of immune agents (cytokine secretion and consumption) with the functional scale of the system, where competition between cells for soluble ligands takes place.

II.3.4 Typical binding rates and some biology

We summed up above the basic and general physical chemistry that drives the molecular association and dissociation of molecules in order to delineate key quantitative parameters in immunological regulation. The precise values of these rates is in fact crucial for the immune system to recognized pathogens.

At its core, the immune system can be considered as a collection of cells, called leukocytes, whose activation registers the presence of “new” molecules — lipid and nucleic acid signatures of viral and bacterial infection for the innate system, or unknown proteins and peptides for the adaptive immune system — and translates into a defensive response — secretion of neutralizing antibodies, killing or phagocytosing of infected targets, etc.

In the innate immune system, non-self recognition is encoded in the structure of the ligand and the receptor. Monocytes, macrophages, dendritic cells, and Mast cells are endowed with genome-encoded receptors called Toll-Like Receptors whose ligands are non-mutable elements of pathogens — peptidoglycans and liposaccharides from the membrane of bacteria, CpG unmethylated dinucleotides derived from viruses. Recognition of non-self in that case is a “simple” lock-and-key proposition whereby pathogenic ligands engage the receptors and trigger a signaling response that activate innate responses. Hence, for the innate immune system, self vs. non-self discrimination is hard-wired at the molecular level and gets triggered in the early moments of an immune response.

The adaptive immune system offers a more challenging issue. Each B and T cell clone expresses its own unique receptor (and one only), whose assembly is driven by random events (described in detail in Sec. VI and Sec. VIII). The ligands of these random receptors, which are called antigens, are not pre-determined — in fact, they may not even exist at the time of birth of the organism, if for instance they are derived from fast evolving strains of viruses. Understanding how the binding and unbinding rates of such ligand-receptor pairs contribute to self vs. non-self discrimination is one of the core issues in quantitative immunology. Here we give some orders of magnitude of the binding rates to be discriminated.

B and T cell differ fundamentally by the type of antigenic ligands they interact with. T-cell receptors (TCR) interact with their antigen on the surface of antigen-presenting cells, which they scan for anomalies. These antigens are complexes which include a short peptide (around 10 amino acid-long) produced within the antigen-presenting cells by chopping up larger proteins expressed by the cell, either from the genome or from foreign pathogens. Each peptide is loaded onto a large protein called the Major Histocompatibility Complex (MHC), forming a peptide-MHC (pMHC) complex. By contrast, B cell receptors (BCR) bind directly to proteins. The exact position where the binding occurs is called an epitope. Antibodies are soluble versions of the BCR with the same antigenic specificity, and also bind directly to the pathogen proteins to neutralize them.

The affinity of antibodies (produced by B-cells) are highly variable. At the onset of an immune responses, B-cells express and secrete Immunoglobulin M (IgM), which is a pentamerized version of antibodies whose affinity for the target is weak ($K_{D}>10\mu M$): the ability of secreted IgM to bind to pathogens and guide them towards eradication is then driven by the multimeric interaction of IgM with its target. Upon engagement of the adaptive immune response, the interplay between T-cell help and B-cell somatic hypermutation drives the maturation of antibody affinity (see Sec. VIII). At this point, B cells increase the affinity of the antibody they express (down to $K_{D}<1nM$) and switch the class of antibodies from IgM to IgG i.e. from a pentameric version to a monomeric version that does not require as much multimerization of the antibody to bind to their target. Such decrease in $K_{D}$ and increase in affinity is not driven by better association rates: IgM binds to their target with typical $k_{\rm on}$ around $10^{6}M^{-1}s^{-1}$ – similar to IgG. $k_{\rm on}$ is essentially driven by the collision rate and the contribution of the energetic barrier in the association rate is very limited. Instead, the improvement afforded by affinity maturation is driven by a better fit between antibody and antigen, resulting in a lower $k_{\rm off}$.

In the case of T-cell antigen discrimination, the difference between ligands that will or not activate the immune response is extremely sensitive. Some peptides can elicit a robust activation with a single pMHC molecule, while mutated versions of this peptide differing by just a single amino-acid are unable to trigger T cells even in large quantities ($>10^{6}$). The striking feature in T-cell biology is that a single mutation in the peptide of the pMHC only has a marginal effect on $k_{\rm off}$, while greatly impacting its functional capacity to activate T-cells Altan-Bonnet2005 ; Feinerman2008b ). Experimentalists can measure the biophysical characteristics of pMHC–TCR interactions in solution, either by surface–plasmon resonance on soluble pMHC interacting with surface–immobilized TCR, by directly observing the surface of cells using single-molecule imaging and fluorescence energy transfer, or by measuring mechanical forces (see Sec. III.1.3). Typically, an agonist ligand, defined as a pMHC that will elicit an immune response, binds to their TCR with $\tau_{\rm off}\sim 1-10$s, while non-agonist ligands bind with timescales that are at least 3-fold shorter (from 0.3 to 3 s), although this varies according to the particular TCR, and also depends on the MHC class (of which there are two, as we will see later). This relatively small difference in binding results in a large fold change ($\geq 10^{5}$) in response, regardless of the concentrations. Understanding how cells can perform such sensitive discrimination is a major challenge of quantitative immunology, which we will discuss in detail in Sec. III.

II.4 Receptor-antigen specificity

T- and B-cell receptors do not interact as set of locks and keys matched to each other. Instead, each receptor can bind many different possible ligands, and vice-versa. Here we discuss numbers, data, and models that characterize this many-to-many mapping.

II.4.1 Cross-reactivity

Despite the huge diversity of BCR and TCR (reviewed in Sec. VI), the diversity of antigens may be even greater. While such an estimate is difficult for BCR epitopes, a rough estimate of the number of antigens can be obtained for peptide-MHC complexes. For each of the two MHC classes that exist, only a few ($\sim 5\%$) percents of peptides can be presented on any of the 6 MHC genes expressed in an individual. For peptides of length 12, this amounts to $0.05\times 20^{12}\approx 2\cdot 10^{14}$ antigens. While this might be a manageable number for humans, who harbor $\sim 10^{12}$ B and T cells, this is not by mice, who have $\sim 10^{8}$ such cells, and this could be an issue even for humans for longer peptides. This argument led Mason Mason1998 to conclude that each TCR and BCR must be able to recognize more than one antigen, a phenomenon called cross-reactivity or polyspecificity. Cross-reactivity has a less discussed counterpart, which is that each antigen must be recognized by a large variety of receptors. If there are $N$ antigens, each of which can be bound by $k$ receptors, and $R$ distinct antigen receptors, each of which can recognize $\ell$ antigens, then one must have $Nk=R\ell$. To fix ideas, the ratio $p=\ell/N=k/R$, the probability that a random antigen and receptor bind together, is thought to be around $10^{-5}$ Boer1993 .

It should be emphasized that finding the sequences of binding pairs of antigens and lymphocyte receptors remains an essentially experimental question, which has received renewed attention lately thanks to the reduced costs of sequencing allowing for high-throughput binding or functional assays Birnbaum2014 ; Adams2016 ; Glanville2017 ; Dash2017 . However, despite increasing amounts of data on these binding pairs compiled in useful databases Shugay2018 ; Tickotsky2017 , and current attempts to leverage these data to make prediction using modern machine learning techniques Jurtz2018 ; Sidhom2018 ; Jokinen2019 , there exists no good predictive model of antigen-receptor specificity. Therefore, the binding models presented in the next paragraph should be viewed as useful toy models for investigating broad properties of cross-reactivity.

II.4.2 Models of receptor-antigen binding

The recognition process of antigens by immune receptors is based on molecular interactions between the two proteins: the receptor protein and the antigen, which can be summarized by $K_{D}=c_{0}e^{\Delta G(a,r)/k_{B}T}$, where $\Delta G(a,r)$ is the interaction free energy of binding between antigen $a$ and receptor $r$, and $c_{0}$ a constant. This energy may be modeled as a string matching problem Chao2005 . In such toy models, each interaction partner is decribed by a string of length $N$ representing the interacting amino acids. The binding energy between these two specific interaction partners is assumed to depend additively on the interaction between pairs of facing amino-acids:

[TABLE]

where the interaction matrix, ${J_{k}}(a_{k},r_{k})$ is an $q\times q$, where $q$ is the size of the space that defines the variabilty of residues at each position. If we describe each residue as one of the $q=20$ amino acids, we need to define an amino acid interaction matrix. The Miyazawa-Jernigan matrix Miyazawa1996 was used in such a model to suggest that thymic selection favors moderately interacting residues in TCR Kosmrlj2008 , or more recently to study the effect of thymic selection on tumor antigenic peptides George2017 . The maximum affinity of a receptor to a large set of antigens, which plays a key role in thymic selection (Sec. VI.3), is amenable to a statistical mechanis analysis of extreme value statistics Kosmrlj2009 ; Butler2013a , leading to universal features that are robust to the details of the models.

Alternatively, reduced models have been considered, where each residue is defined using a projection that attempts to capture the main biophysical and biochemical properties in a generalized “shape space” Perelson1979 . In the initial string model Detours1999a , $a_{k}$ and $r_{k}$ were bounded natural integers, and $J_{k}=a_{k}\oplus r_{k}$, where $\oplus$ is the exlusive OR operator acting on each digit of the binary representations of $a_{k}$ and $r_{k}$. This choice was motivated by algorithmic ease rather than biophysical realism. These models were used by Perelson and collaborators to investigate the effect of thymic selection Detours1999a ; Detours2000 or the immune response Chao2005 .

A more drastic reduction is to binarize the antigen: for each epitope position, the amino-acid is defined either as the one present in the viral wild type ($a_{k}=1$) or another one (a mutant, $a_{k}=-1$) Wang2015a . The selective pressure exerted by each receptor, which acts as a “field” on $a_{k}$, was drawn as a number $r_{k}$ from a random continuous distribution, which reduces the model to $J_{k}(a_{k},r_{k})=r_{k}a_{k}$. In Wang2015a , it was additionally assumed that certain positions of the viral epitope were constrained to take the wildtype value, $a_{k}=1$ for $k>M$, because of strong conservation at these sites. A similar description in terms of binary strings Nourmohammad2015 assumes both the viral epitope and BCR to be binary strings ($a_{k}=\pm 1,\ r_{k}=\pm 1$), with a fixed interaction strength, $J_{k}(a_{k},r_{k})=\kappa_{k}r_{k}a_{k}$, with again conserved viral positions $k>M$ for which $a_{k}=1$ is imposed. Within this description a mismatch between the receptor and the viral “spins” induces an energy penalty. The convention is such that smaller energies imply stronger binding and better recognition. Similar models have been used recently for co-evolution of BCR and HIV Wang2015a ; Nourmohammad2015 , the results of which we describe in detail in subsection VIII.

II.4.3 Data-driven receptor-antigen binding models

To go beyond the toy models described above, one must experimentally map out the binding energies between specific pairs of antigens and receptors. This can be done in a massively parallel way using high-throughput experiments assaying the binding affinity of many pairs in single experiments. Such an approach was applied to a deep mutational scan experiment reporting the dissociation constant, $K_{D}$, of a large number of antibody variants against a fixed antigenic target, fluorescin Adams2016 .

The simplest model assumes an additive contribution of each residue to the binding free energy $E(a,r)=\ln(K_{D}/c_{0})$ as in (9), but with $J_{k}(a_{k},r_{k})=h_{k}(r_{k})$ fixed to a single value of the antigen $a$. However, statistical analysis shows that such an additive model is not able to capture all the variability in the binding energy, accounting for less than $2/3$ of the variability in double and triple mutants. Epistasis, defined as non-additive effects, accounts for $25-35\%$ of variability in the binding energy between antibodies and the antigen and a large fraction of the epistasis was found to be beneficial Adams2019 .

Non-additivity in the binding energy is likely to be a general feature of both TCR and BCR. More sophisticated models including intra-protein or higher-order interactions will be needed to accurately predict receptor-antigen affinity. Ultimately, it would be interesting to reconcile the useful picture of an effective binding shape space with affinity landscapes inferred from data. A major roadblock is that most current experimental techniques only allow for varying one element of the receptor-antigen pair, either testing a library of receptors against a fixed antigen, or a library of antigens against a fixed receptor Birnbaum2014 . Full characterization of the binding landscape would require new methods to test double libraries of antigens and receptors in an ultra-high-throughput manner.

II.4.4 Modeling immunogeneticity

Rather than modeling the full details of the TCR-antigen interaction, an alternative strategy is to focus on the immunogenic potential of particular antigens, without explicitly modeling the TCR. Such an approach was followed to predict the response of cancer patients to immunotherapy. The prediction is based on the knowledge of the pMHC antigens presented by the tumor cells, called neo-antigens Luksza2017 . To estimate the ability of a neoantigen $a$ to be recognized by the TCR repertoire, a score was calculated to evaluate its similarity to a database $D$ of antigens known to elicit an immune response:

[TABLE]

where $s_{a,a^{\prime}}$ is the alignment score between the neoantigen $a$ and antigens $a^{\prime}$ from the database, and $k$ an adjustable parameter. The immunogenicity of $a$ is then predicted to be:

[TABLE]

This quantity was combined with the likelihood $A(a)$ that neoantigen $a$ is presented by class-I MHC (predicted by a neural network model Andreatta2015 ) to form a global score $-I(a)A(a)$. With its minus sign, this score reflects the fitness of the tumor cells carrying antigen $a$. These scores were found to be predictive of patient survival after immunotherapy.

While this similarity-based approach holds great promise for generally predicting the immunogeneticity of antigens, further tests are needed to validate the method for broader use. The link between survival and TCR recognition is indirect, and is complicated by the fact that many neoantigens and tumor clones are involved. Direct functional tests of the immune response against libraries of antigens would help consolidate the foundations of this approach.

III Antigen discrimination

We now move to the cellular scale, and describe how the signal propagates in minutes from antigen-receptor binding to biochemical networks, allowing for fine antigen discrimination. We leave to a later section the discussion how this signal is later integrated to make decisions about cell fate and response on timescales spanning hours to days (Sec. V).

III.1 T cells

III.1.1 Kinetic proofreading for ligand discrimination

To model T-cell antigen discrimination, we remain within the quantitative parameters of the “lifetime dogma,” in which the lifetime $\tau_{\rm off}$ of TCR-pMHC is the sole determinant of the discrimination process. While this dogma was derived from the biophysical measurements on TCR/pMHC interactions Feinerman2008b , as any dogma it must be revisited regularly for exceptions and refinements. For example, recent studies in the mechanics of TCR triggering have highlighted a new mode of signal initiation, as we shall discuss in Sec. III.1.3.

Upon engaging their pMHC ligands, T-cells trigger a cascade of phosphorylation events, which consist of adding phosphoryl groups to conserved residues of TCR-associated chains. The phosphorylations in turn trigger other reactions, ultimately activating transcription factors that modulate gene regulatory responses.

Decision making relative to antigen discrimination can thus be modelled at the phenomenological level by the state of phosphorylation. We assume that, when the TCR complex has accumulated a certain number of phosphoryl groups $N_{\rm P}$, it flips a digital switch (to be specific, the activation of the NFAT or ERK molecules Altan-Bonnet2005 ) which defines the onset of T-cell activition.

One of the first quantitative models that tackled the exquisite specificity of antigen discrimination during T cell activation was proposed by McKeithan in 1995 McKeithan1995 , building upon the classical kinetic proofreading (KPR) scheme proposed by Hopfield and Ninio in the context of protein translation Hopfield1974 ; Ninio1975 . Upon binding its ligand, TCR progresses through the different steps of the cascade, controlled by a phosphorylation rate $k_{\rm P}$ and a de-phosphorylation rate $k_{\rm deP}$ (Figure 3: A). A key aspect of the kinetic proofreading is that, at each step, it is assumed that the complex is de-phosphorylated into unoccupied receptors upon ligand unbinding: this is a reasonable assumption because the CD45 phosphatase — the enzyme that removes phosphoryl groups from the TCR complex — prevents rebinding of TCR by the antigen due to its large size, thus ensuring complete de-phosphorylation before the next binding event. A consequence of this resetting upon unbinding is that antigens forming short-lived complexes with the TCR will typically only accumulate a small number of phosphoryl groups, below the threshold $N_{\rm P}$ necessary to activate the immune response, while long-lived complexes can progress more deeply through the phosphorylation cascade.

Neglecting dephosphorylation ($k_{\rm deP}\ll k_{\rm P},k_{\rm off}$), the number of pMHC–TCR complexes $C_{i}$ that accumulate $i$ phosphorylations scales with $i$ as McKeithan1995 :

[TABLE]

At steady state, the total number of complexes, $C_{\rm total}=C_{0}\sum_{\rm i=1}^{N_{p}}C_{i}$, is given by a second-order equation describing the balance between the rate of binding between free antigens and free receptors, $k_{\rm on}/(V.\mathcal{N}_{A})(L-C_{\rm total})(T-C_{\rm total})$, where $L$ is the total number of antigens, $T$ the total number of TCR, $V$ the volume, $\mathcal{N}_{A}$ –the Avogadro number, and the rate of unbinding events, $k_{\rm off}C_{\rm total}$. In the limit when TCR are not limiting, $C_{\rm total}\ll T$, this balance yields the total number of complexes Francois2013 :

[TABLE]

Further assuming a relatively slow phosphorylation rate, $k_{\rm P}\ll k_{\rm off}$, the number of fully phosphorylated complexes $C_{N_{\rm P}}$ scales as $\propto k_{\rm on}Lk_{\rm off}^{-N_{\rm P}}=Lk_{\rm on}\tau_{\rm off}^{N_{\rm P}}$ as a function of the antigen characteristics — concentration $L$ and binding affinity $\tau_{\rm off}$. This scaling shows how alterations in the pMHC ligand impacting $\tau_{\rm off}$ will get amplified into large differences in the amount of phosphorylation accumulating on the TCR:

[TABLE]

In other words, kinetic proofreading “reads off” the lifetime of the pMHC–TCR complex and amplifies differences in the output. As with the original KPR scheme of Hopfield and Ninio, this amplification implies energy expenditures caused by the phosphorylation steps. Structurally speaking, the TCR complex contains 20 phosphorylation sites (specifically, tyrosine residues), which can potentially participate in the kinetic proofreading cascade. This means that $N_{\rm P}$ can be as large as 20, and a two-fold change in pMHC–TCR lifetime could be amplified into a $2^{20}\approx 10^{6}$ fold difference in the number of phosphorylated receptors $C_{\rm N_{p}}$. This mechanism can account for the fact that TCR ligands with minute differences in affinity may elicit very different signals. In addition, it is consistent with the lifetime dogma, in that $\tau_{\rm off}$ has a much more determining impact on activation than $k_{\rm on}$.

However, KPR is insufficient to capture all aspects of T cells’ ability to discriminate between structurally-related ligands. As pointed out by Altan-Bonnet & Germain AltanBonnet2005 , T cells not only achieve high specificity in ligand discrimination, but they also maintain high sensitivity, as they are be able to trigger activation from a single agonist ligand Sykulev1996 ; Irvine2002 ;. In addition, they respond very fast, typically within minutes of encountering an antigen presenting cell. These requirements for speed, sensitivity and specificity are hard to achieve all at once. For example, a large number of phosphorylation steps in the KPR scheme allows for high discriminability, but also implies a slower response, as each step must be slow enough to discriminate between agonists and non-agonists. It also affects sensitivity, as more steps imply a lower chance of making it the activation step. These considerations demonstrate how quantitative modeling invalidates a bare KPR scheme, and calls for additional mechanisms.

III.1.2 Adaptive kinetic proofreading

To fulfill the conflicting requirements of specificity, sensitivity, and speed, one must expand on the simple kinetic proofreading scheme by adding a mechanism of adaptation to modulate the proofreading steps.

The resulting adaptive kinetic proofreading (AKPR) model Altan-Bonnet2005 ; Francois2013 , which is based on experimental evidence presented by Štefanova et al. Stefanova2003 , quantitatively accounts for key features of ligand discrimination by the TCR.

The adaptation module is introduced through a negative feedback mediated by a phosphatase $S$ (which removes phosphoryl groups), which is itself activated by the engaged TCR in state $C_{1}$, so that its steady-state concentration is $S=C_{1}/(C_{1}+K_{S})$, where $K_{S}$ is a model parameter. The dephosphorylation rate is then enhanced in proportion to the phosphatase concentration, to $k_{\rm deP}+\gamma S$. The equations for steady state have a simple closed form, whose solution is given by the root of a fourth-order polynomial.

Fig. 3 presents a simple graphical argument illustrating how KPR and AKPR schemes perform in their task of discriminating ligands. In the adaptive kinetic proofreading scheme, the negative feedback mediated by $S$ enforces a higher selectivity of antigen discrimination without affecting the high sensitivity of the response.

To better understand how this is possible, consider a large concentration $L$ of non-agonists. Their unspecific engagement with the TCR will cause the first phosphorylation of the complex into $C_{1}$, in proportion to their concentration $L$. This $C_{1}$ state activates the phosphatase $S$, which in turn increases the specificity of the TCR by accelerating its dephosphorylation. When $L$ is large, the amount of phosphatases, and hence desphosphorylation rate, will exactly balance out the amount of non-agonist-bound TCR complexes, preventing these complexes from reaching the activation state $C_{\rm N_{P}}$ (see Figure 3 B). As a result, for non-agonists $C_{N}$ is essentially flat and stays below the threshold for T cell activation, regardless of $L$.

The activation threshold for $C_{N_{\rm P}}$ must to be set at very low values to account for the extreme sensitivity of the TCR signalling cascade (a single ligand is sufficient to trigger activation). Hence, one must introduce stochastic equations to appropriately tackle the low number of molecules involved in the TCR decision making, although this does not affect the qualitative features of this particular model.

To conclude, the adaptive kinetic proofreading scheme satisfies the joint requirements of ligand discrimination and sensitivity. It can explain how a single agonist ligand can trigger T cell activation, while a large number ($>10^{5}$) of non-agonist ligands cannot. Additionally, it can account for the speed of ligand discrimination by T leukocytes, as only two steps can be sufficient to sort ligands Lalanne2013 . It is important to emphasize again the importance of a quantitative approach and of modeling to appreciate why achieving speed, sensitivity and specificity in ligand discrimination is indeed such an amazing feat of the immune system.

III.1.3 Coupling mechanics and biochemistry: the significance of forces for ligand discrimination

One of the most striking events in the T-cell activation process is the dynamic membrane reorganization into what has been called an immunological synapse. Within minutes of signal initiation, the antigen receptors on the surface of lymphocytes congregate at the center of the interface with their antigen-presenting cells (this has been dubbed a central Supramolecular Assembly Complex or c-SMAC), while transmembrane receptors and phosphatases needed for signal propagation accumulate at the periphery of the contact area (p-SMAC and d-SMAC respectively, standing for peripheral and distal). Such a bullseye pattern has been studied quantitatively using high-resolution and single-cell time-lapsed microscopy. The functional significance of immunological synapses has been identified for the long-term activation of T cells, such as directed killing of antigen-presenting cells or directed secretion of cytokines Huse2006 . As the cell biology of immunological synapse formation was being studied, the question arose as to whether such cellular structure could contribute to self vs non–self discrimination. Qi et al. Qi2001 introduced a Ginzburg-Landau model for the mechanical deformation of the T-cell membrane during synapse formation. The model assumes the presence of several receptor-ligand pairs, indexed by $i$ (e.g. the TCR-pMHC pairs is one these pairs). The free energy $\mathcal{F}$ for the membrane system is estimated as:

[TABLE]

and the system evolves according to:

[TABLE]

and

[TABLE]

where $z(x,y)$ is the membrane coordinate for the T-cell membrane, $R_{i}$, $L_{i}$, and $C_{i}$, are the surface densities of receptors, ligands, and receptor-ligand complexes of the $i^{\rm th}$ type, $z_{i}^{(0)}$ is the length of the $C_{i}=R_{i}$-$L_{i}$ complex bond at rest; and the $\zeta$’s are thermal noises for each process (their distributions are not explicitely stated in the original publications Qi2001 ). The free energy $\mathcal{F}$ sums up the elastic energies of each type of complex $i$ (with stiffness $\lambda_{i}$), and the elastic energy of the membrane. $k_{i}^{\rm on}$ and $k_{i}^{\rm off}$ are the binding and unbinding rates of ligand-receptor complex formation.

This model is purely mechanical and passive, in the sense that no energy is injected nor is there any biochemical feedback. Yet it captures salient features of the spontaneous formation of immunological synapses. Importantly, synapse formation strongly depends on the lifetime of ligand-receptor complexes, and thus allows for fine discrimination between agonist and non-agonist ligands. The model also emphasizes the functional relevance of size differences in ligand-receptor pairs. For instance, TCR-pMHC complexes have $z^{(0)}_{i}=15$ nm, while adhesion complexes such as ICAM-1/LFA-1 are much larger, $z^{(0)}_{i}=42$ nm, as it was confirmed experimentally Choudhuri2005 ; Choudhuri2007 .

This model makes several predictions that were validated experimentally. First, it explains how the formation of the characteristic bullseye pattern emerges from the passive sorting of molecules based on their size, upon coupling with the mechanical deformation of the surface membrane. Second, the model predicts an optimum for the lifetime of the pMHC-TCR complex for which the accumulation of antigen in the center of the synapse is maximal. For weak antigens, the membrane energy is dominated by the binding of adhesive complexes, and few antigen-TCR bonds can form because of molecular crowding. For very strong antigens, at intermediate time points ($\approx$ 5min) or with ligands of intermediate affinity, the membranes of the T-cells adopt an inverted bullseye pattern with TCR being in the p-SMAC and integrins being in the c-SMAC. Such patterns can be explained when the binding affinity of the integrin binding is larger than the binding energy of the weak ligands. In this case the integrins populate the center of the c-SMAC and make the deformation of the membrane ”affordable” energetically (see Figure 4). While this “bell-shaped” activation as a function of antigen affinity has been documented experimentally Valitutti1995a ; Kalergis2001 , it has not been consistently observed and there is a sufficient number of exceptions Holler2000 ; Lever2016 to diminish its significance, especially at the time when clinicians are engineering potent chimeric antigen receptors with unphysiologically–large affinity for their ligands Posey2016 ; Schmitt2017 . Third, the state of the membranes of both the T cell and antigen-presenting cell is expected to affect the synapse and thus T-cell activation, as was observed in experiments where the cytoskeleton has been depolymerized Valitutti1995b .

On the other hand, the model also fails to explain more recent observations, such as the ability of T cells to respond to a single agonist ligand while ignoring the presence of many non-agonist ligands. More recently, its interest has been rekindled by new observations emphasizing the importance of mechanics for ligand discrimination.

Using a biophysical setup to measure forces and lifetimes of individual pMHC-TCR pairs under tension, Zhu and colleagues Zhu2013 ; Liu2014 discovered two types of bonds depending on the peptide sequence: slip bonds and catch bonds. The lifetime of slip bonds decreases monotonously under tension: in other words, they break more easily when pulled. By contrast, the lifetime of catch bonds increases under tension. Molecular dynamics simulations Sibener2018 ; Wu2019 suggest that the bonds tension helps expose new residues in pMHC, inducing a better contact with the TCR. At high forces, catch bonds do ultimately break as well, so there exists an optimal force of around 10 pN for which the catch bond lifetime is longest.

This distinction between slip and catch bonds has been suggested to have functional significance in terms of T-cell activation. For instance, it was shown Sibener2018 that TCR-pMHC pairs with similar binding constants elicited different activation modes in T cells, with the stimulatory pair forming a catch-bond and the non-stimulatory a slip bond.

How these catch bonds relate to signalling remains to be explained in detail, but this new effect may be critical for predicting binding pairs, and immunogenic neo-antigens of crucial relevance to immunotherapies (see Sec. II.4). However, the ability of pMHC to activate the TCR signalling pathway might still be set by the lifetime of the pMHC-TCR complex, but with a “catch”: this lifetime would have to be assessed under the tension exerted by the membrane dynamics, e.g. 10 pN. This would also explain the dependence of activation upon membrane properties such as stiffness. The Ginzburg-Laudau model of Eq. 15 will need to be revisited to account for the existence of catch bonds, e.g. by letting the $k_{i}^{\rm off}$ of pMHC-TCR pairs depend on the tension $\lambda_{i}(z-z_{i}^{(0)})$, possibly affecting how ligands get sorted by mechanical forces in the immunological synapse.

Additional measurements in recent years (in particular, using super–resolution microscopy Cai2017 ) emphasize the active role that cytoskeletal rearrangements may play in concentrating TCR in the synapse Dustin2010 and the complex interplay between membrane ruffling, receptor sorting, and mechanical tensions. A more complete quantitative model would help identify key limiting steps deciding how strongly T cells get activated in different contexts — different antigen-presenting cells, or different co-stimulatory contexts.

III.2 B cells

As mentioned before, the affinity of B-cells with their cognate antigen undergoes a Darwinian selection process that improves their affinities from $K_{D}=10\,\mu$M down to 100pM. Simple counting of the number of occupied receptors in equilibrium could be sufficient to enforce antigen discrimination in B-cells. As for TCR binding to pMHC, the association rate of BCR to antigens is constant for all antigens. But unlike TCR-pMHC binding, which is subjected to an activation barrier, this rate is essentially diffusion limited, $k_{\rm on}>10^{6}M^{-1}s^{-1}$. Strong and weak antigens thus only differ by their binding lifetimes $\tau_{\rm off}$, which varies between $0.1s$ and $10^{4}$ s, suggesting a possible kinetic proofreading scheme as for TCR Tsourkas2012 . However, since the number of potential phosphorylation sites in the BCR complex is small compared to TCR, such a mechanism may not be as important.

The previously mentioned monomeric version of the BCR, IgG, is itself actually a dimer (so that IgM is a pentamer of dimers), with two binding sites. By analogy with other dimeric receptors on the surface of cells (e.g. Epidermal Growth Factor Receptors or Insulin Receptors Blinov2006 ; Lemmon2010 ), one could assume that this dimerization would help concentrate kinases (the enzymes responsible for phosphorylation) around the receptors, causing them to form clusters. This “cross-linking model” of BCR activation is appealing because it is consistent with the formation of microclusters of BCR on the surface of B cells in response to antigens, as observed by super-resolution microscopy.

An alternative model proposed by Reth and colleagues Yang2010 assumes that receptors cluster even in the absence of antigens. When in clusters, BCR inhibit each other and do not signal, but they can assemble and disassemble dynamically. Monomers coming out of these clusters are more prone to signaling, and antigen-binding stabilizes this monomeric state supposedly because of steric hindrances preventing the return into the clustered/inhibited state. This “dissocation-activation” model explains how B cells limit spurious activation of its $\sim 120,000$ surface BCR through their clustering, with signaling made only possible by isolated, antigen-bound monomers. The measurements of Reth and colleagues pose a theoretical challenge in terms of understanding the role dynamic clustering-release-re-clustering of BCR during B cell activation. A full theoretical model of such process and its impact on antigen discrimination remains to be proposed.

One observation that must be taken into account when considering B cell activation is the role of membrane spreading and contraction. Upon initiation of the signaling cascade, B cells rapidly ($<$ 10 min) reorganize their cytoskeleton and membrane, first by spreading to capture as many antigens on the presenting cells as possible, second by contracting to “concentrate” the active receptors. The quantitative model proposed in Fleire2006 accounts for binding and unbinding events driving cell spreading and contraction. The model predicts that such a dynamic process can increase dramatically the number of agonist ligands that get captured, compared to a static interface, enhancing the difference between weak and strong ligands, while weaker affinity antigens fail to concentrate. A critical aspect of the model is that, if B cells fail to trigger a sufficient number of BCR by 1 min (a hard cut-off), they terminate the process and shut down their signaling response. Alternatively, the activation response switches to a contraction phase with the surface area decaying according to the phenomenological law $A_{max}t^{-0.35}$ with $A_{max}$ being the largest area the B cell spreads to, $t$ is the time, and 0.35 is an experimentally-determined exponent. This model is phenomenological as it does not model explicitly the biochemical mechanism driving the spreading and contraction, and makes ad hoc assumptions about their behaviour. Yet it illustrates quantitatively how such membrane dynamics can help discriminate ligands.

To conclude, although the issue of antigen discrimination in B cells may not be as stringent as for T cells, it poses interesting quantitative issues for a different parameter range of binding constants ($K_{D}=100$ pM - 10 $\mu$M), over longer timescales ($>$10min), and using different cellular mechanisms: receptor clustering and membrane dynamics. It is worth recalling that this initial recognition process connects to additional processes over longer timescales, such as validation by helper T-cells, and dynamics within germinal centers. We will come back to these questions in Sec. VIII.

III.3 Coarse-graining of molecular details and model reduction

As we go deeper into the molecular details of immune recognition, the number of molecular species, reactions, intermediates, and therefore parameters explodes. This poses a challenge for a number of reasons. First, the effective behaviour of the system as a whole may still be relatively simple, suggesting that simpler phenomenological models may describe these processes equally well. While the variables of such models may be hard to relate directly to molecular entities, they are easier to interpret and allow for better analytical progress and predictions. Second, even assuming that the full complexity of all interactions is needed, large numbers of parameters are likely to lead to overfitting problems, meaning that many parameters or combinations of parameters are underdetermined, undermining the accuracy of predictions. And even when they can be determined, it is not always clear which ones need to be fine-tuned to ensure proper function, or what are the broad design principles presiding over their choices. In Sec. XII.5.2, we briefly review the principles of model selection – the classical approach of reducing model complexity from statistics. However, that approach requires to have first defined a hierarchy of models to test, from least to most complex. Besides, model selection relies on goodness of fit as a criterion to evaluate models, while in many cases we may be more interested in capturing the principal features of a biological function, rather than fitting all the data.

In another approach to reducing complexity, François and Hakim hakim-2004 developed an method for generating simple molecular networks in silico that realize a desired biological function, simply based on a genetic algorithm that selects the “best” solutions. Applied to the problem of absolute ligand discrimination reviewed in Sec. III.1, this method infers a class of network motifs, called “adaptive sorting”, that recapitulates known features of T cell activation Lalanne2013 . In particular, it predicts the emergence of kinetic proofreading and biochemical adaptation. However, the chemical species and reactions of the networks produced by that method may generally not be directly related to the known actors of the phenomenon under study.

To reduce model complexity while keeping close to the details of actual biochemical reactions, one can instead start from a complex biochemical network described by many ordinary differential equations, and “prune” its parameters by setting them, individually or by their combinations (ratio or products), to 0 or infinity Proulx-Giraldeau2017 . Applying this strategy to a complex model of T cell recognition Lipniacki2008 , which contains close to 100 parameters, shows that its behaviour can be boiled down to just three coupled differential equations highlighting its main features of adaptation and discrimination, and reveals the broad design principles that implement these features.

Beyond this particular example, this kind of approach has great potential for helping to make sense of complex biological systems with many entities and interactions, as is often the case in the immune system.

IV Cell-to-cell communication through cytokines

Cytokine communication is critical to synchronizing the activation of various immune cells and to bridging multiple spatial-temporal scales in immune responses, from local or individual cell activation to global, systemic responses. Cytokines are small glycoproteins that get produced and secreted by all cells, immune or not, with varied dynamics, amplitude and frequency. These cytokines then diffuse and bind to receptors present on the surface of adjacent cells to elicit a signaling response that trigger a gene regulatory response. In short, cells use cytokines to communicate between themselves. In this section, we discuss quantitative models that have been introduced to model how individual leukocytes respond to cytokines (*e.g. *using the JAK-STAT pathway) over multiple time- and length-scales.

IV.1 Cytokine signaling and the JAK-STAT pathway

Cytokine receptors are heterodimers, or more rarely heterotrimers, so that they carry at least two intracellular signaling domains. The JAK-STAT pathway is the dominant pathway engaged by cytokine signaling. Upon binding by its cognate cytokine, the receptor undergoes a conformational change that presents phosphoryl attachment sites on the JAK (Janus Kinase) receptor to face the kinase domain that induces phosphorylation, which in turn leads to activation. Then a protein called STAT (Signal Transducer and Activator of Transcription) interacts with the intracellular domains of the receptors, gets phosphorylated by the activated JAK, and dimerizes upon release from the receptor. Dimers of STAT then translocate into the cell nucleus, bind to the chromatin in specific sites and elicit a transcriptional response. Such signal transduction is one of the simplest biological pathway connecting the extracellular environment and its messenger cytokines to a transcriptional response (Figure 5). Its complexity is encoded in the large number of pairs of cytokines ($>40$) and their receptors that cells can express, as well as the 7 forms of STAT and $7^{2}=49$ homo- or hetero-dimers that they can form upon activation Altan-Bonnet2019 . Additionally, the dynamics of cytokine signaling enriches the biology of these pathways, with the existence of positive and negative feedback regulations: cytokine signaling inducing the expression of additional cytokine receptors, cytokine degradation, cytokine receptor endocytosis, expression of negative regulators such as Suppressor of Cytokine Signaling (SOCS). Here we present simple derivations that can help understand the quantitative regulation of cytokine communication within the immune system, which underlies the coordination and orchestration of the immune response.

IV.1.1 Cytokine binding and signaling at equilibrium

In this section, we discuss the contribution of the field of quantitative immunology to understanding the JAK-STAT pathway as a signal transduction cascade in well-mixed conditions. Spatial heterogeneities in cell-to-cell communication via cytokines will be tackled in section IV.2.

Within an individual cell, the biochemical reactions of the JAK-STAT pathway can be treated as well-mixed: the typical concentrations of molecules are fairly high, with $\sim 1000$ molecules per cell of diameter 10$\mu$m, which translates into concentrations $\sim 100$nM. We also assume that most reactions are essentially diffusion-limited, as discussed before: $k_{\rm on}>10^{7}M^{-1}s^{-1}$. The signal transduction cascade for JAK-STAT cascade can generically be modelled in a step-wise manner, as sketched in Figure 5.

The kinetics of cytokine binding to its receptor activating JAK into ${\rm JAK^{*}}$ can be modeled as a simple bimolecular process:

[TABLE]

where ${\rm\#Receptor^{*}}$ is the number of cytokine-engaged receptors, $[{\rm Cytokine}]$ is the cytokine concentration. Because binding is very strong (typically $K_{D}=k_{\rm off}/k_{\rm on}\approx 10pM$ with $k_{\rm off}^{-1}\approx 3600$ s and $k_{\rm on}\approx 3\cdot 10^{7}M^{-1}s^{-1}$), most cytokine molecules bind at a diffusion-limited speed with very strong binding. We assume a large extracellular volume, so that cells do not deplete cytokines as they bind them. Receptors are usually pre-loaded with JAK, waiting for the conformational change associated with cytokine binding, to induce phosphorylation. From Eq. 18 the number of activated JAK (JAK∗) at steady-state is:

[TABLE]

where ${\rm JAK}_{\rm total}$ is the total JAK receptor concentration. The kinetics of STAT phosphorylation into pSTAT is described using a classical Michaelis-Menten equation:

[TABLE]

where ${\rm[JAK^{*}]}={\rm JAK^{*}}/(V\mathcal{N}_{A})$ is the concentration of activated JAK within the cytoplasmic volume $V$, ${\rm pSTAT}$ is the phosphorylated form of STAT, $k_{\rm cat}$ and $k_{\rm dephos}$ are phosphorylation and dephosphorylation rates, and $K_{m}$ is the dissociation constant of STAT to the receptor, which controls the specificity between different receptors and different variants of STAT. Note that, for most signaling networks, the identity of phosphatases that take care of dephosphorylating receptors and kinases remain often undetermined (because of overlapping and pleiotropic activities): their activity is modeled phenomenologically as a single rate reaction.

This system reaches steady-state within minutes Vogel2016 ; Altan-Bonnet2019 such that one can solve $d{\rm[pSTAT]}/dt$ in Eq. 20 for pSTAT (with the conservation of matter condition [pSTAT]+[STAT]=$S_{\rm total}$) as a quadratic equation. Combined with Eq. 19, the result gives a direct expression of the response STAT as a function of the input cytokine concentration. This model serves as the basic building block to tackle the dynamic complexity of cytokine responses.

IV.1.2 Tunability of cytokine responses.

In the previous section we assumed that cytokines, which are formed of several sub-units, were in a pre-assembled configuration. When that is not the case, cytokine binding to the cytokine receptor can be modeled as a two-step scheme.

[TABLE]

The cytokine ligand L binds to the cytokine receptor composed of two parts R1 and R2. A simple equilibrium model can be solved, with the dissociation constants $K_{i}=k_{i}^{\rm off}/k_{i}^{\rm on}$:

[TABLE]

that yields the total concentration of fully-engaged and signaling receptors:

[TABLE]

This two-step model was shown to be valid for one of the most studied cytokines, IL-2 Feinerman2010 . IL-2 binds weakly to the abundant $\alpha$ chain of the IL-2 receptor (${\rm R_{1}=IL2R}\alpha$), before being “locked” into a stable and signaling complex by association with the two other chains (${\rm R_{2}=IL2R}\beta$-$\gamma_{C}$) to form the full IL-2$-$IL2R complex. In that case, $K_{1}[L]\ll 1$ such that equation 22 simplifies to:

[TABLE]

This equation reveals key insights for IL-2 and other cytokines Cotari2013 : the amplitude of cytokine signaling is proportional the number of $\beta-\gamma$ part of the IL2 receptor (the limiting part), and reaches half of its maximum at $[L]_{50}=K_{1}K_{2}/[R_{1}]_{\rm total}$, inversely proportional to the $\alpha$ chain of the IL-2 receptor (the non limiting-part).

Such tunability of both the amplitude and sensitivity of the dose response of cytokines based on the number of cytokine receptor chains can be critical to achieve immune plasticity, as demonstrated in the context of competition for limited amounts of cytokines. For example, the tug-of-war for IL-2 between effector T-cells and regulatory T-cells is mediated by the exact levels of IL2R$\alpha$ on the surface of these cells. Stimulated effector T-cells initiate an immune response, whereas stimulated regulatory T cells (Treg) downregulate the response to suppress auto-immune responses (effector response to self-antigens). In general, Treg cells are thought to bypass negative thymic selection despite their strong recognition of self-antigens. These cells then act as pre-activated sentinels that respond to inflammatory cues (such as the upregulation of self-antigens) while not contributing to inflammation itself through cytokine secretion: Tregs constitutively express transcription factor FoxP3, a general downregulator of cytokine production.

Whenever an effector T-cell initiates a spurious response to self-antigens, Tregs can extinguish it by consuming cytokines (*e.g. *IL-2) and by downregulating their inflammatory impact. In this tug of war, whichever cell type, effector or regulatory, expresses more receptors boosts its ability to capture the cytokine and deprives the other cells of this key cytokine for anti-apoptosis and proliferation in the T-cell compartment. Quantitative models of such competition for IL-2 based on differential expression of IL-2R$\alpha$ receptors have illuminated how self vs. non-self discrimination can emerge from such IL-2 tug–of–war Busse2010 ; Feinerman2010 ; Hofer2012 .

Some cytokine receptors have an extracellular domain that gets secreted extracellularly, in a soluble form. These can “pre-bind” the cytokine in the extracellular milieu, and deliver it to the complementary chain bound to the cell membrane: such molecular event can have positive or negative effects on cytokine signaling depending on the context. If the pre–formed soluble cytokine/receptor complex binds to the cell to form an incomplete receptor which lacks the intracellular signaling domain of the soluble cytokine receptor, the JAK misses its trans-phosphorylation partner and the cytokine/cytokine receptor fails to signal— “no clap from one hand.” Such soluble complexes act as decoys that antagonize cytokine response (as with viral analogues of IFN receptors, or IL-1 receptors), limiting inflammation to the benefit of viruses Levine2004 .

Alternatively, secretion of soluble cytokine receptors can trans-activate cytokine signaling and result in a boost in cytokine response (as with IL-1, IL-2 or IL-6). For example, in the case of IL-6, the soluble portion of the IL-6 receptor (sIL-6R) gets secreted to stabilize IL-6 in the extracellular medium (most cytokines are very small proteins with short half-lives) and to accelerate the assembly of a complete IL-6/IL-6R/gp30 signaling complex by binding to gp130 dimers on the surface of cells: such activation of the IL-6 signaling response thus does not require the presence of the IL-6 receptor in the membrane of the receiving cells, but only the ubiquitous presence of gp130: this is called cytokine trans-activation Rose-John2004 .

Hence, depending on the exact molecular details, secretion of cytokine receptors in the extracellular environment can trigger or antagonize the cytokine signaling response. The role of quantitative immunology in that context is then to tease out the physico-chemical parameters of such regulation to better understand in which regime inflammatory cues are regulated.

IV.1.3 Regulation by cytokine consumption

Immune responses need tight regulation to avoid spurious auto-immune activation. This is particularly important for cytokine regulation as overabundant cytokines can be extremely deleterious to the organism as a whole. For example, high concentrations of cytokines induce a “capillary leak syndrome,” whereby tissues lose their barrier against blood serum: this results in septic shocks or viral hemorrhagic fevers. One simple mode of regulation of cytokine signaling by the JAK-STAT pathway is simply to limit the availability of cytokines in the extracellular medium. Most cytokines are very small proteins (with molecular weights in the 10-20kDa range) such that their half-life in the bodily fluid is short. For example, the half-life of IL-2 injected intravenously for the treatment of renal cell carcinoma is at most 3 h in the serum. More significantly, even in the extracellular environment of lymphoid organs, cytokines have a very short half-life due to rapid binding and consumption.

Given the strong binding of cytokine to their receptors (typically $K_{D}<100$ pM or $k_{\rm off}<(3000s)^{-1}$), immune cells can rely on binding to switch off the cytokine signaling response through buffering. This resetting process relies on cytokine consumption by endocytosis of the cytokine with its receptor, trafficking them towards lysosomes, release of the cytokine because of low pH within the lysosomal compartment, and degradation. Note that the signal transduction may persist after receptor endocytosis as long as the cytokine-receptor pair remains engaged Cendrowski2016 . What happens to the cytokine receptor in this process varies based on the cytokine identity. For example, the endocytosed IL-2 receptor dissociates and has its chains sorted towards different compartments: the IL2R$\gamma_{C}$ chain goes to the lysosome and gets degraded, while the IL2R$\alpha$ chain goes to early endosomes and gets recycled. This process of endocytosis, sorting, recycling and degradation can be very complex at the molecular level. Yet, a coarse-grained model of this process as a single-biochemical step can be sufficient when modeling IL-2 availability Tkach2014 ; Voisinne2015 : the typical rate for this step has been measured to be of the order of (900 s)-1.

One functional consequence is that cytokines get rapidly consumed while triggering signaling. At the more global level, one can integrate a dynamic equation for production/consumption that accounts for the rise and decay of inflammatory signals in the immune system:

[TABLE]

In the most simple case, $\kappa_{\rm prod}=N_{\rm prod}k_{\rm prod}/(V{\mathcal{N}}_{A})$ is fixed by the number $N_{\rm prod}$ of activated cells at a given time, the rate of secretion $k_{\rm prod}$ per individual cells, and the extracellular volume $V$. The rate $\kappa_{\rm consum}$ is determined by the rate of binding of the cytokine to its receptor on the surface of the cells:

$\kappa_{\rm consum}=k_{\rm on}N_{\rm consum}N_{\rm R}[{\rm Cytokine}]/{V{\mathcal{N}}_{A}}$,

where $N_{\rm consum}$ is the number of consuming cells and $N_{\rm R}$ the number of receptors per cell. The time dependency of $N_{\rm prod}$, $N_{\rm consum}$, $\kappa_{\rm prod}$, and $N_{\rm R}$ can be arbitrarily complicated and must be parametrized for each immunological setting under consideration. For example, in the case of IL-2 in the early events of an immune response within a lymph node, $V\approx 50\mu\ell$ (free volume), $k_{\rm prod}=10{\rm s}^{-1}$, $k_{\rm on}=3\cdot 10^{7}M^{-1}s^{-1}$, $N_{\rm consum}=10,000$ Treg cells, each endowed with typically $N_{\rm R}=3,000$ receptors. For IL-2 to accumulate and reach a significant concentration (typically, $10^{-11}M$ to trigger STAT5 phosphorylation), one must have $\kappa_{\rm prod}>\kappa_{\rm consum}$, i.e. $N_{\rm prod}>1,000$. This estimate illustrates that there exist thresholds of activation for immune responses, whereby the systems needs a critical mass of activated, cytokine-secreting cells to overcome consumption and drive activation and differentiation Polonsky2018 . The estimate above does not take into account some of the intricacies of IL-2 regulation: positive feedback in IL-2 secretion, recycling of IL2R$\alpha$ chains, upregulation of IL-2R in Treg cells have been documented. To account for these, Eq. 24 must be solved numerically in more complex settings Tkach2014 ; Voisinne2015 .

IV.1.4 Other regulations

Cytokine consumption is one major mechanism to limit the duration of availability of cytokines in the extracellular medium, but there exist additional, cell-intrinsic mechanisms that also limit or expand the duration of JAK-STAT signaling in cells, as illustrated by the following examples.

Certain cytokine pathways rely on intracellularly stored pools of receptors that get recruited upon receptor engagement and JAK-STAT signaling. This recruitment further fuels JAK-STAT signaling by avoiding saturation of the initial receptors on the membrane surface. This process was analyzed quantitatively in Ref. Becker2010 to reveal how a certain type of leukocyte called erythroid progenitors can sense a cytokine called Epo over a large range of concentrations ($>1000$-fold range), where classical ligand-receptor binding would predict fast saturation with increasing cytokine concentration.

Another positive feedback in JAK-STAT signaling exists in the IL-2 signaling pathway: The $\alpha$ chain of the IL-2 receptor gets upregulated upon IL-2 signaling, lowering the cytokine concentration at which activation is half-maximum (by virtue of Eq. 23), hence driving further IL-2 signaling Busse2010 ; Feinerman2010 ; Cotari2013 . Such positive feedback was shown to be critical to commiting cells to long-term JAK activation and proliferation Voisinne2015 .

Alternatively, JAK-STAT pathways are also endowed with negative regulators, such as SOCS and Cytokine Inducible SH2-containing (CISH) proteins, which compete with STAT for JAK binding, preventing activation Raue2011 . This can be quantitatively modeled with SOCS or CISH binding competitively with a stronger $K_{m}$ than STAT onto the receptor (Eq. 20). Varying the exact value of $K_{m}$ allows SOCS and CISH to negatively regulate JAK activity in two separate regimes of concentrations of cytokines Bachmann2011 . Such dual feedback regulation of signaling was also shown to extend the range of cytokines that trigger STAT5 phosphorylation while avoiding saturation, and ultimately control cell survival. The modeling of such complex regulation of the JAK-STAT pathway involves adding biochemical steps that can be integrated numerically. The challenge in this context is to acquire a large number of system-specific biochemical parameters. This task has become more amenable as the field progresses and quantitative approaches have been delivering more and more estimates: protein abundances, biophysical parameters, enzymatic rates Cotari2013 ; Karr2015 ; Shi2016 ; Mitchell2018 .

To conclude, in this section we have reviewed the basic equations governing the regulation of the cytokine-activated JAK-STAT pathway in the immune system. While JAK-STAT signaling can be solved analytically in the adiabatic limit, the long-term dynamics of cytokine accumulation and consumption are complicated by the multitude of feedback regulations triggered in leukocytes. Such rich dynamical complexity will need to be tackled quantitatively and systematically to deliver a more comprehensive understanding of cytokine communication in the immune system Altan-Bonnet2019 .

IV.2 Communication across space and time, and cytokine niches

Cytokines allow immune cells to communicate and to modulate their response collectively. In section V.1, we will discuss how lymphocyte differentiation can be decided by toggle switches encoded in gene regulatory networks. Most of the positive feedbacks in these toggle switches are in fact associated with a response to cytokines, and are thus intrinsically collective. Tackling the nonlinearities and spatio-temporal aspects of cytokine communication at a more quantitative level is critical to expanding our understanding how antigen recognition leads to such a collective response. We now focus on the dynamics of cytokine secretion and capture in space and time, over the relatively short timescales (minute to hours) over which cytokines convey information between cells.

IV.2.1 Cytokine dynamics as diffusion-degradation

The propagation of cytokines within tissues is governed by reaction-diffusion equations which generalize Eq. 24:

[TABLE]

where $c(\vec{r},t)$ is the spatio-temporal profile of the concentration of cytokine, $D$ is the diffusion coefficient for the cytokine in the extracellular medium, $\kappa_{\rm prod}$ and $\kappa_{\rm cons}$ are volumic rates of cytokine production and consumption respectively. Cytokines are produced by discrete cells, and the production rate $\kappa_{\rm prod}(\vec{r},t)$ concentrates on the cell surface. The consumption rate can be estimated by taking into account that cytokines bind tightly to cytokine receptors on the cell surface and get consumed by endocytosis of the engaged receptors. The limiting step for this cytokine consumption is then the binding to receptors, so that

[TABLE]

where $n_{\rm consum}$ is the density of consuming cells. This expression is valid at low concentrations, when cytokine receptors are not saturated and the kinetics of binding remain linear with the concentration of cytokines. It also assumes that consuming cells are uniformly distributed, in a mean-field manner, rather than modelling the precise location of their surfaces. In more general cases, $\kappa_{\rm prod}$ and $\kappa_{\rm consum}$ may have additional, nonlinear dependencies on $c(\vec{r},t)$, for instance because of spatial heterogeneity or receptor saturation. Little can be done analytically in these cases and one must resort to numerical solutions Thurley2015 . In the linear regime, Eq. 25 becomes:

[TABLE]

with a characteristic length $\xi=\sqrt{{D}/{k_{c}}}$ over which cytokines diffuse before being consumed. In the following section, we focus on the simple homogenous and linear case to derive insights about the quantitative aspects of cytokine communication.

IV.2.2 Screening by cytokine-consuming cells

Diffusion and consumption of cytokines regulates how the cytokine gradient spreads i.e. the lengthscale for cell communication. It was demonstrated that the short timescale to reach stationary concentration profiles for soluble proteins are explained by the first arrival time of the cytokine ligands, rather than the characteristic diffusion timescale for an individual molecule Berezhkovskii2010 ; Kolomeisky2011 . This result implies that the characteristic timescale to reach steady state at a distance $r$ from the secreted cell scales linearly with $r$ rather than with $r^{2}$. The analytical expression for this local relaxation time $\tau(r)$ is Berezhkovskii2010 :

[TABLE]

Two relaxation processes are at play in this system: the diffusion timescale ($\tau_{d}(r)=r^{2}/2D$) and the reaction timescale ($k_{c}^{-1}$). The relaxation timescale is dominated by the reaction time at short distances ($r\ll\xi$) where cytokine consumption dominates. However, at large distances ($r\gg\xi$), both diffusion and consumption are important and the relaxation timescale is the geometric mean of the two, $\tau(r)\sim({{\tau_{d}(r)k_{c}^{-1}}})^{1/2}\propto r$. Numerical simulations in the context of immune responses Oyler-Yaniv2017 validated this theoretical argument and estimated that it takes only a few minutes for the cytokine gradient to reach its steady state.

The timescale at which immunological regulation by cytokine exchange operates is typically in the range of hours for the activation of gene regulatory elements, to days for the decision to proliferate or to die. The typical timescale at which cells move is also very slow. In the first phase of an adaptive immune responses, antigen-responding T-cells are essentially stuck on their antigen presenting cells. One can thus assume steady state for the distribution profile of the cytokines, $c(\vec{r},t)=c(\vec{r})$, and solve for $dc/dt=0$ in Eq. 27. Using spherical symmetry $c(\vec{r},t)=c(r)$ for the Laplace operator, the steady state profile for the cytokine distribution is

[TABLE]

and the boundary value $c(R_{\rm cell})$ can be determined by estimating the flux of molecules secreted and consumed at the cell surface:

[TABLE]

Typically, cells secrete cytokines at a constant rate between $k_{\rm prod}=$10 to 1000 molecules per second (e.g. IL-2 and IFN-$\gamma$ cytokines) and consumption by the secreting cell is given by $k_{\rm consum}=k_{\rm on}N_{\rm R}c(R_{\rm cell})$. Solving for $c(R_{\rm cell})$ yields:

[TABLE]

The characteristic scale

[TABLE]

is analogous to the screening length of electrostatic interactions and is inversely proportional to the square root of the density of consuming cells. By this analogy, consuming cells “screen” the diffusion of cytokines from the secreting cells and end up determining the extent of cell-to-cell communication in a dense tissue. In other words, the simple diffusion-consumption of cytokines can account for the heterogeneous accessibility of cytokines in dense tissues and the formation of cytokine “niches,” defined as the portion of space in where cytokines can be sensed.

$\xi$ is controlled by the density of cytokine-consuming cells, and the number $N_{\rm R}$ of receptors on the surface of consuming cells. In close-packed tissues the density of consuming cells is typically 1 cell every 10$\mu$m, implying $n_{\rm consum}\sim 10^{-3}\mu{\rm m}^{-3}$; for IL-2 $D$ is approximately 100$\mu m^{2}/s$ (even in close-packed tissues, there is sufficient free space in the extracellular medium for cytokines to diffuse freely as in solution), and each cell consumes typically 10 cytokine molecule per second Oyler-Yaniv2017 , implying $k_{\rm on}N_{R}=10^{11}M^{-1}s^{-1}$, so that $\xi\sim 25\mu{\rm m}=2.5$ cell diameters. In more physiological conditions, where cells are more dilute, occupying only 10% of the space and cells express lower level of receptors, $n_{\rm consum}\sim 10^{-2}\mu{\rm m}^{-3}$, this length scale rises to $\xi\sim 250\mu{\rm m}\sim 25$ cell diameters: there is negligible screening and all the cells within a lymph node or a spleen have access to the cytokine. Note that this 10-fold increase in the screening length $\xi$ implies an increase by 1000-fold in the volume of the cytokine niche and a 1000-fold increase in the number of cells that potentially respond to the cytokine. Such tuneability can be critical to decide which cells (helper T cells or regulatory T cells) win the tug of war for IL-2 and ultimately whether the immune response expands or gets extinguished Busse2010 ; Feinerman2010 ; Hofer2012 .

IV.2.3 Probability of autocrine capture

Cells that secreted cytokines often express the associated receptors, and can in principle communicate with themselves. To quantify this effect, we must estimate the non-zero probability $P_{\rm auto}$ that cells capture their own cytokines rather than let them diffuse and interact with neighbouring cells. This phenomenon, called autocrine signaling, is particularly significant when cytokine-secreting cells are surrounded by cytokine-consuming cells and compete for their consumption.

At steady state, there is a balance between cytokine production $k_{\rm prod}$, cytokine consumption by the secreting cell $k_{\rm consum}$, and diffusion away from the cell followed by consumption by other cells, as given by Eq. 30. The probability of absorption by the secreting cell is given by:

[TABLE]

We can define the characteristic number of receptors at which this probability reached 50%, $N_{R}^{*}=4\pi DR_{\rm cell}(1+R_{\rm cell}/\xi)/k_{\rm on}\approx R_{\rm cell}/\sigma_{\rm R}$, where $\sigma_{\rm R}=R_{\rm receptor}e^{-\Delta G_{\rm assoc}}$ is the effective cross-section of the receptor-cytokine interaction (Eq. 2 and 6), and where we assumed $\xi\gg R_{\rm cell}$.

Leukocytes are typically $R_{\rm cell}=5$ to 10 $\mu$m in radius, and $\sigma\sim 0.5$nm, hence $N^{*}_{R}\approx 10^{4}$. If cells have much fewer than $10^{4}$ receptors (a common situation), $P_{\rm auto}\ll 1$ and most of the secreted cytokines diffuse away without being captured in an autocrine manner. Alternatively, if cells have more than $10^{4}$ receptors, $P_{auto}\sim 1$ and secreting cells capture a large fraction of their own cytokines. The functional relevance of autocrine capture remains debated in the field: why would a cell need to respond to the cytokine it secretes. There are also molecular mechanisms limiting its impact. For instance, T cells secreting IL-2 upon antigen recognition express low levels of IL-2 receptors as long as the antigen response lasts to limit autocrine capture Tkach2014 ). In the tug of war between cell types for cytokine consumption, competition between autocrine and paracrine signaling is crucial, and thus the amount of autocrine capture is significant Hofer2012 . Discriminating between autocrine and paracrine signaling is often difficult and it has been suggested that dilution experiments in vitro can be used to distinguish the two modes Marcou2018c .

IV.2.4 Size of cytokine niches

Eqs. 29-31 can also be used to estimate the functional size of a cytokine niche established by a single secreting cell. In the limit of no autocrine capture and long dispersal $\xi\gg R_{\rm cell}$, we have:

[TABLE]

The radius of the niche $r_{\rm niche}$ can be defined as the $r$ at which $c(r)$ falls below the characteristic concentrations $K_{D}$ needed to trigger cytokine signaling, and is the solution to the implicit equation: $r_{\rm niche}=\xi\ln(k_{\rm prod}/4\pi rK_{D})$. The logarithmic dependence on the rate of cytokine secretion $k_{\rm prod}$ is significant because that rate may vary greatly depending on the cytokine or depending on the state of differentiation of cells. For example, CD8+ T cells can secrete up to 1000 IFN-$\gamma$ per second, CD4+ T cells secrete 10 IL-2 per second upon activation from a naive state, and 1000 IL-2 per second as they mature Tkach2014 or when get activated from a memory state Huang2013 .

In Figure 6, we present the numerical solution for $r_{\rm niche}$ for different secretion rates for cytokines with different screening lengths. Note how variable the volume of the cytokine niche can be, based on realistic parameters of immune responses.

To conclude, in this section we discussed how simple equations for the diffusion and consumption of cytokines in dense immunological tissues can quantitatively account for the tunability of cell-to-cell communication in the immune system. It is particularly striking that such active modulation of cell communication plays a functional role in modulating immune responses. In the first moments of activation, there are very few cytokine-consuming cells, and the field of cytokines extends to the entirety of the lymphoid organ ($\xi>100\mu$m), reaching blood vessels and making communications global. Within 24 hours, many cells migrate closer to cytokine-secreting cells. As they do so, they increase their cytokine consuming capabilities by upregulating their receptors, thereby “screening” the diffusion of cytokines ($\xi\sim 10\mu$m) and limiting the extent of cell communication to nearest neighbours, within a tight niche around secreting cells.

V Cell fate

There is a large variety of types and subtypes of immune cells, making the field of immunology sometimes difficult for the non-specialist. These fates are acquired either at the very beginning during haematopoiesis, the process by which all types of blood cells are generated by differentiation from stem cells, or following an immune response or challenge, during which cells further specialize to better fight infections. We first describe the physical mechanisms of gene regulation that lead to stable distinct cell fates. We then review modeling strategies for understanding and inferring the hematopoiesis differentiation tree. Finally, we discuss the problem of differentiation during an immune response.

V.1 Gene regulation and cell differentiation

Many cell decision making processes in the immune system are carried out at the level of gene regulation, whereby expression of key transcription factors decide cell fate and subsequent immune effector function. For example, T-cells sense their inflammatory environment as defined by the combinations of cytokines in the extracellular environment, and drive their signaling response towards expressing key transcription factors that in turn decide which cytokines get produced and secreted.

To consider a concrete example of gene regulation, we use the classical example of the differentiation of CD4+ helper cells, a particular type of T cells whose function is to regulate the action of other immune cells through the production of cytokines — the messenger molecules used by immune cells to communicate with each other. In their naïve state, T cells are essentially a blank slate which can later differentiate into more specialized cell states to best match the pathogenic threat. Upon antigen stimulation, helper T cells can orchestrate three types of immune responses. They can elicit a ‘Th1’ response, further unleashing T- and NK cells’ cytotoxic response against infected cells. This is particularly relevant when eradicating an intracellular infection such as a viral infection. Alternatively, helper T cells can elicit a ‘Th2’ response to engage with B cells, drive affinity maturation and class switching in antibody production towards secretion of antibodies that can annihilate extracellular (e.g. bacterial) infections. Finally, T-cells can elicit a ‘Th17’ response that protects mucosal barriers from infection. These three characteristic types of T-helper cell differentiation respond to and amplify distinct cytokine environments. Note that similar decisions are made by other immune cells: for example macrophages differentiate into ‘M1’ (pro-inflammatory) and ‘M2’ (anti-inflammatory) types, which mirror the decisions and inflammatory outputs associated with the Th1 and Th2 helper T-cell types.

Immunologists are interested in dissecting the interplay and feedback between inflammatory environments and immune cell differentiation Zhu2010 . In that context, quantitative immunology brings the tools of nonlinear dynamics to explain how sharp and cross-inhibitory decisions can be made by activated T-cells during differentiation. Many immunological systems involve nonlinearities with self-reinforcing feedback loops, whereby the protein of interest acts as a positive feedback that drives further expression. Additionally, most gene regulation involve mutltimerization of transcription factors onto one gene locus, as in the case with phosphorylated STAT, a transcription factor induced by cytokine signaling. Multimerization implies a nonlinearity because the association rate of an n-mer scales with the nth power of the monomer’s concentration. These nonlinearities and their functional consequences for cell fate in the immune system can be studied using nonlinear stability analysis. Many transcription factors auto-amplify their production in a direct loop or, more generically, in a more convoluted cytokine-mediated manner. For example, the transcription factor Tbet drives the secretion of a cytokine called interferon $\gamma$ (IFN-$\gamma$) that signals through pSTAT1 to produce more Tbet in the context of Th1 T-cell differentiation. In section XII.2.1 of the Methods, we recall the basic formalism of nonlinear analysis, which explains the generic occurence of bimodal distributions of expression of transcription factors, cytokines and surface markers that immunologists encounter in their single-cell measurements (cytometry or single-cell transcriptomics). We discuss classical models of gene regulation and feedback regulation that account for the molecular programs enforcing sharp cell differentiation in the immune system.

V.1.1 Th1/Th2 differentiation

Early experimental evidence demonstrated that CD4+ helper T-cells commit to two distinct and incompatible states of differentiation, Th1 and Th2, upon antigen activation and response to inflammatory cues. A model describing this system is based on two inhibitory loops (see Sec. XII.2.1) whereby transcription factor $P_{1}$ inhibits the production of $P_{2}$ and vice versa (Figure 7A & B).

[TABLE]

where $R_{1}$ and $R_{2}$ are number of mRNA leading to the production of $P_{1}$ and $P_{2}$, $\gamma$ decay rates, $k_{i}$ and $v_{i}$ are production rates, $K_{i}$ are binding dissociation constants, and $n_{i}$ are Hill coefficients. At steady state,

[TABLE]

to which some background expression may be added to account for leaky gene regulation.

The nonlinearities, which are controlled by the Hill coefficients $n_{i}>1$, originate from the multimeric control of transcription by the transcription factors. A graphical inspection of the nullclines (Eq. 35) for this dynamical model (Figure 7.C) reveals that there exist three fixed points for this gene regulatory system, two stable ones, $(P_{1}^{high},P_{2}^{low})$ and $(P_{2}^{high},P_{1}^{low})$, and one unstable fixed point $(P_{1}^{high},P_{2}^{high})$. This simple dynamical system reveals how cross-inhibition of transcription factors yields a classical toggle-switch that yields two incompatible states. Note that such a model predicts hysteresis in T cell differentiation (Figure 7D) as has been observed experimentally Chaouat2004 ; Yates2004 ; Antebi2013 .

V.1.2 Other differentiation switches

A similar formalism has been applied to account for the differentiation of CD4+ T-cells in Th17 and regulatory T-cells (Treg). Th17 cells are pro-inflammatory and critical to eradicate infections in mucosal tissues; regulatory T cells are anti-inflammatory and limit overzealous immune responses, such as auto-immune disorders caused by a strong response against self antigens, or runaway immune responses such as septic shock. These two states of differentiation for CD4+ T cells can be accessed upon exposure to the cytokine TGF$\beta$. However, that cytokine has been shown to be both pro- or anti-inflammatory depending on the context. Tyson & coworkers Hong2011 proposed a model to account for such bimodality in TGF-$\beta$ action. In that model, the concentration of each protein species $P_{i}$ generically follows:

[TABLE]

where $\omega_{i}=\omega_{i}^{0}+\sum_{j}\omega_{j\rightarrow i}P_{j}$, $P_{i}^{\rm eq}$ is the characteristic concentration of $P_{i}$ in differentiated cells at steady state, $\gamma_{i}$ are degradation rates, $\omega^{0}_{i}$ is the constitutive level of expression of regulation of $X_{i}$, $\omega_{j\rightarrow i}$ is the (possibly negative) influence of $X_{j}$ on $X_{i}$, and $\sigma_{i}$ is a parameter modulating the steepness of the regulation. Note that proteins can regulate their own production, $\omega_{i\to i}\neq 0$.

There are three main species involved in Th17/Treg differentation: $P_{1}$=FoxP3, $P_{2}=$ROR$\gamma$t (the transcription factors controlling Treg and Th17 differentiation, respectively), and $P_{3}=$TGF$\beta$ (a cytokine that regulate T cell differentiation), which acts as an externally fixed stimulus. Their interactions are represented in Figure 8A. For instance ROR$\gamma$t+ is activated by TGF$\beta$ ($\omega_{3\to 2}>0$), as well as by itself ($\omega_{2\to 2}>0$), and but is repressed by ROR$\gamma$t ($\omega_{1\to 2}<0$).

Graphically, one can represent the nullclines ($dX_{i}/dt=0$) for both transcription factors as shown in Figure 8. There are three stable fixed points where the nullclines meet: ROR$\gamma$t+FoxP3- for the Th17 cells, FoxP3+ROR$\gamma$t- for the Treg cells and ROR$\gamma$t+FoxP3+ for a mixed state whose existence was later confirmed experimentally Hong2011 . This formalism is interesting because it handles coupling and feedback phenomenologically, through the effective regulation parameters $\omega_{i\to j}$, rather than through detailed mechanisms of how these regulations are implemented.

V.1.3 Experimental test of bistability in cell differentiation

The elegant models presented in the two previous subsections (Th1/Th2 and Th17/Treg differentiation) account qualitatively for the observation that many activation conditions lead T-cells to commit to clear, separate and stable cell fates. The crux of these models is the existence of toggle switches based on cross-inhibition of transcription factors and, optionally, self-reinforcing positive feedback loops for transcription factors. Similar models have been proposed in the context of hematopoietic differentiation Laslo2006 ; Kueh2013 , and they provide a quantitative framework to account for the bimodal distribution of transcription factors in immune cells. How do these models compare with experiments?

A theoretical inspection of the Th1/Th2 model of T-cell differentiation by Antebi et al. Antebi2013 highlighted how the bistable solution of the model presented in Figure 7 does not generically produce a bistable solution. A simple scan of parameters within the physiological range for induction of Tbet and GATA-3, the candidate $P_{1}$ and $P_{2}$ in Eq. V.1.1, demonstrated that the ratios of rates for their self-reinforcement and inhibition most commonly predicted a single stable fixed point, and only rarely more than two. In addition, the induction of transcription factors is known to be inherently stochastic, because of intrinsic noise in gene induction Walczak2005b ; Walczak2005c ; Friedman2006 ; Elowitz2002 or extrinsic noise in the levels of expression of signalling molecules in T cells Feinerman2008 . Overall, the “simplest” model of a toggle switch in T-cell differentiation, though elegant in its inception, is not consistent with observations. This case is typical in quantitative immunology, where initial insight is often derived from genetic manipulations with severe consequences as in the case of the Th1/Th2 system where knocking out Tbet leads to a reduction/maintenance of induction of Th1/Th2 cytokine secretion, respectively Szabo2000 . These knockout experiments must be balanced and reinterpreted in nonperturbative settings, e.g. by monitoring the spontaneous induction of Tbet in differentiating T-cells. Similar reassessments are happening in many models of leukocyte differentiation, such as in the M1/M2 model of macrophage differentiation Palma2018 , or PU-1-controlled differentiation in the myeloid compartment Kueh2016 ). Future work will require better parametrization an innovative modeling approach that can embrace the combinatorial and dynamic complexity of cytokine communications and gene regulation in the immune system.

Such “murkiness” in immunology must be embraced as it opens up functional possibilities. For example, Peine et al. Peine2013 tracked the formation of Tbet+GATA-3+ mixed phenotypes, i.e. with characteristics of both Th1 and Th2, so in apparent contradiction with bistability. Their presence was demonstrated to limit the deleterious impact of all-out, non-mixed Th1 or Th2 inflammation. Thus, a better understanding of noise and stability in immune cell differentiation will be key to understanding virtuous and pathologic inflammation.

V.2 Hematopoiesis

V.2.1 Timescales

Where do immune cells come from? A human being has $\sim 5\cdot 10^{13}$ cells in the body and more than half of them are made by hematopoietic stem cells (HSC) found in the bone marrow. There are an estimated order of magnitude $\sim 10^{16}$ cell divisions per human, which gives $\sim 10^{6}$ cell divisions per second  Moran2010 ; MiloPhillips , and most of them are linked to the HSC cells. HSC number $\sim 10^{4}$ in mice (there are no reliable numbers for humans), many of them already made in the fetus. Through a series of differentiation and phenotypic commitment events (HSC$\rightarrow$ Short Term HSC (ST-HSC) $\rightarrow$ Multipotent progenitor (MPP), these cells give rise to different kinds of cells found in the blood, including red blood cells and all immune cells (see Fig. 9). The first branching decision, whose precise timing is currently being questioned, is about becoming a myeloid progenitor, that will give rise to red blood cells, mast cells, thrombocytes (so not immune cells) but also macrophages, granulocytes (neutrophils, basophils and eosinophils) and dendritic cells – so cells of the immune system that eat up cells and proteins non-specifically. They function as cells of the innate immune system, simply eliminating bacteria and other pathogens, but some of them (e.g. dendritic cells) also play an important role in adaptive immunity as antigen presenting cells (APC). The remaining branch of the differentiation tree leads to lymphoids that include cells of the adaptive immune system: B and T-cells as well as natural killer (NK) cells. The decision about becoming a myeloid or a lymphocyte is widely assumed not be an autonomous decision but is believed to be influenced by external signal (see below for a detailed discussion).

General quantitative questions about differentiation apply to the differentiation of cells in the immune system. For example, given that in mice there are of the order of $10^{4}$ HSC, what kind of dynamics results in $\sim 10^{7}$ very short lived granolocytes (of the order of a day) while at the same time producing $\sim 10^{11}$ lymphocytes that can live even for years? How long does it take to produce these cells? What kind of differentiation process can produce this kind of diversity, and how is it regulated to produce the right numbers of cells? A lot of information about these differentiation processes has be gained from two types of experiments. The first involves exposing cells at a given upstream stage with died or radioactive markers (i.e. bromium, deuterium) that get taken up — and then diluted — upon cell division. Analysis of the decay curves is informative about cellular lifetimes. The second type of experiment — an adoptive transfer experiment — is experimentally harder, since it involves transplanting new marked cells into an animal (typically mouse) and then tracking them. To study cells during the early stages of hematopoiesis the mouse must first be cleared of its natural cells, and recent results suggest that the dynamics of differentiation after transplantation may be very different from regular dynamics Busch2015 . Lymphocyte divisions can be studied without killing the host’s own immune system. The principle of these experiments is simple — after some time $\delta t$ the marked cell will appear in a given compartment.

At the population level, given that cells at stage $i$ of differentiation can proliferate with rate $\lambda_{i}$, differentiate with rate $\alpha_{i}$ or die with rate $\delta_{i}$, the dynamics of a given cell type $n_{\rm R}$ in the differentiation process follows Hofer2016 ; Busch2015

[TABLE]

where $n_{\rm u}$ describes the upstream (pre-differentiation) cell type. Cells can also migrate, which introduces a spatial component to the equations. In general, all these rates are functions of the concentrations of the different populations $n_{i}$, resulting in non-linear equations. $\alpha_{\rm R}+\delta_{\rm R}-\lambda_{\rm R}=k_{\rm R}$ define a collective timescale of decay of the $n_{\rm R}$ population. This highlights the identifiability problem that without specially planned out experiments guided by theory, it may be hard to tease apart the timescales of biological interest.

Commonly in this kind of analysis it is assumed that the cellular differentiation of the adult individual is in steady state. We will return to this assumption when we learn the timescales for the process. Within this assumption the steady state ratio of cells in the compartment of interest, compared to the upstream compartment is ${n_{\rm R}}/{n_{\rm u}}=\alpha_{\rm u}/(\alpha_{\rm R}-\beta_{\rm R})$, with $\beta_{\rm R}=\lambda_{\rm R}-\delta_{\rm R}$. This kind of approach has been used to learn proliferation rates for T-cells Zilman2010 . Stochastic versions of such models have also been considered in detail Perelson2002 ; Yates2007 ; Zilman2010 , showing that theories of age-structured cell populations Perelson2002 are equivalent to branching processes Zilman2010 (see section XII.3.6). For hemapoietic stem cells, a HSC in the bone marrow divides once every 50 days ($\lambda_{R}=1/110{\rm days}$) and a HSC becomes a granulocytes after $\sim 1$ year (and then lives for a day). This is known from so-called fate-mapping experiments, consisting of tracking labeled cells, which is easier for granulocytes since they are known not to proliferate. It is also known that granulocytes have no other cell source, unlike other cell types. In this case the label from the HSC appears in a granulocyte after a year. Since a lifetime of a mouse is $\sim 2$ years, a $1$ year timescale to produce a new cell may put the steady state assumption into question. The fetal dynamics is very different, the steady state assumption cannot be used, and a new granulocyte is produced in a few weeks Busch2015 .

Technically, measuring cell counts directly is not reliable because of sampling bias. Instead, one should measure cell ratios, $f_{i}=n^{i}/n^{\rm total}$, as was shown in the previous example. In steady state, the fraction of cells of a given type evolves according to $d{f}_{i}/dt=k_{i}(f_{i-1}-f_{i})$, with $k_{i}=\alpha_{i}-\beta_{i}$. By looking at subsequent compartments it is possible to disentangle differentiation and proliferation. This gives an estimate of how self-renewing the compartment is ($\beta_{i}$), and how much comes from upstream differentiation ($\alpha_{i}$). For a perfect stem cell we expect $k_{i}\rightarrow 0$. Real stem cells (e.g. HSC) are essentially self-renewing with $\tau$ (in mice) on the order of the lifetime of the mouse. Other differentiation compartments are clearly transitional, with $k_{i}^{-1}\sim 1$ day. HSC divide once every $110$ days but they also differentiate once every $110$ days, so they seem to have a lifetime of (they replace themselves every) $110$ days (this remains to be verified with more direct means). Short Term HSC (ST-HSC, see Fig. 9) die at a rate of 1 per month and Multipotent progenitor (MPP) of 1 per day. This produces a strong self-renewal gradient from HSC to MPP. At the same time the proliferation rates of these cells go up from HSC to MPP with a similar gradient. Some have speculated that these inverse trends of proliferation and self-renewal protect the organism from malignancy by making sure the right amounts of cells are produced. Myeloid lymphomas — a condition where the organisms produces too many myeloid cells — have been linked to disturbances in this balance and the condition is only visible in aged individuals when the overproduced cells have accumulated.

Most HSC become myeloid cells, with only one in 100 cells becoming a lymphoid cell. Lymphoid cells are very long lived and the thymus can function for long time without input from new HSC. In general, due to the slow timescales of the whole process, the system has a lot of inertia. If all HSC are killed in a young adult mouse, there are no visible phenotypic effects for $\sim 5$ months, unless the hemopoietic pathway itself is stressed Schoedel2016 .

There are a priori no biological reasons for downstream cells in the hematopoietic pathway not to influence upstream cells, however this has not yet been directly observed. Probably certain feedback mechanisms are at play.

V.3 Inferring the hematopoiesis differentiation tree

Analyzing experiments from mice with progenitor cells carrying unique barcodes with a maximum likelihood approach (see section XII.5.1) coupled to a stochastic model, Perié et al Perie2014 , uncovered a more complex picture of myeloid and lymphoid differentiation. They considered the lymphoid primed multipotent progenitors (LMPPs) differentiation process during which each LMPP is called an MDB, since it has the potential to become a myeloid (M), dendritic (D) or B cell. They built a stochastic branching process model of cell differentiation (see section XII.3.6), where in each step of the decision tree a cell can loose one or two of its potentialities: an MDB can become an MD, MB, DB or directly an M, D, B. Similarly an MD can become an M or a D etc. They contrasted this with a traditional differentiation model, where the MB phenotype is forbidden, and an MDB cannot loose two potentials at once. Using mice with progenitor cells carrying unique barcodes, they build lineage trees, traced each offspring cell to its parent and counted how many offspring of each type ($N_{i}$ for $i\in\{M,D,B\}$) were produced by each MDB cell. The generating function $G(z_{M},z_{D},z_{B})=\sum_{m,d,b}z^{m}_{M}z^{d}_{D}z^{b}_{B}P(m,d,b)$ for the probability of a cell to give rise to a given combination of offspring cells $P(N_{M},N_{D},N_{B})$ can be calculated from a convolution of simple branching process generating functions Perie2014 ):

[TABLE]

where $H(p,z)$ is given by the generating function, in Eq. 130 in section XII.3.6, and

[TABLE]

$p_{\rm MDB},p_{\rm MB},p_{\rm MD},p_{\rm BD}$ are the probabilities of not differentiating from (staying in) the $MDB,MB,MD,BD$ states and $p_{MDB\rightarrow ij}$ describes the transition probabilities from the $MDB$ state to the $ij\in\{MD,MB,DB\}$ states. The remaining functions are defined recursively as in Eq. 39, e.g. for $ij=MD$: $G_{\rm MD}(z_{M},z_{D})=p_{MD\rightarrow M}z_{M}+p_{MD\rightarrow D}z_{D}$ etc.

The marginals of this distribution give the probability that a given barcode appears in each of the 7 possible cell types. Following the equal loss of potential model (ELP) Ogawa1983 ; Tsuji1989 , the rate of loss of the cells’ ability to make different cell types was assumed to be constant and independent of current cell type, which gives three rates $\alpha_{M},\alpha_{B},\alpha_{D}$, that parametrize the transition probabilities between the cell types $p_{i\rightarrow j}$ (for all allowed combinations of $i,j\in MDB,MD,MB,DB,M,D,B$) and reduces the number of parameters. The parameters are determinied by carrying out a maximum likelihood fit $\vec{p}_{\rm max}={\rm argmax}_{\vec{p}}\left(\sum_{i}n_{i}\pi_{i}(\vec{p})\right)$ , where $n_{i}$ is the number of observed barcodes in state $i$ and $\pi_{i}(\vec{p})$ is the model probability (from Eq. 38) that an ancestral barcoded cell gives rise to offsprings in states $i$, where $i\in MDB,MD,MB,DB,M,D,B$ calculated from the generating function: $\pi_{MDB}(\vec{p})=\sum_{m,d,b\geq 1}P(m,d,b)$, $\pi_{MD}(\vec{p})=\sum_{m,d\geq 1}P(m,d,b=0)$, $\pi_{M}(\vec{p})=\sum_{m\geq 1}P(m,d=0,b=0)$ etc.

This approach showed that the classical model of sequential loss forbidding the MB state is not consistent with data. In fact, the MB state is very unlikely but is essential to explain the observed barcode distribution. Simulation of the dynamics of these models showed that the differentiation process between LMPPs and the committed cell types requires $\sim 20$ rounds of differentiation, which translates into about $2$ weeks. Recent barcoding experiments combined with similar inference methods that looked at earlier stages of differentiation from HSC cells Pei2017 showed that most HSC cells give rise to many cell fates.

Similar techniques of branching processes coupled to likelihood inference methods were initially used by Yates et al Yates2007 , to analyse fluorescent dye experiments of differentiating T-cells and estimate cell division and death rates.

V.4 Cell fate during the immune response

V.4.1 Choice and timing of cell fate under stimulation

Lymphocytes, both T and B-cells, start growing and dividing upon stimulation during an infection. When they divide they also acquire new specialized functions, make many decisions at the individual cell level, e.g. whether to divide further, whether to switch class for immunoglobulins (Ig) molecules, whether to become memory cells for both T-cells and B-cells, or to become plasmablasts for B-cells. As their number increases, they eventually make another decision to stop dividing and start dying. At the population level we observe a decrease of the population size, after the large expansion due to proliferation, and a return to pre-infection cell counts. Since each cell can in principle adopt a different decision path, at the population level we see a large combinatoric diversity of cell states: from cells that have divided a small number of types and become memory cells, to cells that have divided many types to leave no offspring.

This diversity is even more staggering for B-cells, since they undergo a phenomena called class (or isotype) switching recombination (CSR), where the variable region (discussed in Sec. VI) of the antibody does not change but the constant region of the receptor is modified. This process does not change the affinity of the antibody for the antigen, but it means the antibody cannot interact with signaling molecules (see sec. III.2). B-cells do this by expressing a specific gene located in the heavy chain locus, organized in a specific order. During CSR, segments of genes in this locus are removed, and the remaining DNA is recombined to encode a different isotype. Since the double stranded breaks occur at conserved nucleotide motifs, the identity of the expressed isotype-encoding genes are conserved between cells and individuals. Since the non-expressed isotype-encoding genes are deleted from the locus, in general the order of cycling through the classes is conserved, although inter-chromosomal translocation from the other allele can add deleted loci Laffleur2014 .

How are these decisions made by cells? One view is that cells integrate signals from the environment, which act as cues to trigger decisions. While this is certainly true Murphy2007 ; Murphy2000 , and experimentalists can induce certain cell fates in vitro using cytokine cocktails (see sec. IV), an alternative (yet not contradictory) idea was proposed in the “cyton” model Hawkins2007 . The basic idea behind this model is that each cell in the population makes a stochastic decision, choosing from one of the accessible cell fates. To do this it uses an intrinsic pre-programmed timer, telling it when it should divide and die. However since the timer in each cell picks a time from a distribution of times, each cell will go through a different decision scenario. Putting together the different decisions each cell can make this leads to a heterogeneity of cell fates in one population. Yet, for each trait (division times, becoming a memory cell), the population level distribution is reproducible.

The model can be and has been adapted to different cells Hawkins2007 ; Marchingo2014 ; Marchingo2016 ; Duffy2012 but let us present it on the simple example of three possible options: to divide, to die or do nothing even if stimulated. Given an experimentally characterized distribution of division time $\phi_{i}(t)$ and death times $\psi_{i}(t)$ for each division round $i$, and assuming that a fraction $F_{i}$ of cells will divide when stimulated in division round $i$, the number of cells dividing or dying per unit time is a given round $i=1,...,m$ of divisions is:

[TABLE]

where the terms in parantheses account for cells that would have divided (died) at time t but had previously died (divided). The total number of cells in each division round is then calculated by summing gains and losses in each round:

[TABLE]

for $i=1,...m$. These integrals can be performed analytically for exponential distributions but since the experimental data are better described by log-normal distributions Deenick2003 ; Tangye2003 , the integrals are performed numerically.

The cyton model and its generalizations have been shown to be able to fit the number of cells in experiments with and without stimulation Hawkins2007 . Of course, one could note that the cyton model takes the experimentally observed distribution of division and death times and simply calculates the number of surviving cells. Yet the power of the model lies in the paradigm shift that assumes there are rules and, while the immune system functions stochastically and relies on heterogeneity, we can nevertheless predict and understand its behaviour. The idea of a programmed stochastic system, as opposed to a black box integration of cues, is extremely powerful.

Using data from transgenic mice infected with the influenza virus, Marchingo et al Marchingo2014 showed that T-cell receptors and costimulatory signals impose an intrinsic division fate on each cell, and cells then count through generations, as defined by the cyton model, before returning to the pre-stimulation state. This initial heritable priming can be later modified by dose dependent cytokine signling, which is also integrated by the cells and influences their the response. As a result different combinations of costimulatory signal and cytokines can generate signals of similar magnitude. These experiments also show that the cyton model is not in contradiction with signaling-based decision making.

V.4.2 An aside on cell types

Cell types in immunology are defined using surface markers. Cells communicate by binding and unbinding signaling molecules (cytokines), for which they need specialized receptors. A cell that commits to participate in a given communication channel, expresses a surface receptor and the combination of surface receptors gives the spying experimenter an idea of how this cell is bound to behave when triggered – this is what we call a cell type. This phenomenological approach has been very successful in the history of immunology and provides a way to classify different cells and give us some idea about their function and properties. In practice, FACS (Fluorescence-activated cell sorting) sorting experiments segregate cells according into ones that have a high concentration and a low concentration of a given surface marker on their surface. By doing this in many dimensions, one can zero in on a very specific cell type. Technically, this leaves the problem of deciding where the boundary between high and low is – in practice cells are heterogenous and the experiment produces a distribution of marker concentrations, which, if we are lucky, is bimodal. This problem of finding the separatrices between cell types is called gating, and is an art implemented in analysis softwares. Advances in machine learning will alleviate the need for expert-guided gated and make the definition of cell types more automatic and easy-to-validate. Being aware and using the heterogeneity of the population data is often useful. Looking at the population distribution has made people realize that often there is continuum of phenotypes and a binary approach is not valid. This has been amplified by recent high throughput cell sorting experiments by mass cytometry –so-called CYTOF Bendall2012 , which make it clear that cell types are more continuous and dynamic than traditionally assumed. With that in mind, we can now try to learn something about how cells acquire and switch identity.

V.4.3 Inferring cell fate timelines during the immune response

The problem of infering the timeline of cell differentiation during T-cell response can also be treated using the cyton theory Hawkins2007 ; Marchingo2014 ; Marchingo2016 ; Duffy2012 . Here we present another more data-driven approach to the problem described in Buchholz, Flossdorf et al Buchholz2013 . We will concentrate on the example of CD8 T-cells. A clone, here defined as all the TCR that respond to the same antigen, starts with $\sim 10-100$ long-lived naive cells (lifetime of about 80 days in mice). During an infection this clone grows to $\sim 10^{7}$ short-lived effector cells. After infection clearance, there remains $10^{3}$ memory cells with an intermediate lifetime (about 15 days in mice). These memory cells can then respond more robustly in a subsequent infection. Memory and effector cells also proliferate within their classes to maintain their pool. This is done through small antigenic signals they constantly receive for memory, and strong antigenic signals for effector cells.

The above textbook description provides a simplified picture, but is unable to put measurable numbers behind each of the described populations. More importantly, this picture was painted by looking by eye at bulk population data summing over many cells. Can inference and data discrimination approaches help us infer the differentiation pathways? Population-level adoptive transfer experiments can be used to infer the parameters of known models, but it is much harder to infer the topology of the network, especially from single snapshots. Flossdorf et al Buchholz2013 ; Flossdorf2015 analyzed single cell adoptive transfer (single cell fate mapping) experiments to discriminate between different differentiation tree topologies and find the underlying model parameters. They considered a model with four cell types: naive (N), effector (TEF), effector memory precursor (TEMp) and central effector memory precursor (TCMp) and wrote down the most general network diagram, in which every non-naive cell type could differentiate into any other cell type, all of which also could proliferate and further differentiate. Naive cells could become any of the other cells, but no other cell could become a naive cell. This led to 302 possible models, defined by their network topologies, as well as their parameters. The precise definitions of the particular cell types in the experimental data was done using surface markers (see discussion above) can be found in the original paper. The immunological meaning of these cell types is under debate – but we are mainly interested here in the power of the theoretical method.

The considered models where simple stochastic dynamics for each cell type that considered differentiation ($\vec{\alpha}$) and proliferation ($\vec{\beta}$) reactions. As mentioned above, a lot of heterogeneity was observed within naive cells and in vitro experiments show heterogeneity in the reaction rates, which had implications on the model. Differentiation rates where chosen from underlying distributions: an exponential and gamma distribution of differentiation rates were considered. The data resolution did not allow to discriminate between the two, with an exponential giving a sufficiently good fit. The proliferation rates were taken to be fixed, since longitudinal correlations of proliferation events were previously shown to decay rapidly Dowling2014 . The considered models ignored the heterogeneity in the phenotype of naive cells, but since the heterogeneity in the number of naive cells would render model selection impossible, the experiment was tuned to start with one naive cell ($N_{\rm naive}=1$).

The model consisted of a Master equation for the evolution of the joint probability of the different cell types:

[TABLE]

where is F a linear operator describing the proliferation and differentiation dynamics of the different cells types, $\vec{N}=\left(N_{\rm TEMp},N_{\rm TCMp},N_{\rm TEF}\right)$ are the numbers of cells of each type, and $\vec{r}$ is the death rate vector of all cell types. The analytically calculated first and second moments of these equations are fit to the variances ($CV_{i}(t)$) and covariances ($\Sigma_{ij}(t)$ and mean concentrations ($\langle{N}\rangle{}_{i}(t)$) for each cell type $i$ from cell fate mapping data using a mean squared error minimisation. The model selection is then finalized by comparing the AIC (see section XII.5.1) of the best fit for each topology. The necessity to use correlations in the inference procedure, and not just population averages, was demonstrated on synthetic data generated with known parameters and topologies. For instance, both a linear ($A\rightarrow B\rightarrow C$) and tree ($C\leftarrow A\rightarrow B$) topology can give rise to the exact same evolution of the mean number of cells as a function of time, for well chosen parameters. However covariances between these numbers give clearly different signatures, allowing us to distinguish the different topologies. These data-based, model-selection approaches are similar to those developed for stem cell hemapoiesis Perie2014 , described in section V.3.

As a result, from the 302 possible models, two strongly stand out as potential candidates to describe the cell differentiation of CD8 T-cells responding to a listeria expressing chicken ovalbumin (to which these T-cells were specific) in mice. These two models are:

[TABLE]

where each cell type except for the naive one proliferates. The models have well defined parameters that can easily be identified. This bacterial response is an acute infection and generates exponential cell growth, and the inferred proliferation rates form a gradient from ${\rm TEF}\rightarrow{\rm naive}$ cells with naive cells proliferating the least. Both of these models correctly predict the time dependent dynamics of the mean concentrations of the different cell types. Even if learned on single-cell progeny data at only one time point, the first linear diversification model is able to predict the phenotypic composition of the expanding population at earlier time points. These kinds of approaches shed light on the possible differentiation dynamics and allow us to rule out possible scenarios. To give just one example, for this system the cell differentiation is largely ($\sim 90\%$) symmetric: the two offspring of a cell are usually the same, although not the same as the mother, as opposed to asymmetric divisions where the two daughter cells have different cell types which happen about $\sim 10\%$ of times. Also, it turns out that the differentiation and proliferation rates must depend on time to explain the data. These time-dependent rates are the same in a response to both bacterial (listeria) and viral challenges, suggesting universal response dynamics.

Using similar methods to those proposed by Flossdorf et al Buchholz2013 ; Flossdorf2015 , Miles et al Miles2018 showed that the model is consistent with the data of Buchholz Buchholz2013 et al, but applied to the data of Kinjyo et al KinjoNatComm2015 suggests that memory precursors are produced before effector cells. This shows that this type of inference is robust but the biological interpretation of the results depends both on the experimental conditions and model assumptions.

V.4.4 Quorum sensing

Recent experiments have revisited the question of naive T-cells differentiation into different memory cell types upon antigen stimulation, by considering collective effects through cell-cell communication. Polonsky et al Polonsky2018 tracked live cells tagged with antibody markers over time in microwells with different initial cell numbers. By continuously tracking individual cell differentiation states and proliferation, they are able to overcome the limitations of bulk experiments where local cell densities are hard to control. Live-cell imaging showed that the decision to differentiate into progenitor central memory (pTCM) is collective: in equal medium conditions, cells surrounded by more cells had a larger differentiation rate. That rate depended solely on the instantaneous number of cells in the well, as shown by data collapse of the differentiation rate onto a single curve as a function of the number of cells, for various conditions (different media, different initial number of cells in the well). The curve shows that cells are more likely to differentiate into pTCM if there are more than $N_{c}\sim 30$ cells in the cluster. The process can be described by a stochastic differentiation model, in which cells divide and die, and differentiate with a rate $R$. However, to explain the data $R(N)$ must be made to depend on the number of cells $N$. Such a stochastic model also captures the observed well-to-well variability. The collective decision could be linked to the IL-2 and IL-6 cytokine communication: blocking IL-2 reduced the maximum of the universal differentiation rate curve, while blocking IL-6 increased $N_{c}$ without significantly altering the plateau value. Interestingly, IL-2 and IL-6 receptors where found to cluster on cell surfaces, with receptor patches directed towards neighbouring cells. This experiment shows that the percentage of cells that differentiate into memory cells and effector cells depends on the instantaneous local T-cell density, independently of the additional post-stimulation influence of T-cell receptor signalling strength, or the effect of antigen presenting cells.

Polonsky et al Polonsky2018 give a threshold in terms of absolute number of cells in a micro-well needed for collective decision making, $N_{c}\sim 30$. However, that number must depend on the size of the well, if cells do communicate through cytokine diffusion. As we have discussed in section IV.2 in the context of CD4+ T-cell differentiation into Tregs and T-helper cells, signal propagate to short distances as they are taken up by neighboring cells Oyler-Yaniv2017 . For pTCM differentiation, this is facilitated by the IL-2 and IL-6 orientation of the receptors towards secreting cells. But at large distances, diffusion starts playing a role. Active clustering seems to be a mechanism that puts cells in the non-diffusive regime.

The scale of cytokine communication has also been shown to depend on the type of signal: an anti-inflammatory IL-10 signal produced by dendritic cells is long-ranged, while the pro-inflammatory TNF$\alpha$ signal produced by the same cells is short ranged. These cells are also known to cluster in cultures and in-vivo. A two step model for the integration of these two opposing signals was proposed to explain the nonlinear response of dendritic cells when integrating pro- and anti-inflammatory information Marcou2018c : an initial intermediate molecule integrates both signals in an indiscriminate way until a threshold value, providing a bottleneck. This way either the pro- or anti-inflammatory signal can stimulate the system. But the output of the this bottleneck is later down-regulated by only the anti-inflammatory signal. This combination of collective decision making with a modulated bottleneck safeguard was proposed as a way to control possibly excessive immune responses.

Collective decision making, also called “quorum sensing”, was previously proposed theoretically as a way for immune cells to help solve the self non-self discrimination problem Butler2013a . According to that theory, which we discuss in more detail in Sec. VI.3, communication allows cells to make a census-based decision by integrating the signal read out by many clones in order to correct mistake made by individual cells. This is a cellular implementation of an error-correcting code, similarly to kinetic proof reading discussed in section III.1.1. The context of that proposal is different from that of Polonsky2018 , where the T-cell population is monoclonal, while the idea of error correction relies on taking a census of different T cell clones that make independent decisions. This difference in clonality does not change the general similarity in the nature of the collective decision. The experimental results are likely to hold for different T-cell clones with with similar affinities for the stimulating antigen. The T-cell response is multiclonal (see section VI.7) and any quorum-sensing mechanism is likely to involve multiple distinct clones. Understanding the details of how the memory and effector cells pools are controlled in a multiclonal setting remains an open question.

VI Repertoires

VI.1 Size of immune repertoires

For the adaptive immune system to protect us against all the different pathogens we may encounter, including ones that may not exist when we are born, we need a large set of different immune receptors.

There are of the order of $4\cdot 10^{11}$ T cells circulating in the human body Jenkins2010 , and of the same order of B cells Glanville2009 , each expressing a single type of receptor to a first approximation Han2014 ; Dupic2018 . However, T and B cells divide, meaning that many cells can carry the exact same receptor, defining a “clonotype”. The repertoire size, or number of clonotypes, is therefore smaller than the number of cells. Early work Arstila1999 based on a subset of the repertoire gave an estimate of $\sim 10^{6}$ unique $\beta$ chains in one human individual.

However, these estimates were indirect, and have been updated recently thanks to the advent of repertoire high-throughput sequencing (RepSeq) techniques in 2009 Weinstein2009 ; Robins2009 ; Boyd2009a , discussed in a number of recent reviews Benichou2012 ; Six2013 ; Robins2013a ; Georgiou2014 ; Heather2017 ; Minervina2019 . This method, which we briefly discuss in section XII.6, allows one to obtain the sequences of a non-exhaustive but fairly substantial portion of all receptors in a biological sample. As different technologies have been emerging, the need for standarisation of RepSeq data reporting has appeared, and the Adaptive Immune Receptor Repertoire (AIRR) community NatureAIRR has been organizing meetings and producing web-based and journal publications to promote data reporting unification.

RepSeq experiments on human blood samples, combined with statistical extrapolation estimators, report unique TCR$\beta$ numbers ranging from $4\times 10^{6}$ Robins2009 (using Fisher’s Poisson abundance model Fisher1943 ) to the order of $10^{8}$ Qi2014 (using the Chao2 estimator Chao2002 ). For BCR heavy chains, estimates obtained with the Poisson abundance model yields $1$-$2\cdot 10^{9}$ DeWitt2016 . One should take all these empirical estimate with great caution, as they make unverified, and probably wrong, assumptions about the shape of the clonotype abundance distribution when doing the extrapolation from small blood samples to the entire organism (more on that in Sec. VI.4). Theoretical estimates based on population dynamics models assuming that naive cells divide little or not at all logically give much higher estimates, closer to the number of cells itself, $10^{10}-10^{11}$ Lythe2015 ; Mora2016e .

How large is large enough? De Boer and Perelson Boer1993 proposed that the minimal repertoire size is ultimately determined by self-tolerance. We revisit their argument in the light of modern estimates, with simplified notations, starting with the easier case of B-cells. B-cells interact directly with antigens, unlike T-cells that interact with short peptide fragments presented on the multi-histocompatibility complex (MHC). Call $p\ll 1$ the probability that a randomly chosen receptor recognizes a randomly chosen epitope. During the self-tolerance selection process, the initially generated repertoire of size $R_{0}$ is reduced to $R$, $R=fR_{0}$, where $f$ is fraction of receptors that survive negative selection, i.e. that do not recognize any self-epitope, so that

[TABLE]

where $n\sim 10^{5}$-$10^{6}$ is the number of self-epitopes (in humans). The probability that an epitope escapes the immune system is then given by

[TABLE]

This gives the repertoire size $R=-(1/p)\ln p_{e}$ as well as the pre-selection repertoire size

[TABLE]

Assuming that evolution has optimized the recognition probability $p$ so as to minimize $R_{0}$, solving for $\partial R_{0}/\partial p$ an optimal $p^{*}=1/n$, and an optimal repertoire size

[TABLE]

with an optimal selected fraction $f^{*}=\exp(-1)\approx 0.37$.

This estimate suggests that the number of self-epitopes is the main determinant of the minimal repertoire size, with which it scales linearly. Another prediction is that the recognition probability should be inversely proportional to the number of self-epitopes. We can intuitively understand this result by thinking about the repertoire as a covering problem: viewing the set of epitopes recognized by a single receptor as a ball in an abstract phenotypic space, how many such balls does one need to cover the entire space of foreign epitopes, while avoiding a finite number of self-epitopes? The simple argument given above, as well as the schematic of Fig. 10, tells us that the volume of each ball, $p$, should be inversely proportional to the number of self-epitopes, $n$, implying in turn that their number ($R$) should scale with $n$.

Empirical estimates suggest $p\sim 10^{-5}$ and $n\sim 10^{5}-10^{6}$ Boer1993 , consistent with the theoretical prediction $p\sim 1/n$. The prediction of $f\approx 37\%$ of cells passing negative selection is consistent with recent experimental estimates $25\%-45\%$ Wardemann2003 . The predicted B cell receptor diversity crucially depends on $p_{e}$, which is hard to estimate. De Boer and Perelson estimated that 99% of antigens, presenting each $10$ epitopes, were recognized by at least one antibody of the immune system. This implies $p_{e}^{10}\sim 10^{-2}$, and $R^{*}\sim 5\cdot 10^{4}--5\cdot 10^{5}$, which is much smaller than current experimental estimates of BCR diversity ($\sim 10^{9}$).

These estimates become a bit more involved for TCR since epitopes are presented by the Multi-Histocompatibility Complexes (MHCs), which translates into an additional level of sampling. As already mentioned in previous sections, there are 2 types of MHCs relevant to adaptive immunity: MHC type I complexes present peptides from inside the cell to CD8+ (also called killer) T-cells, whose goal is to kill infected or cancerous cells; MHC type II complexes present peptides from outside the cell to CD4+ (helper) T-cells whose role is to stimulate B-cells during affinity maturation process (see section below).

Each epitope, or peptide derived from the antigen of interest (self or foreign) must both be presented by one of the $m=6$ MHC genes for each MHC type, and be recognized by one of the $R$ TCR expressed in the body. This implies that the probability of escape reads:

[TABLE]

where $q\ll 1$ is the fraction of peptides that can be presented by a given MHC molecule. Eq. 43 is still valid, with $n$ being the total number of possible MHC-self-peptide complexes. The repertoire size is then given by $R_{0}=-({e^{pn}}/{p})\ln[1+({\ln p_{e}}/{mq})]$, and the optimal size is obtained through $\partial R_{0}/\partial p$:

[TABLE]

with $p^{*}=1/n$ and $f^{*}=\exp(-1)$. Thus, the scaling and selection probability are predicted to be the same as for B cells.

The probability of presenting a given peptide on MHC class I complexes is relatively well known: deep learning algorithms have been very succesful in building algorithms that predict which peptide is likely to be presented. These methods are implemented in software packages such as netMHC Buus2003 ; Andreatta2015 , which account for known biological features to guide the learning. The same prediction task has been implemented for MHC class II Jensen2018 , but the prediction is much harder, partly because unlike MHC class I presented peptides, which have a relatiely fixed length of 9-11 amino acids, MHC class II presented peptides have varying length, from 13 to 25. More data than is currently available is needed to successfully train the neural networks, and to improve prediction. The models predict that about 1% of peptides are presented by any given particular MHC allele, i.e. $q=0.01$.

The number of MHC-self-peptide complexes is given by $n=qmn_{\rm p}$, where $n_{\rm p}$ is the number of peptides in the human proteome. Each individual has 6 MHC alleles for each class, i.e. $m=6$. $n_{\rm p}$ is roughly the number of proteins, $\sim 3\cdot 10^{4},$ times their average length, $\sim 400$, i.e. $n_{\rm p}\sim 10^{7}$, and $n\sim 6\cdot 10^{5}$. The recognition probability, $p$, has been estimated to be in the range $10^{-4}$-$10^{-6}$ Moon2007203 ; moon-2012 . These numbers are consistent with the theory $p\sim 1/n$. The fraction of negatively selected cells, predicted to be $f^{*}\approx 37\%$ by the theory, is estimated to fall in the range $20$-$50\%$ Yates2014 . As for B cells, predicting the size of the repertoire is hard without a good estimate of $p_{e}$, but current estimates of $R=10^{8}$-$10^{10}$ are probably much larger than predicted by (48) regardless of the estimate of $p_{e}$.

In summary, these theoretical estimates still seem to be relevant today for the relationship between recognition probability, $p$, and number of self-epitopes, $n$, both for T and B cells. Recent estimates of TCR and BCR diversity, however, are much larger than predicted by these simple theories. However appealing, a major limitation of these estimates and predictions is that Eq. 43 assumes that each receptor is tested against all possible self-epitopes, while in practice this is impossible due to the limited duration of lymphocyte maturation. We will come back to this point in Section VI.3 devoted to modeling thymic selection.

VI.2 Inference of the stochastic repertoire generation process

Antigen receptors are proteins, which must be encoded as genes in the DNA. Humans have of the order of $10^{4}$ protein-coding genes. Directly encoding the whole diversity of immune receptor genes ($\sim 10^{8}-10^{10}$) in each genome would make it impossible for the DNA to fit in the nucleus. The immune system has solved that problem by stochastically creating receptors in each cell, in a process called V(D)J recombination that combines combinatorics and randomness to generate diversity. B- and T- cell receptors are made of two chains, light and heavy for BCR, and $\alpha$ and $\beta$ for TCR. The genome encodes a certain number of gene templates, called V, D and J for the heavy and $\beta$ chains, and V and J for the light and $\alpha$ chains. Upon creation of a chain, DNA is edited and one of each of these gene templates is chosen per receptor. The combinatorics of templates typically results in $\sim 10^{3}$ different receptors. To obtain the (much larger) observed diversity, nucleotides are randomly inserted in a non-templated way and deleted at the junctions between the V and D and D and J genes (or V and J genes for $\alpha$ and light chains). This process of generating “junctional diversity” in fact accounts for most of the diversity of the repertoire Murugan2012 ; Elhanati2015 .

Characterizing the above described process of receptor generation in quantitative detail has been made possible thanks to the development of RepSeq methods (see Sec. XII.6). This can be done using out-of-frame sequences (with a frameshift due to the random number of additions and deletions at the junctions), which are nonproductive and hence a raw product of recombination, as they are free of selection effects. These sequences survive selection because they are in the same cell as an in-frame sequence. Since each cell has two chromosomes, out-of-frame sequences come from cells where the first rearranged chromosome resulted in an out-of-frame sequence and the second rearranged chromosome resulted in an in-frame sequence. Due to the random insertions and deletions of nucleotides, it is impossible to reliably determine how a given receptor was formed (V, D, J assignments, as well as distinguishing inserted from templated nucletoides) from its mere sequence. Instead one can consider a list of scenarios (which includes V, D and J gene choice plus a number of insertions and deletions at each of the junctions) and sum over them weighted by their probabilities as determined self-consistently using a probabilistic model. Concretely, in the case of the simpler $\alpha$ or light chains, the probability of a given recombination scenario $r$ is:

[TABLE]

where $\text{del}V$ and $\text{del}J$ denote the number of deletions at the V and J ends, and ${\rm ins}VJ$ is the list of inserted nucleotides. A similar expression can be written for the $\beta$ or heavy chains by adding the D gene and the related deletions and insertions. The model is the most general factorizable distribution that is consistent with the data and the known biological constraints (such as the relative positioning of the genes in the genome). The overall generation probability of a given sequence $s$ is obtained by summing of the probabilities of all scenarios that could have given rise to this sequence:

[TABLE]

The learning of the model parameters, encoded in the probability distributions in the factorized form (49), is performed by maximizing the likelihood of the sequences, using the Expectation-Maximization algorithm to deal with the sum over the hidden variable $r$ (see Sec. XII.5.1 for details about Maximum Likelihood and Expectation-Maximization).

This inference procedure has been applied to variety of immune receptors and species, TCR $\beta$ Murugan2012 and $\alpha$ Pogorelyy2017 chains, BCR heavy Elhanati2015 and light chains Toledano2018 in humans, as well as TCR$\beta$ in mouse Sethna2017 and BCR heavy chains in trout Magadan2018 . It is implemented in the IGoR software Marcou2018 which can be used to learn additional models from other species and locus combinations. IGoR can also be used to generate synthetic sequences from the $P_{\rm gen}$ distribution, and to estimate $P_{\rm gen}(s)$ of an arbitrary nucleotide sequence $s$. Another method, OLGA Sethna2019 , was designed to estimate the probability of amino-acid sequences, using dynamic programming to deal with the enormous sums involving the enumeration of all possible nucleotide variants. Other methods relying on hidden Markov models (see Sec. XII.5.1), were also proposed to handle the high dimensionality of hidden variables Munshaw2010 ; Elhanati2016 ; Ralph2016 ; Ralph2016a .

Repertoires generated in silico using the learned model confirm that the inference is able to call the correct recombination scenario in at best $\sim 25\%$ of cases, making the deterministic annotation of sequences impossible and necessitating a probabilistic inference approach. A perhaps counterintuitive consequence is that the model parameters can be inferred with arbitrary accuracy despite our inability to annotate any sequence reliably.

The results of the inference show that the probability of generation of receptors is incredibly reproducible between individuals of the same species. The distribution is dominated by insertions and deletions, whereas gene choice contributes relatively little to the overall probability. Among the features of V(D)J recombation, the gene choice distribution varies the most between unrelated individuals Murugan2012 , even when corrected for single nucleotide polymorphisms (SNP) in V, D and J alleles Wang2011 . As far as can be determined, the generations of the $\alpha$ and $\beta$ chains are largely independent, and the overall generation probability of a receptor is well approximated by the product of the generation of its two chains Dupic2018 . Although no similar analysis has been performed for BCR, it is likely that the formation of its two chains are also independent.

Although these recombination models ignore selection effects, to be discussed in the section below, it is a good predictor of the abundance of particular TCR$\beta$ in a human population Sethna2019 , and it has been used to detect signatures of immune responses as deviations from this baseline distribution Pogorelyy2018a ; Pogorelyy2018b .

VI.3 Thymic selection and central tolerance

After the receptors are generated, they undergo an initial selection step, known as thymic selection in T-cells. A similar process called central tolerance occurs in B-cells maturation Wardemann2003 , but let us focus on T-cell for concreteness. For a more comprehensive survey of models of thymic selection, we refer readers to the review by A. Yates Yates2014 , and focus here on recent development involving quorum sensing and data-driven models of sequence specific selection.

T-cells mature in the thymus, an organ that contains only proteins that are native to the host organism, termed self-proteins. The newly generated receptors are expressed on cell surfaces and their binding properties are tested against the peptides from the self-proteins, presented by the MHC complexes. If the receptor fails to bind any self-peptide, even weakly, it will probably fail to bind any protein and the cell carrying these receptors is discarded – a process called positive selection which removes $\sim 80\%$ of immature cells. Conversely, if a receptor binds any one self-protein too strongly, it is also discarded as a result of negative selection, since it is likely to bind self-proteins later on and trigger auto-immune diseases. Negative selection was at the core of the argument for the optimal recognition probability and repertoire size presented in Sec. VI.1, and it is estimated to remove $50\%-80\%$ of cells. The detailed process of T cell selection, complete with the timing of the $\alpha$ and $\beta$ chain recombination events, is summarized in Fig. 11. During positive selection, cells with TCR that recognize peptides presented by the MHC class I commit to becoming CD8+ cells, whose main function is to kill infected cells, while cells with MHC class II TCR specificity become CD4+ cells, whose main function is to help B cells in their affinity maturation, and to regulate the immune response. There are about 4 times as many CD4+ cells as CD8+ cells. A small subclass of CD4+ T cells, called regulatory T cells (Tregs), actually suppress the immune response, and play an essential role to prevent auto-immunity Wing2010 . Tregs are selected for higher affinities to self-epitopes than regular T cells, which allows them to selectively suppress immune responses to self-antigens, although the picture seem to be quite complex (see Bains2013b ).

To be able to ensure self tolerance, each receptor should in principle be tested against all possible self-epitope, which would take an impractically long time. In practice, experiments suggest that each T cell may encounter around 500 antigen presenting cells LeBorgne2009 , while theoretical arguments estimate the number of presented self peptides around a few thousands Butler2013a (these numbers are not inconsistent because each antigen-presenting cell may present many different peptides in single encounter with a T cell). These numbers are much lower than the estimated diversity of presentable self peptide-MHC complexes, $5\cdot 10^{5}-5\cdot 10^{6}$ (see Sec. VI.1). Calling $x$ the fraction of presented self peptides during negative selection, the probability that none of the selected TCR are self-reactive is $(1-p)^{RNx}$, which quickly goes to 0 as $R$ becomes even moderately large, even if $p$ scales as $1/N$.

A recently proposed solution to this problem is to assume that a minimal number $t$ of TCR must recognize a peptide to trigger an immune response Butler2013a . By virtue of the law of large numbers, the number of T cells responding to a self-peptide is distributed as a Gaussian with mean and variance $=Rp(1-x)$, while the number of T cells responding to a foreign peptide is of mean and variance $=Rp$. Good discrimination is achieved when the two distributions are well separated, giving the condition $Rpx^{2}\gg 1$, with an optimal discrimination threshold $t^{*}=Rp(1-x/2)$ in the limit $x\ll 1$. This condition is barely satisfied by $x=10^{-2}$, $R=10^{9}$ and $p=10^{-5}$, for which the optimal threshold for the number of TCR involved is $t^{*}\sim 10^{4}$. This argument suggests another, stricter lower bound on the diversity $R$ of the T cell repertoire that is necessary to make such a collective decision. It also implies a “quorum sensing” mechanism by which responding T cells have a way of estimating how many other cells are involved in order to commit to a response. Recent work suggests that such quorum sensing does occur locally, probably using cytokine signaling Polonsky2018 .

To model thymic selection in more quantitative detail, a common strategy has been to use additive models of binding free energy similar to (9) Detours2000 ; Kosmrlj2008 ; George2017 . In these models, the only receptors that survive are those whose maximal affinity $E^{*}$ to any self peptide falls within a range $(E_{p},E_{n})$ corresponding to the positive and negative selection thresholds. This condition can be mapped onto an extreme value statistics problem, allowing for a statistical mechanics treatment Kosmrlj2009 . Under this framework, it was shown that the sequence composition of negatively selected TCR was biased towards weakly binding residues Kosmrlj2008 , and the theory was subsequently applied to explain clinical data on “elite controllers” of HIV expressing a particular type of MHC class I molecule, HLA-B27 Kosmrlj2010 . A similar theory was used to study the sensitivity of TCR that target tumor cells George2017 . The parameters of these models are not inferred empirically, but are instead picked from popular but unrealistic interaction potentials between amino acids such as the Miyazawa-Jernigan matrix Miyazawa1996 . Yet many conclusions of these studies are relatively insensitive to the details of the interaction matrices in the extreme value statistics regime.

Statistical inference methods based on immune repertoire sequencing can also be used to estimate the probability that a particular receptor passes selection based on its sequence. We can define a selection factor corresponding to the ratio of probabilities to find a sequence in a selected repertoire, $P_{\rm sel}(s)$, relative to its probability in the unselected repertoire, as given by the recombination model $P_{\rm gen}(s)$: $Q(s)=P_{\rm sel}(s)/P_{\rm gen}(s)$. In practice, it is impossible to evaluate $P_{\rm sel}(s)$ directly, as the number of possible sequences to be considered is too large, spanning more than 20 orders of magnitude in generation probabilities. However, simplifying assumptions can be made on the form of $Q(s)$. Specifically, Elhanati et al Elhanati2014 considered:

[TABLE]

where $(a_{1},a_{2},\ldots,a_{L})$ is the amino-acid sequence of the Complementarity Determining Region 3 (CDR3, which forms a loop important for antigen recognition, and covers the most variable part of the receptor ranging from the end of the V to the beginning of the J segments). Selection is assumed to act independently on the length of the CDR3 region $q(L)$, the V gene, $q(V)$, and the J gene, $q(J)$. The $D$ gene is not taken into account seperately, but selection on amino acids in the CDR3 is considered explicitly. The parameters of the model are inferred by maximizing the likelihood using Expectation Maximization (Sec. XII.5.1). This approach has been applied to the initial selection process of $\beta$-chain sequences Elhanati2014 , $\alpha$-chain sequences Pogorelyy2017 , and BCR heavy chain Elhanati2015 and light chain Toledano2018 sequences. For many of these cases, the $Q(s)$ selection factors did not differ substantially from individual to individual. Instead, they reflected an global selection for general biophysical and biochemical features of amino acids (suggesting positive selection for proper protein function), rather than individual-specific removals in the repertoire caused by negative selection. Similar observations were made for BCR by direct comparison of pre-mature and mature repertoires Kaplinsky2014 . Modeling negative selection is much harder, as it requires to learn the “holes” that selection pokes into the repertoire, whose post-selection landscape looks like a dense golf course. A simple way to model negative selection is to remove a random fraction $1-q$ of the sequences from the selected repertoire, boosting the likelihood of surviving ones by a factor $1/q$ Elhanati2018a . In practice, in simulations one can use a hashing function, which associates a quasi-random number to each sequence, and select sequences whose hashing number falls below a chosen threshold, so that selection is reproducible, yet agnostic to the features of the sequences.

Note that models of the form (51) are not restricted to thymic selection, and can be used to describe the sequence-wide selection pressure during any process, by comparing the initial distribution before selection to the final distribution after selection, as was for instance done to characterize responsive receptors to yellow fever vaccination Pogorelyy2018b . For instance, different T cell subsets (e.g. CD4+ and CD8+ cells), different BCR isotypes, in different organs or environments, are probably characterized by distinct selective pressure, each of which could be modeled using similar $Q(s)$ factors.

VI.4 Diversity and the clone size distribution

The distributions of TCR and BCR span a very large space, and are highly skewed. It is often useful to summarize these distributions by a single summary statistics that quantifies their diversity. Diversity measures in the context of immune repertoires have been thoroughly discussed in a recent review Mora2016e . Diversity is important for understanding how well our immune repertoire prepares us against a wide range of pathogenic threats. As we will see, diversity measures are deeply linked with the distribution of sizes of immune clones (number of cells with the same receptor), which contains important information about the dynamics of immune cells in response to environmental challenges (see Sec. VII). In addition, it determines the amount of overlap one expects between repertoires of distinct individual, and underlies the concept of “public repertoire” shared by all individuals (see Sec. VI.5). All these aspects can be formalized within the common langage of statistical mechanics through the definition of a density of state, which we introduce below.

The most general family of diversity measures for a distribution $p(s)$ is given by Renyi entropies:

[TABLE]

This definition reduces to the Shannon entropy for $\beta=1$, $H_{1}=-\sum_{s}p(s)\ln p(s)$ (see Sec. XII.5.5), to the total number of receptors for $\beta=0$, $R=\exp(H_{0})$, and to the probability of drawing the same receptor twice with replacement for $\beta=2$, $\exp(-H_{2})=\sum_{s}p(s)^{2}$, also equal to the inverse of the the Simpson index. More generally, the family $H_{\beta}$ recapitulates the entire clone size distribution, defined in its cumulative form as the number $G(E)$ of receptor sequences with $-\ln p(s)<E$, through a Laplace transform:

[TABLE]

Note that $G(E)$ is formally equivalent to a cumulative density of states in statistical physics. $E=-\ln p$ is sometimes called the “surprise” and is also formally similar to an “energy” by analogy with Boltzmann’s law, $p\sim\exp(-E)$, making $\beta$ analogous to an inverse temperature. The mapping between $G(E)$ and $H_{\beta}$ is formally analogous to the equivalence of the canonical and micro-canonical ensembles, and is amenable to a statistical physics analysis Mora2016d . The lower the temperature (the higher the $\beta$), the more the Rényi entropy concentrates on low-energy, high probability sequences.

$G(E)$ can also be interpreted as the rank (ordered from most probable to least probable) of sequences with probability $p=e^{-E}$. The distribution of clone sizes is often presented in terms of a clone-size frequency rank distribution, where clones are ranked ordered by their sizes, and their normalized frequency is plotted as a decreasing function of its rank. $G(E)$ precisely encodes this rank-frequency relation.

It is important to distinguish between the potential diversity, corresponding to a theoretical distribution $P_{\rm gen}(s)$ or $P_{\rm sel}(s)$, and the realized diversity, corresponding to the actual distribution of a finite number of receptors in a particular sample, $p(s)=n(s)/N_{\rm tot}$, where $n(s)$ is the number of molecules or cells with receptor sequence $s$, and $N_{\rm tot}=\sum_{s}p(s)$.

Let us start with the potential diversity derived from the probability distribution $P_{\rm gen}(s)$ (see Sec. VI.2). For human TCR $\beta$, the total number of possible sequences, $e^{H_{0}}$, is infinite for all practical purposes ($>10^{39}$ Mora2016e ), as it is larger than the total number of TCR$\beta$ receptors having ever been produced by a human being. For this reason, it makes more sense to compare Shannon entropies. These can be computed using (53), where $G(E)$ is evaluated from the probability density distribution of $E$ among randomly drawn sequences according to $p(s)$, $\rho(E)=(dG/dE)e^{-E}$, so that $G(E)=\int_{0}^{E}dE^{\prime}\,\rho(E^{\prime})\exp(-E^{\prime})$.

Entropy estimates for both nucleotide and amino acid sequences are reported in Sethna2019 . Human nucletoide TCR $\beta$ diversity is $H_{1}\sim 44$ bits (bits refer to $\ln(2)$ units) for nucletoides and $30$ bits for amino-acids, TCR $\alpha$ divesity is $\sim 30$ bits for nucleotides and $25$ bits for amino acids. The BCR heavy chain entropy is $67$ bits for nucleotides and $53$ for amino acids. These numbers are very large. If the distribution were uniform, that would correspond to $2^{55}\approx 3.6\cdot 10^{16}$ distinct amino acid TCR$\alpha\beta$ sequences. Because the distribution is non-uniform, the pooled repertoire of all humans having ever lived do not exhaust the potential diversity. As noted before, most this diversity is due to random insertions and deletions.

Thymic selection and central tolerance effectively reduce these potential diversities in $P_{\rm sel}(s)$, by around $9$ bits for TCR$\beta$ Elhanati2014 , $4$ bits for TCR$\alpha$ Pogorelyy2017 , and $12$ bits for BCR heavy chain Elhanati2015 , but the corresponding diversity number remain very high. This diversity loss does not mean that e.g. $2^{-13}\approx 10^{-4}$ of TCR must be discarded; the probability of thymic selection $Q(s)$ is in fact correlated with generation probability $P_{\rm gen}(s)$, meaning that diversity is reduced by selecting already likely sequences.

Let us now turn to the realized diversity in a sample or an individual. This diversity is best described by the whole clone-size distribution encoded by $G(E)$. Based on recent high-throughput sequencing experiments, the distributions are often long-tailed, spanning several orders of magnitude (see Fig. 12) Mora2010 ; Desponds2016 ; Desponds2017 ; Mora2016e . These tails seem to be mostly due to the memory fraction of the repertoire Greef2019 . Initial worries that these long tails may arise from noise is the Polymerase Chain Reaction, which amplifies differences exponentially Best2015b , were put to rest thanks to the introduction of unique molecular barcodes associated to each initial mRNA molecule Vollmers2013 ; Shugay2014a ; Heather2017 . Specifically, the clone size distribution follows a power-law:

[TABLE]

with $\beta\leq 1$ but close to 1. The probability distribution of clone frequencies,

[TABLE]

is also a power law of exponent slightly below 2. Another way to look at the distribution is to consider the inverse relationship,

[TABLE]

which also follows a power law, called Zipf’s law in this particular context when its exponent $\hat{\beta}^{-1}$ is close to 1. Since $G$ is the rank of the clone (from most frequent to least frequent), this relation is called the rank-frequency relationship.

Because of these long tails and limited sampling, very few of the diversity measures can be estimated reliably. One way to overcome these issues is to build a statistical model that is then parametrised by the experimental data. In Sec. VI.1 we already mentioned Fisher’s Poisson abundance model and the Chao estimators, which use this idea. More sophisticated methods have been proposed for immune repertoire, e.g. using Expectation Maximization (see Sec. XII.5.3) Kaplinsky2015 or using rarefaction curves Laydon2015 , but all these methods are susceptible to huge errors when the clone size distribution follows a power law Haegeman2013 . We expect the reported numbers to grossly underestimate the true diversity of both BCR and TCR. To make progress, models of lymphocyte population dynamics should be leveraged to reliably extrapolate the tails of clone-size distributions. We will discuss some of these models in Sec. VII.

VI.5 Repertoire sharing

Receptor sharing is at the heart of the public vs private debate in immunology. Public receptors are those shared between individuals in a given cohort, while private are ones that are seen only in individuals. If one studies a group of individuals with a specific condition (e.g. cytomegalovirus or CMV Emerson2017 or Ankylosing Spondylitis Faham2017 ), receptors that are shared between these people can be expected to be linked to a disease. This justifies the search and investigation of “public” sequences. Both humans and mice do share a non-negligible numbers of receptors in functional repertoires, regardless of their lifestyles, environmental factors and family ties. However, receptors are also expected to be shared by chance, because some sequences are likely to recombine independently in different individuals Venturi:2006hk ; Venturi:2011co ; Elhanati2014 ; Madi2014a . One can try to correct for this convergent recombination by using an outgroup — a group of individuals that are as similar to the condition-specific group but do not have this condition — and identify sequences shared more in the studied group than in the control group, as was done for CMV Emerson2017 . Alternatively, one can use the knowledge of immune repertoire diversity encoded in $P_{\rm gen}(s)$ to estimate our null expectation.

All sharing properties are recapitulated by the cumulative density of states $G(E)$. Given $n$ independent repertoires of sizes $N_{1},\ldots,N_{n}$, the expected number of sequences shared by exactly $m$ repertoires, $M_{m}$, is given by the following generating function Elhanati2018a :

[TABLE]

For example, for $n=2$, in the limit of small samples, $N_{1,2}e^{-E}\ll 1$, the expected number of shared sequences is

[TABLE]

reducing to simple biased birthday problem related to Simpson’s diversity and thus Rényi entropy of order 2. More generally, sharing between $m$ small samples is governed by $\int dG(E)\,e^{-mE}$ and to Rényi entropies of order $m$. The ensemble of sequences shared between $m$ samples thus corresponds to sampling the generation probability at temperature $1/m$, focusing on more and more likely sequences as $m$ increases.

This calculation shows that in such large ensembles of cells as the B or T cell repertoire, which contain upward of $10^{8}$ unique receptors, we expect the most common receptors with a high probability to be independently generated multiple times in different people. The number of shared TCR$\beta$ sequences among sampled repertoires from a cohort of 658 human donors Emerson2017 was well predicted by the model in Eq. VI.5. By itself, the generation probability $P_{\rm gen}$ was enough to build a classifier that can determine whether a given sequence will be shared by a minimal number of people, with a $>95\%$ accuracy for $m>10$. This analysis also brings our attention to the definition of public repertoires. A receptor that is shared between $m=2$ and $m=1000$ people has a very different notion of publicness. For this reason it is better to talk about degrees of publicness. In addition, sharing and publicness depends very sensitively on the collected sample size $N_{i}$. A reasonable choice for a single definition of receptor publicness could that the sequence must be present on average once in each individual, implying $P_{\rm gen}(s)>1/R$.

Interestingly, the theory can also identify receptors that are shared for other reasons than sheer chance. These receptors have a high propensity to be shared, but a relatively low generation probability. The discrepancy between chance prediction and observed sharing was used to discover shared T cell clones in identical twins Pogorelyy2017 , to find candidate TCR sequences associated with conditions such as CMV or diabetes without a control group Pogorelyy2018b , or to assess the public BCR response to vaccination in trout Magadan2018 .

Identifying common, even unusual receptors does not necessarily mean they are responding to the same antigen or pathogen, and even when it does it does not tell us what the target antigen is. As discussed in Sec. II.4, immune response is a complex phenomenon that involves the binding different epitopes to immune receptors, followed by signal propagation, cell commitment, and cell-cell communication through messenger molecules. Yet tools based on the statistical expectations of an unbiased repertoire allow us to identify deviations from this baseline, and propose candidate responding receptor sequences to be tested in other experiments.

VI.6 Optimal receptor distribution

In Sec. VI.1 we discussed theoretical arguments for an optimal repertoire size, but recent data suggests that actual repertoire sizes are much bigger than these theoretical predictions. However, as stressed in Sec. VI.4, the distribution of receptors is itself highly skewed, suggesting that the relative abundance of receptor types, rather than their absolute number, is an important factor of repertoire design. Can we find an optimization principle over the composition of the repertoire? Take a distribution $Q_{a}$ of antigens to be recognized, where $a$ lives in a theoretical high-dimensional “shape space” encompassing both antigens and receptors, and in which proximity reflects receptor-ligand affinity Perelson1979 ; Perelson1997 , through a recognition or “cross-reactivity” matrix $f_{ar}$ between antigen $a$ and receptor $r$ which encodes the probability that a random encounter between the two results in an immune response. For a given distribution of receptors, $P_{r}$. The expected cost of an infection linked to $a$ is a decreasing function $c(\cdot)$ of the probability that a chance encounter with a receptor is successful, $\sum_{r}f_{ar}P_{r}$, so that the expected cost for a random infection reads Mayer2015a :

[TABLE]

Minimizing that cost with respect to the receptor distribution $\{P_{r}\}$ gives an interesting multi-peak structure (Fig. 13B). Because of the degeneracy of recognition allowed by cross-reactivity, similarly well protecting repertoires in different individuals may exist, even if they sample essentially the same pathogenic environment (Figure 13B). As a result, we should not be surprised by great differences at the phenotypic level, let alone at the sequence level, when looking at repertoires responding to specific threats. In other words, even if convergent selection is at play, each immune repertoire may find a different molecular solution to it, which may explain the observation that most repertoire sharing seem to occur by chance.

VI.7 Repertoire response to an immune challenge

RepSeq technologies allow one to track changes in the repertoire in time and in response to environmental changes, infections, or vaccinations. In many animal models such as mice or fish, it is generally not possible to sample the repertoire in the same individuals. Instead, the repertoires of distinct but isogenic individuals are sampled at the different times before and after an immune challenge. By construction, this strategy restricts the analysis to public features of the repertoire response. To identify common sequences features of repertoires that responded to a particular challenge, machine learning techniques based on low-level features Thomas2014b or support vector machines Cinelli2017 have been used in mice to classify responding repertoires. Enrichment of particular V classes and sequences features have also been reported in the BCR repertoires of trout following immunization Castro2013 ; Magadan2018 .

In humans, the BCR response to influenza vaccination Vollmers2013 ; Laserson2014 , BCR and TCR response to varicella-zoster vaccination Wang2015 ; Qi2016 , and the TCR response to yellow fever vaccination (a model for an acute infection) DeWitt2015 ; Pogorelyy2018 have been studied at the repertoire level using RepSeq. These studies are mostly descriptive, but they pave the way towards a better characterization of the specificity, reproducibility, and publicness of the immune response at the repertoire level. A common observation is that the response is mostly private, although responding clones are more shared than would expected by chance. Responding sequences also tend to be clustered in sequence space Pogorelyy2018 . This clustering gives a criterion for identifying responding sequences in a single repertoire snapshot, based on the density of similar sequences (differing by at most one amino acid) in the repertoire, relative to this density expected by the recombination model, $\sum_{s^{\prime},|s-s^{\prime}|\leq 1}P_{\rm gen}(s^{\prime})$ Pogorelyy2018b .

VII Lymphocyte population dynamics

T and B cells are organized in clones of cells that express the same immune receptor. These clones grow and decay as individual cells divide or die either spontaneously or in response to external signals — e.g. an antigen recognition event triggering cell proliferation, or cell-cell communication through messenger molecules called cytokines. The correct way to approach and to model these complex lymphocyte dynamics, as well as their implication for experimental observables such as clone size distributions or clone expansion, is still a largely open question. There is a vast body of literature on modeling lymphocyte dynamics for broad subpopulations with no clonal information, using models ordinary differential equations. Ref. Perelson2002 provides a useful entry point into that literature. Here we will focus on stochastic models population dynamics, where the identities and sizes of individual lymphocyte clones are tracked. This approach is in line with recent developments in high-throughput repertoire sequencing (Sec. VI), which allows for such a fine-grain description of population dynamics, and for which much work remains to be done, both experimentally and theoretically. While progress on this topic is still in its infancy, it promises to give us better insight into the collective decision-making and dynamics of populations of immune cells, and to extract important in vivo parameters from data, at a moment when most studies are based on cultured cells or mouse models.

We also briefly review the application of differential equation models to HIV dynamics, in what constitutes perhaps the most successful application of modeling to translational immunology.

VII.1 Neutral dynamics

Let us fist consider a simple stochastic model of repertoire evolution in the absence of any antigenic proliferation. Such a model is relevant for the naive pool, which is often believed to be unaffected by external stimuli (although this is also debated, see below), but can also serve as a “neutral” baseline against which to compare other dynamics or longitudinal data. While some of the results discussed in this section were presented in Desponds2017 , here we expand on them and add a few original results, notably on the establishment of the steady state, and on estimates of the total number of clones.

In this simple model, new clones of size $C$ come out of the thymus (for T cells) or bone marrow (for B cells) with rate $\theta_{C}$, so that the total cell output $\theta=\sum_{C}C\theta_{C}$. Then each cell can divide with rate $\nu$, and die with rate $\mu>\nu$.

The total number of cells $C_{\rm tot}(t)$ follows the simple differential equation:

[TABLE]

This simple equation is at the basis of many studies aiming at quantifying lymphocyte population dynamics, using experiments with isotope labeling Borghans2007 ; DeBoer2013 or other markers of cell divisions Bains2009a . When thymic output, division and death rates are constant, Eq. 60 is solved by (assuming no cell at $t=0$):

[TABLE]

At steady state, the total number of lymphocytes $C_{\rm tot}=\theta/(\mu-\nu)$ reflects the balance between thymic output, cell division and death. It is believed that the division and death rates ajust themselves so as to keep $C_{\rm tot}$ constant, a process called homeostasis. For T cells, $\theta\sim 10^{6}$/day for mice, and $\theta\sim 10^{8}$/day for humans, although that number varies with age. Taking $C_{\rm tot}\sim 10^{11}$, this gives an effective decay time of $(\mu-\nu)^{-1}\sim 1,000$ days, which is consistent with lifetime estimates of T cells in human using deuterium water DeBoer2013 .

However, these numbers are not informative about the repertoire structure and its clone size distribution. These calculations are also unable to disentangle division and death, lumped in the single parameter $\mu-\nu$. The expected number of clones of size $C$, $N_{C}(t)$, is governed by the following equations Desponds2016 :

[TABLE]

At steady state, $dN_{C}/dt=0$, the solution for $C>\max\{C:\theta_{C}>0\}$ is a power-law of exponent 1 with an exponential cutoff:

[TABLE]

When birth and death are balanced, $\nu\sim\mu$, the exponential cutoff disappears. Comparing to the power law Eq. 55 in the distribution of frequencies $\rho(p)$ with the correspondance $N_{C}=N_{\rm tot}\rho(p=C/C_{\rm tot})$, this model would predict a exponent $\hat{\beta}=0$, which is not supported by repertoire data on unsorted or memory cells. This suggests that the dynamics of the memory repertoire cannot be simply neutral as modeled by Eq. 62, but must be affected by selection, such as expansions events following antigen recognition, as we will see in the next section.

The dynamics of Eq. 62 can in fact be solved analytically using generating functions and the method of characteristics Anandprivate . Assuming for simplicity that clones have initial size 1, $\theta_{C}=\theta\delta_{C,1}$, we obtain (see Sec. XII.4.3):

[TABLE]

This relation is still a power law with an exponential cut-off. The power law exponent remains 1 at all times, but the cut-off gets larger and larger as the system reaches steady state. The total number of clones, $N_{\rm tot}(t)$, then evolves according to:

[TABLE]

and the total number of cells, $C_{\rm tot}(t)=\sum_{C}CN_{C}(t)$ follows Eq. 61. At steady-state, the total number of clones, $N_{\rm tot}=(\theta/\nu)\ln[\mu/(\mu-\nu)]$ depends on both the division and death rates, and not just their difference. In the limit of small division rate, $\nu\ll\mu$, which is often assumed for naive cells, we get $N_{\rm tot}=\theta/\mu\approx C_{\rm tot}$. Conversely, in the limit $\mu-\nu\ll\mu$, one has $N_{\rm tot}\sim C_{\rm tot}\epsilon\ln(1/\epsilon)$ with $\epsilon=(\mu-\nu)/\mu$, meaning that the number of clones may be arbitrarily smaller than the number of cells. Starting with initial clone sizes of $k_{0}$ does not substantially affect this picture, with $N_{\rm tot}\sim(H_{k_{0}}/k_{0})C_{\rm tot}$ (where $H_{n}=\sum_{i=1}^{n}1/i$) in the limit $\nu\ll\mu$, and with unchanged scaling $N_{\rm tot}\sim\epsilon\ln(1/\epsilon)C_{\rm tot}$ in the limit $\mu-\nu\ll\mu$ (see Sec. XII.4.3).

These models can be refined by accounting for convergent recombination, by which clones sizes can also be increased by thymic exports that have the exact same sequence, which happens with probability $P_{\rm gen}(s)$ for a particular sequence $s$. With this correction, Zheng and collaborators Zheng2017 found fair agreement between the results of a neutral model with a source Eq. 62 and the bulk distribution of clone frequencies observed in the naive repertoire of mice. However, they also found many outliers, i.e. large naive clones that cannot be explained by the neutral assumption. A major difficulty of such a comparison is that, as we already mentioned in Sec. VI.4, the clone frequency distribution of the full repertoire is heavy-tailed. Any small contamination of non-naive cells into the studied repertoire is likely to introduce spurious outliers in the distribution, confounding the analysis.

VII.2 A note about “neutral processes”

In this review we discuss four different types of processes that are referred to as neutral in different domains: (i) the neutral model of population genetics (see section XII.3.8), (ii) the neutral clone size distribution model (see section VII.1), (iii) the Yule-Simon process of speciation (see section XII.4.5) and (iv) the neutral ecological model (see section XII.4.1). Since these models are often confused with one another, sometimes simply for semantic reasons, otherwise because many of them give power law distributions for some observable, we thought it is pedagogical to discuss how they are different and summarize the exponents they give and what kind of distributions they predict.

Many of these models predict a power law in the distribution of clone or species frequencies, but details differ. Table 1 gives a summary of these distributions and differences between the model assumptions.

Neutral models of ecology come in different flavours Hubbell2001 , but the simplest one is essentially equivalent to Eq. 62, and the corresponding clone size distribution, called the species abundance distribution in that context, has long been known as Fisher’s logseries Fisher1943 , althouth Fisher derived it through purely statistical means using unjustified assumptions.

The neutral model of population genetics is similar to the neutral model of Eq. 63, with a major difference: the population size $C_{\rm tot}$ is fixed. Most often, novel types arise from mutations arising in existing types. In population genetics the “site-frequency spectrum” corresponds to the distribution of allele frequencies, and is thus the equivalent of the clone size distribution $N_{C}$. In the neutral model, the site-frequency spectrum is a straight power law of exponent $-1$. The difference with the result of Eq. 63, beyond the absence of an exponential cutoff, is that mutants that reach the (fixed) population size $C_{\rm tot}$ become the wildtype, and are thus removed from the dynamics of mutants, while in lymphocyte dynamics no single clone ever takes over the whole population.

The Yule process is a model of an ever expanding population with mutations. It assumes that individuals of all types divide with rate $\mu$, with a mutation probability $\alpha$ that generates new types. The crucial difference with Eq. 63, beyond the absence of death, is that the rate of novely increases linearly with the population size. Also, the population expands exponentially, and the clone size distribution only reaches a quasi steady state, after renormalizing the total population size, $C_{\rm tot}$. The renormalized size distribution asymptotically follows a power law of exponent $-2$ in the limit $\alpha\ll 1$. Although this exponent is tantalizing close to those observed in data, the assumption that new types arise with a rate proportional to the population size is incompatible with the biology of T cells, as new clones originate from the thymus, and not from already circulating cells. On the other hand, the Yule process may be appropriate for modeling population-wide affinity maturation of B cells, as their receptors hypermutate upon expansion.

VII.3 Population dynamics model with external signals

To explain the power laws observed in the clone size distribution from sequenced repertoires, one needs to go beyond neutral models and introduce external signals affecting the division and death rates — or fitness — of cells as a function of their phenotypic state or of their immune receptor. Several models of clonal dynamics have been proposed to describe the growth and decline of clones in various populations (see section XII.4.1 for a general discussion of competition models). We will not describe all of these models, but try and give the reader an idea of their general features. The key ingredient is that cells divide in response to antigenic stimulation, which depends on the concentration or frequency of available antigens that are susceptible to be recognized. Cells in a clone expressing receptor $r$ die with constant rate $\mu$ as in Eq. 62, but they divide with clone-specific rate

[TABLE]

where $f_{r,a}$ is the cross-reactivity function between an antigen and a receptor, $Q_{a}$ the frequency of antigen $a$, and $A(x)$ is a decreasing function describing the availability of a given antigen as it is being bound by other receptors than $r$; for example, $A(x)=(1+\epsilon)/(1+\epsilon x)^{2}$ with $\epsilon$ setting the strength of competition. In this model, receptors interact only indirectly through this competition factor.

As a historical note, models with direct interactions between receptors, namely antibodies binding to each other, were once proposed as a way to generate interesting dynamics on so-called idiotypic networks Jerne1974 . While idiotypic networks were often used to explain phenomena in the 1980’s, and became very popular with physicists because of their link to network theory and the physics of disordered systems Perelson1997 , it is not currently considered a dominant paradigm in immunology (and already was not by the time Ref. Perelson1997 was written, as discussed therein), as its predictions can be explained by clonal selection theory alone without invoking additional elements.

In the linear noise approximation, the size of the clone of cells expressing receptor $r$ is described by a stochastic differential equation:

[TABLE]

with $\langle{\xi_{i}(t)\xi_{i}(t^{\prime})}\rangle{}=\left(\nu_{r}\left(\{Q_{a}\},\{C_{r^{\prime}}\}\right)+\mu\right)C_{r}\delta(t-t^{\prime})$ with Itô’s convention. This general form of the model can describe a wide variety of situations. It can be used as a model for the naive T cell repertoire, where $Q_{a}$ is the distribution of self-antigens for which T cells compete for survival cues. Such models of competitive exclusion DeBoer1994 have been studied and used to explain biological facts for both T- DeBoer1997 and B cells DeBoer2001 repertoires, such as the fact that a larger diversity is beneficial in terms of overall repertoire fitness. Interestingly, the same dynamical equation emerges as a way to approach an optimal repertoire Mayer2015a , as defined in Sec. VI.6.

Eq. 67 can also be used to model lymphocyte population dynamics in response to transient antigenic stimulation. Antigens appear with some rate $s$, and decay as they are being cleared by the immune system as $Q_{a}(t)=Q^{0}_{a}e^{-\lambda(t-t_{a})}$, where $t_{a}$ is the time at which antigen $a$ appeared in the population. In this context, the antigen is not limiting, and is fully available, $A(x)=1$. When $\lambda$ is very large compared to $\mu$ and $\nu$, each clone experiences a series of short spikes of expansion of magnitude $\exp(\nu\langle Q^{0}\rangle/\lambda)$ with rate $sp$, where $p$ is the probability that $f_{ar}=1$ for a random choice of $a$ and $r$ (and $f_{ar}=0$ with probability $1-p$, in an all-or-nothing approximation). With the change of variable $x=\ln C$, each expansion event causes an average jump $\Delta x=\nu\langle Q^{0}\rangle/\lambda$.

In Ref. Desponds2016 , expansion events were assumed to be both small and frequent compared to other time scales, which led to an effective diffusion equation for $x$. Here, as an original result of this review, we present the solution to the general process. A more detailed solution using Laplace transforms is presented in Sec. XII.4.4). The density of clones with logarithmic size $x$, $\rho(x,t)=(dC/dx)N_{C=e^{x}}(t)=e^{x}N_{e^{x}}(t)$, follows the simple jump process with constant negative drift:

[TABLE]

where $\tilde{\theta}(x)=e^{x}\theta(C=e^{x})$ is the density rate of new clones entering the population with initial logarithmic size $x$. Assuming a power law Ansatz for the clone size distribution,

[TABLE]

translates into exponential decay in $x$, $\rho(x)=\rho_{0}e^{-\alpha x}$, giving the consistency equation for $\alpha$:

[TABLE]

The total numbers of cells and clones can be calculated at steady state (see Eqs. 204,205 in Sec. XII.4.4), and read in the case of fixed introduction size $\theta_{C}=\theta\delta_{C,C_{0}}$:

[TABLE]

If we assume many small expansion events, $\Delta x\ll 1$, expanding the exponential at second order allows us to recover te power-law exponent of Desponds2016 :

[TABLE]

where the numerator is the net decay rate of the clone size obtained from the difference between death and expansion, while the denominator corresponds to an effective diffusion coefficient stemming from the random arrival of expansion events.

The result of Eq. 69 suggests that a wide range of models with clone-dependent, random expansion or selection events produce power laws with arbitray exponents, in agrement with the data on memory or unfractioned lymphocyte repertoires. To distinguish between different types of dynamics, more details of the dynamics of individual infection events, as well as of the kinetics of unstimulated clones should be studied in terms of their impact on the clone size distribution.

The hypothesis that the expansion factor $\Delta x$ is constant, or at least does not depend on the clone size, is questionable. For instance, in secondary immune responses, the expanded cells already have a memory phenotype, and these cells tend to expand less than naive cells upon a primary infection. In addition, the amount of inflammatory signals that promote expansion should in general depend on the number of cells involved in the response. From a design perspective, it would seem beneficial not to expand clones that are already big, as the organism is already well protected against the pathogens that these clones are specific to. The optimal expansion strategy upon each pathogen encounter can be calculated within the simplified framework described in Sec. VI.6. Using Eq. 59 with a logarithmic cost $c(x)=-\ln(x)$ and a uniquely specific Kernel $f_{ar}=\delta_{a,r}$, combined with Bayesian prediction theory, the optimal dynamics can be shown to approximately follow Mayer2018 :

[TABLE]

where $\tau_{m}=2\tau/(\chi^{-1}C_{\rm tot}-1)$ (with $\tau$ the effective timescale of the pathogen dynamics), $E_{r}$ is the set of infection events $i$ occuring at time $t_{i}$ in response to which lymphocytes carrying receptor $r$ expand, and $\chi$ is a scaling parameter setting the total number of lympthocytes. The main difference between these and previously considered dynamics is that the growth of $C_{r}$ is not exponential, meaning that regulation or homeostasis mechanisms must exist to tune the magnitude of the expansion as a function of the size of the antigen-specific repertoire subset.

So far we have assumed that selection manifests itself at the level of clones: all cells in the same clone has the same average growth rate. However, selection could be cell dependent, which is the case if cells respond differently to non-antigenic stimulation signals, such as cytokines. Growing experimental evidence shows the heterogenous response to cytokines, due to differential expression of cytokine receptors, signalling molecules and their diffusion. Models with a cell-dependent fitness do not give power law behaviour, but do produce long tails Desponds2016 . With the current experimental cut-offs, it may not be possible to discriminate between clone-specific and cell-specific fitness models using repertoire data only. Further studies with longitudinal tracking of clone sizes may help settle these questions.

VII.4 In host HIV dynamics.

Perelson and collaborators Perelson2002 ; PerelsonRibeiro considered a modified SIS model of an HIV infection. This is one of the most influential examples of how computational models influenced medicine. It was used to predict the effects of anti-viral drugs on HIV: reverse transcriptase (RT) which blocks the ability of HIV to infect a cell, and protease inhibitors (PI) that result in the production of non-infectious viruses. Target cells $T$, that correspond mainly to CD4+ T-cells, become infected ($I$) at rate $\beta$, are born with rate $\alpha$ and die with rate $\mu$. Infected cells die at rate $\nu$, but they also help the virus ($V$) reproduce and produce $p$ new virions that are cleared with rate $c$.

[TABLE]

To study the effect of an antiviral drug, the virus equations are modified to account for two viral species: non-infectious viruses ($V_{\rm NI}$) and infectious viruses ($V_{\rm I}$):

[TABLE]

where $\beta=1-\epsilon_{\rm RT}$ and $0<\epsilon_{\rm PI}<1$ and $0<\epsilon_{\rm RT}<1$ are efficacies of RT and PI drugs. The total viral populations is held constant $V_{\rm NI}+V_{\rm I}=V=\rm{const}$. For $\epsilon_{\rm PI}=1$ and $\epsilon_{\rm RT}=1$, assuming $T={\rm const}$, the total viral load will decay according to PerelsonRibeiro ; Perelson2002 :

[TABLE]

This decay was fit to data from an HIV-infected patient on anti-viral therapy by adjusting the $c$ and $\nu$ parameters using least-squares regression, showing very good agreement Perelson1996 . Explaining viral dynamics in the case of combination therapy (tri-cocktails) requires introducing long-lived infected cells that act as a secondary viral source Perelson1997a . Fitting these more complex, yet still extremely simple models to data, allowed Perelson and collaborators to estimate the half-lives of the different species in Eq. VII.4. This in turn made it possible to effectively administer a tri-therapy cocktail treatment in intervals that made it very hard for the virus to escape. This extremely simple calculation saved lives.

VIII Affinity maturation

Upon antigenic stimulation, B-cell receptors acquire somatic hypermutations (SHM) that can help them explore the binding landscape. Based on high-throughput repertoire sequencing data, combinations of clustering and tree building methods Mccoy2015 ; Yaari2015 have been proposed to characterize lineages. Since affinity maturation is a fascinating example of a Darwinian evolutionary process, a lot of effort is going into extracting the details of the evolutionary processes from data, along with more theoretical efforts. We very briefly review tree building, clustering and lineage reconstruction approaches that are available in different software in section XII.3.11. A recent review more thoroughly sumarizes these attempts Hoehn2016 , arguing for methods departing from the traditional assumptions of population genetics, and tailor-made for the hypersomatic mutation scenario. BCR repertoire sequencing provides large amounts of data that allow software packages like Immcantation Gupta2015 ; VanderHeiden2014 ; Yaari2015 ; Yaari2015a ; Cui2016 ; Yaari2013 , Partis Ralph2016 (see section XII.3.11) SPURF Dhar2018 or IGoR Marcou2018 ; Mccoy2015 to learn hypermutaiton models (as described below), and understand the evolutionary process. As mentioned in section XII.3.11 the problem of lineage reconstruction is quite general to many areas in immunology. Incorporating the observation that more abundant clones are likely to have more offspring allows for reconstructing B-cell lineage ancestries from germinal center imaging experiments DeWitt2018d .

Affinity maturation is a sped-up Darwinian evolution process. Upon recognizing an antigen even with low affinity Jacob1991 , naive B-cells migrate to germinal centers (GC) in lymph nodes and about 3 days after antigen injection the germinal center reaction starts Nieuwenhuis1984 , initially dividing without hypermutations (see section V.4.1) for about 7 days and reaching a population size of $\sim 1500$ cells Jacob1991 . About a week after the infection starts, they first enter the dark zone region of the germinal center, where the Activation-induced cytidine deaminase (AID) enzyme introduces hypermutations as the cells continue to divide.

VIII.1 Modeling the hypermutation process

Repertoire sequencing analysis of unselected mutations Yaari2015 ; Elhanati2015 ; Marcou2018 has shown that AID acts highly non-uniformly, with hypermutation hotspots and strong context dependent motifs. The difficulty in learning these models stems from eliminating selection that acted on memory sequences. A 5-mer model, called S5F, was learned from synonymous mutations (that are in principle free of selection) from long heavy chain sequences Yaari2015 . Only half of possible 5-mers were observed from data with synonymous mutations, the remaining ones were inferred by averaging over related observed pentamers. Alternatively, out-of-frame memory sequences were used to learn heavy V and J segment di- and tri-nucleotide context dependent hypermutation models from low throughput sequencing experiments shapiro1999predicting .

The S5F model learns the probability $P_{\rm mut}(\vec{\pi})$ of mutating the central position of a 5-mer motif $\vec{\pi}=(\pi_{-2},\ldots,\pi_{2})$ by counting occurences of synonymous mutations in this motif centered at position $i$ along the sequences $\vec{\sigma}^{j}$ of a large dataset. Because only synonymous mutations must be taken into account, ad hoc heuristic normalization rules were used. In practice, mutation rates are given by:

[TABLE]

with $\mu$ an adjustable baseline mutation rate (not learned within S5F), and

[TABLE]

is the expected number of mutations of each 5-mer in a sequence. This number is obtained by normalizing for being synonymous in each sequence, and then averaging over sequences weighted by their number of synonymous mutations. In that expression, $C^{j}_{\vec{\pi}}$ is the observed number of synonymous mutations in sequence $\vec{\sigma}^{j}$ in motif $\vec{\pi}$, and $B^{j}_{\vec{\pi}}$ is a background normalization:

[TABLE]

where: $\vec{\tau}^{j}$ is the germline (unmutated) ancestor of $\vec{\sigma}^{j}$; $I_{\vec{\tau}^{j}}(i,\vec{\pi},b)$ is an indicator function that is equal to 1 if $\vec{\tau}^{j}$ matches $\pi$ over the position range $(i-2,\ldots,i+2)$, and $\tau^{j}_{i}\to b$ is a synonymous mutation; finally $s^{\vec{\pi}}_{b}=N^{\vec{\pi}}_{b}/\sum_{b^{\prime}\neq\pi_{0}}N^{\vec{\pi}}_{b}$ and $s^{\vec{\pi}}_{\pi_{0}}=0$, where $N^{\vec{\pi}}_{b}$ is the total number of synonymous mutations to $b$ observed in 5-mer motifs $\vec{\pi}$ across all sequences in the training dataset.

The S5F model provide a profile of hypermutation hot and cold spots in the absence of selection, which is widely used for analyzing BCR data. In particular, it can used as a baseline for quantifying selection from synonymous mutations Uduman01072011 ; Yaari2013a . Note that a 5-mer model was also learned from VH and JH genes sequences in rearrangements engineered not to be productive, in both the heavy and light chains of mouse, without having to use the trick of using synonymous mutations Cui2016 .

To avoid for over-fitting, one can further assume that the impact of the motif on mutability is additive. This strategy was applied to non-productive heavy-chain sequences in humans Marcou2018 . Mutability of $(2m+1)$-mer motifs $\vec{\pi}$ is given by:

[TABLE]

where $(\pi_{-m},...,\pi_{m})$ is the $2m+1$-mer sequence context around the mutation $i$, and $e(\pi_{i})$ are the inferred elements of the Position-Weight Matrix (PWM), learned within the IGoR framework Marcou2018 (see Sec. VI.2) using the Expectation-Maximization algorithm (see Sec. XII.5.1).

The actual hypermutation process is yet more complicated than captured by all these models. It has been pointed out that because of the context dependence of mutations, the order in which they appear matters. Since considering all possible histories is computationally not feasible, Feng et al proposed to Gibbs sample the orderings to learn a hypermutation model that was validated on simulated data Feng2017 . More importantly, hypermutations are not independent of each other, and tend to cluster along the sequence Marcou2018 . Hypermutations result from a lesion of DNA followed by error-prone DNA repair over extended regions ($\sim$ 20 nucleotide) of the sequence Unniraman2007 . These events may induce several simultaneous mutations at close-by positions along the sequence. It also suggests that context can influence the mutation rate over fairly large distances, since the erroneous repair can occur far from the original damage causing the lesion. Finally, context models cannot explain all of the variance of hypermutation rates, suggesting that other factors, such as length dependence, may be at play.

VIII.2 Cycles of selection

After acquiring hypermutations B-cells pass to the light zone of the germinal center where they undergo the selection step of Darwinian evolution. In this step a B-cell receptor must bind the antigen presented on a follicular dendritic cell (FDC) strongly enough for the B-cell to internalize the antigen and present them to helper T-cells. B-cells compete for binding of their pMHC to these helper T-cells that give them an essential survival signal. Therefore each B-cell passes a two-fold selection step: selection for recognizing an antigen, and selection for competitive binding to a helper T-cell. Binding to T-cells is specific and depends on the peptide presented on the MHC. B-cells that have a higher affinity for the antigen, have a larger internalization rate and a larger probability to succesfuly solicit T-cell help. Helper T-cell stimulation is essential for survival. One unverified hypothesis is that helper T-cells that have previously been trained in the thymus against self-proteins, now discriminate against hypermuations that lead to self-reactive B-cells. B-cells that do not receive a survival signal die.

Theoretical work oprea1997 ; Oprea2000 that was later verified by experiments has shown Victora:2012gx show that a single cycle of passing through the light and dark zones is not enough to generate the numbers of observed high affinity B-cells. This leads to a hypothesis called recycling where $\sim 90\%$ of B-cells go back to the dark zone. At the end of the affinity maturation process, B-cells increase their affinity to the antigen by $\sim 1000$ fold Berek:1987tf ; Oprea2000 (or even $\sim 10000$ in rabbits Eisen1964 ) and the B-cell population is always oligoclonal Tas2016 , with a large diversity of different B-cell receptors  Kuraoka2016 . Interestingly, a germinal center is binary: either it produces many high-affinity cells, or it fails to produce high-affinity B-cells Victora2012 .

The details of the heterogeneity of GC have been discovered through so-called brainbow experiments Victora:2010kd ; Tas2016 , in which individual B-cells were permanently tagged during BCR acquisition with one of 10 possible combinations of four different colored fluorescent proteins. This technique combined with imaging allowed for multicolor fate mapping of a cell’s progeny since all cells that come from the same ancestral cell have the same color. By injecting cells in germinal centers with different dyes and tracking their movement in mice, Victora et al Tas2016 have shown that in one lymphnode there are many germinal centers active at the same time, producing different diversities of cells. Most GCs probed at both 6 and 15 days showed a lot of different colored clusters, with lower estimates of $\sim 50$ clones per GC, going up to hundreds. By delaying tag formation until the cells are in the GC, selection in GCs was shown to keep multicolored clusters, showing that diversity of lineages is always maintained even if a certain dominance of one specific lineage (typically less than $40\%$ of all cells belong to the dominant lineage) is observed at later times. Notably the heterogeneity was huge: some GC had one lineage making up to $80\%$ of cells, others as little as $\sim 20\%$) two weeks after the infection. The rate of diversity loss was also very heterogenous: some GCs converged to a single lineage within a few days, and others took over two weeks, or never converged. Even in adoptive transfer experiments (see section V), where high affinity clones were introduced and the mouse was challenged with either bacteria or virus antigens, the GCs showed great heterogeneity regardless of the type of antigen - from nearly monoclonal to ones with an effectively neutral distribution of lineage diversity. The GCs expressing one lineage can be traced back to one single ancestral BCR (with possible SHM) that expands in a clonal burst over a short period of time, leading to a loss of diversity. Interestingly, the affinity of BCRs coming from GCs with varying levels of lineage diversity is similar. The oligoclonality of antibody repertoires is also confirmed by RepSeq experiments that show that antibodies responding to a specific antigen found at physiologically relevant concentrations can have between $3-147$ distinct CDR3s Lee2016 .

VIII.3 Evolution of broadly neutralizing antibodies

Many pathogens such as viruses (influenza, VIH) come in a variety of strains. Typically, antibodies mature to recognize just one type of strain, by targeting easily accessible epitopes on the surface of its proteins. However, viruses can often easily escape immunity afforded by these antibodies by mutating these epitopes, with little fitness cost. Epitopes that are more conserved (i.e. in which mutations carry a significant fitness cost to the virus) would be better targets, but they are usually harder to access by antibodies, as viruses have evolved to hide those conserved regions. To mature antibodies that bind strongly to more than one strain (called broadly neutralizing antibodies — BnAbs), the immune system needs to be trained with different antigens from different strains. In the context of a vaccination strategy, a crucial question is in what antigens to present, and in what order.

To address this question in the context of HIV, Wang et al Wang:2014vb simulated different temporal immunization schemes, by considering three antigens (wildtype HIV antigen and 2 mutants) and modelling their interaction with antibodies using a BCR-antigen binding model such as described in section II.4. The virus was modeled as consisting of a conserved part, in which mutations are deleterious, and a variable part. The three considered immunization schedules were: all three variants together (scheme 1), WT with mutant 1 at the same time, followed by mutant 2 (scheme 2), WT followed by mutant 1 followed by mutant 2 (scheme 3). Running these simulations many times, they looked at the antibodies that come out of these in silico germinal centers and scaned them against a standard panel of antigens that are used in affinity maturation experiments. The breadth of each antibody is defined by how many different test antigens it recognizes. Scheme 1 does not produce any broadly neutralizing antibodies (bnAbs). In fact, it hardly produces any germinal centers that have any output as B cells die quickly during the process. The frustration of not being able to discriminate between conserved and non-conserved residues results to a situation where the antibodies are not really selected since a BCR is likely seeing very different antigen in each round of selection. Scheme 2 does produce bnAbs, but with very low probability. This is likely what happens during a normal infection because you get infected with a single strain that later diversifies. In this case, the affinity maturation process in the first round does produce binders to the conserved residues, breaking the frustration. Finally, scheme 3 produces bnAbs with high probability ($\sim 69\%$). The sequential application allows the evolving BCR to focus on the only non-moving part of the target — the conserved residues. These predictions have been tested experimentally, showing that mice that were sequentially immunized focus their immune response.

There are effectively two parameters that govern the B-cell germinal center distribution: the probability of surviving selection $P_{s}$ and the probability of mutation $\mu$. With probability $P_{s}(1-\mu)$ the existing B-cells will expand, with probability $P_{s}\mu$ it will mutate, and with probability $1-P_{s}$ is will die. Wang et al Wang:2014vb suggest this system is optimally frustrated and the probability of survival puts the system in a special regime: too stringent selection (small $P_{s}$ will kill all cells), while too lenient selection (large $P_{s}$) will fill the germinal centers with existing cells, since $\mu$ is small. So only intermediate levels of $P_{s}$ lead to mutant cells surviving, allowing for many rounds of selection that results in acquiring many mutations that give large breadth. Support for this intermediate level selection comes from simulations that show that intermediate antigen numbers result in the largest number of mutations in a BCR. However experimental validation is still lacking.

Repeated antigen exposure helps increase the affinity of cells through selection, as in Darwinian evolution. In addition, new naive cells are recruited, adding novelty in the selection process other than through mutation. Murugan et al Murugan2018 compared the impact of selecting naive clones to that of producing new mutants by SHM. Looking at the response of BCR to the malaria parasite in humans coupled with string models of affinity maturaion, they showed that antigens with low complexity are much more efficient in generating good binders by means of SHM affinity maturation compared to high complexity antigens, which rely on recruitment of new naive cells. Both routes (SHM and selecting naive cells) where observed experimentally. Influenza vaccine studies show that, upon secondary immunization (booster vaccines), $\sim 7\%$ of BCR contained no SHM Neu2019 (although these experiments showed strong clonal dominance with $\sim 10\%$ of clones representing $\sim 90\%$ of sequences Lee2016 ). The model of Murugan et al Murugan2018 predicts that long exposure to small amounts of antigen leads to selecting for antibodies with a lot of SHM, whereas short exposure to high antigen concentrations results in selection on existing naive cell variation, which is confirmed in experiments. These results, which hold for relatively short timescales ($\sim$ few months), are consistent with the argument of Wang et al Wang:2014vb in the case of chronic infections. Interestingly, both Murugan et al Murugan2018 and Neu et al Neu2019 experimentally found that the majority of SHM do not necessarily improve antigen binding, and that both routes can lead to equally high affinity mutants, which supports the heterogeneity observed by Tas et al Tas2016 in the brainbow experiments. Yet, it is still not clear why some individuals fail to produce any good binding antibodies to the very strong antigenic challenge of malaria.

Additionally, the BCR repertoire changes with age, with both positive and negative selection acting differently in older people Ohare2018 , producing longer CDR3s that are more promiscious and likely to bind self-proteins. The rate of observed BCRs with SHM also changes from infancy to adulthood Wendel2017 . In influenza vaccine studies Neu2019 , people born before 1977 show less inter clonal diversity but achieve the same affinity as people born after 1977. Adaptation via the SHM route is less likely in older people, and they are more likely to use their larger existing memory pools than younger people whose antibodies acquire a lot of SHM. As a result, the antibodies of the older group target a conserved part (the stalk) of the influenza protein as do BnAbs, whereas younger people’s antibodies target the head (aiming for a strain-specific antibody). Younger people can produce BnAbs when hit with a strong challenge (vaccine combined with a pandemic), as predicted by the analysis of Wang et al Wang:2014vb , but the response to subsequent infections is still specific. In summary, these results highlight the complexity of BCR maturation, which involves both naive cells (with no SHM) and highly mutated cells. The details of the response depend on the dose concentration and timing, past infections, age and antigen complexity.

VIII.4 Population genetics approaches to affinity maturation

VIII.4.1 Evolutionary analysis of repertoire dynamics

Recent RepSeq data was used to quantify the evolutionary regime of in-host HIV evolution Nourmohammad2018 and in response to influenza vaccines Horns2019 . The flu vaccine is a model of an acute infection, which occurs on short timescales. BCRs from 9 timepoint blood samples, before and after the infection over the course of $\sim 2$ weeks, resulted in very skewed lineage trees and a U-shaped site-frequency spectrum (SFS). SFS corresponds to the distribution of frequencies of mutations in a population, and its shape is often informative about the underlying evolutionary process (see section XII.3.8). U-shaped SFS, as observed in BCR lineages after vaccination, are characteristic of strong selection. However, 3 out 5 of the strongly selected BCRs recognized neither the flu vaccine epitope nor the full virus, suggesting bystander evolution, by which non-responsive clonotypes expand in response to a multitude of signals. This could be a mechanism for upkeeping memory clones between infections with the same or similar pathogens.

Skewed trees and U-shaped SFS were also observed in the BCR lineages of HIV-infected indivuals Nourmohammad2018 . However, the same observations were made in out-of-frame sequences of healthy individuals. This is not surprising because the repertoires of healthy people are under constant selection, and out-of-frame sequences evolve in cells with functional receptors, and their evolution reflects the hitchhiking and selection on the whole cell. However, a finer analysis exploiting the dynamic trends of synonymous versus non-synonymous mutations is able to distinguish selection due to the chronic HIV infection from overall selection patterns. More interestingly, it finds that the CDR3 region of the BCR of untreated HIV carriers evolves according to a regime known as “clonal interference” (see section XII.3.9). In that regime, several new beneficial mutations arise at similar times and compete with each other. This competition slows down adaptation, as only one the beneficial mutations can survive. This analysis is based on estimating the probability that a (beneficial) mutation first rises to a threshold frequency $x$, but then is driven to extinction by a competing mutation, $H(x,x_{i})=G(0|x)G(x|x_{i})$, where $G(x|y)$, called a propagator is the probability of ever reaching $x$ starting at $y$, and $x_{i}$ is the initial frequency of the mutant. Both $H(x,x_{i})$ and $G(x|x_{i})$ can be estimated from the data and for various models of evolution. Neutral models cannot explain empirical observation, and neither does a simple model of selection with a fitness advantage. Instead, data can be fit by a model of varying selection (see section  XII.3.5), where the fitness advantage fluctuates as competing mutants come and go.

Analysis of phylogentic trees build from RepSeq data from HIV patients combined with S5F hypermutation models Yaari2013 also led to identifying the constraints on BCR adaptation Vieira2018 . Ancestral sequences in lineages are more likely to mutate their CDR3s than the framework (FWR). The propensity of a given residue to mutate was decreased more in framework (FWR) regions than in CDR3s, but in both cases a decrease in mutability was much more likely than an increase. Although most of the constraints on the residue to mutate in CDR3 that are under strong positive selection to increase their binding affinity, come from nonsynonymous mutations, in all regions up to $21\%$ of loss in residue mutability was caused by synonymous mutations, which are the result of neutral evolution. These results also point to a slow down in adaptation from clonal interference, as expected in the clonal interference regime.

VIII.4.2 Models of co-evolution of phenotypic traits

The fitness function for the co-evolution between HIV and BCR depends on the ability of the BCRs to recognize the virus, which is a (possibly nonlinear) function of the binding affinity. By focusing on a single antibody lineage and a viral population, Nourmohammad et al Nourmohammad2016 derived an effective stochastic equation for the rescaled mean binding energy $\epsilon$ between antibodies and viral epitopes:

[TABLE]

where $\sigma_{i}$ are the variances of the binding energy across the viral and antibody populations, $\theta_{i}$, $s_{i}$ and $N_{i}$ are the mutation, selection coefficients and population sizes of antibodies ($i=A$) and viruses ($i=V$); $\xi_{\epsilon}$ is Gaussian noise and times is measured in units of antibody coalescence time. This equation was derived by assuming a linear relationship between binding energy and fitness, and using a additive string model of binding energy as a function of genotype (see section II.4.2). The approximate method to get from the evolutionary dynamics, which are defined on the genotype, to an effective equation for the phenotype (binding energy), is described in a simpler context in  XII.3.10.

Intuitively, in general mutations on both antibodies and viruses reduce recognition. On the other hand, antibody diversity ($\sigma_{A}$) increases binding energy by selecting the best binders, while virus diversity ($\sigma_{V}$) decreases it by selecting the best escapers on the viral side of co-evolution. The same formalism can be applied to the evolution of BnAbs, which target epitopes on the conserved regions of the virus. In that context, the equations simplify to $\theta_{V}=0$ and $\sigma_{V}=0$, because the viral epitope is constant.

However, tracking the binding energy is not alone a good signature of co-evolution. Time-shifted statistics, such as the time delayed viral fitness (see Fig. 15), can be used to identify strongly co-evolving populations. The time delayed viral fitness is defined by the interaction between the viral population at time $t$, with distribution $y^{\gamma}(t)$ of genotypes $\gamma$, and a population of antibodies $\alpha$ taken at a later time $t+\tau$ distributed according to $x^{\alpha}(t+\tau)$: $F_{V,\tau}(t)\sim-s_{V}(t)\sum_{\alpha\gamma}E_{\alpha\gamma}y^{\gamma}(t)x^{\alpha(t+\tau)}$, where $E_{\alpha\,\gamma}$ is the binding energy from the string model between antibody $\alpha$ and viral epitope $\gamma$. Note that same-time fitness $F_{V,0}$ is equivalent to $-\epsilon$ up to a selectivity constant.

If the selective effects on the phenotype are comparable in the two populations, $s_{A}\theta_{A}\approx s_{V}\theta_{V}$, the current virus has highest fitness against antibodies from the near past since it has acquired mutations to escape them, and smaller fitness against antibodies from the future, which have caught up with these escape mutations (Fig. 15). For long times, mutations randomize the genome and $F_{\tau}(t\rightarrow\infty)=0$. If the selection pressure of viruses is larger than that of the antibody $s_{V}\theta_{V}\gg s_{A}\theta_{A}$, the virus is effectively evolving against a neutral antibody population. It still has a high fitness against previously seen antibodies, and genome randomization decreases its fitness advantage against future antibodies, but without a penalty.

The left slope of the time-shifted viral fitness quantifies the adaptation of the viral population to the existing antibody population and defines the a “fitness flux” $\partial_{\tau}F_{V,\tau}(t-\tau)|_{\tau=0-}=\phi_{V}(t)=\sum_{\gamma}\partial_{y^{\gamma}}F_{V}(t)(dy^{\gamma}(t)/dt)$. The right slope measures the pressure on the viral population from the adapting antibody population, which defines the transfer flux of fitness from antibodies to viruses $\partial_{\tau}F_{V,\tau}|_{\tau=0+}={\cal T}_{A\rightarrow V}(t)=\sum_{\alpha}\partial_{x^{\alpha}}F_{V}(t)(dx^{\alpha}(t))/dt)$. At stationarity both fluxes sum to zero, and the derivative is continuous, as in Fig. 15.

The characteristic S-curve shown for $s_{A}\theta_{A}\approx s_{V}\theta_{V}$ is indicative of two competing populations. Neutralization measurements Blanquart2013 ; Richman2003 ; Frost2005 ; Moore2009 show that viruses are more resistant to past antibodies and more susceptible to future antibodies, which results in the S-curved time delayed viral fitness Nourmohammad2016 . Interestingly, Blanquart and Gandon Blanquart2013 also showed that antibodies that evolved in one HIV positive patient are better at targeting the virus found in another person than this person’s current antibodies. This result remains unexplained.

While this analysis only involves a single antibody lineage, in germinal centers multiple lineages can compete with each other. Extending their analysis based on the dynamics of quantitative traits, Nourmohammad et al Nourmohammad2016 considered the competition of different lineages during affinity maturation. The change in the frequency of a given antibody lineage (each composed of many genotypes $\alpha$) of size $N_{A}^{C}$ in the population $\rho^{c}=\sum_{\alpha\in C}x^{\alpha}$ is driven by its mean fitness $F_{A}^{C}$ compared to the mean fitness of antibody lineages $F_{A}=\sum_{C}F_{A}^{C}\rho^{C}$:

[TABLE]

where $\xi_{C}$ is a Gaussian white noise. The lineage fitness is an average of the genotype $A^{\alpha}$ fitnesses $f^{\alpha}_{C}(t)$ that make up the lineage C, weighted by their frequencies in that lineage $x^{\alpha}_{C}(t)$, $F_{A}^{C}=\sum_{\alpha}f^{\alpha}_{C}(t)x^{\alpha}_{C}(t)$. The mean fitness of a lineage depends on the genotypes within this lineage and also on the frequency of virus $y^{\gamma}$ through $f^{\alpha}_{C}(t)$.

The probability of fixation $P^{C}_{\rm fix}$ of lineage $C$ is not only a function of the mean fitness advantage and population size, as in Eq. 114 (Sec. XII.3.4), but also of the ability of the population within the linage to adapt. Its value is given by:

[TABLE]

where $P^{0}_{\rm fix}=\rho^{C}(t=0)$ is the neutral fixation probability, and where $\phi^{C}_{A}$ and ${\cal T}^{C}_{V\to A}$ are the antibody counterpart of $\phi_{V}$ and ${\cal T}_{A\to V}$ for a specific lineage, and $\phi_{A}$ and ${\cal T}_{V\to A}$ their average over lineages. The first corrective term is standard and comes from the mean fitness advantage. The second corrective term accounts for the adaptation of each lineage in response to selection of antibodies for viral recognition within lineages. The third corrective term corresponds to the adaptation or escape capability of the viral population and its effect of the lineage.

For typical antibodies targeting a variable region of the virus, a broad viral diversity is detrimental as it allows the viral population to escape immunity (large ${\cal T}^{C}_{V\to A}$. For BnAbs, which bind to constant regions of the virus, this effect disappears. In fact the effect is the opposite: higher diversity increases the probability of BnAb fixation compared to a non-BnAbs antibody, whereas at low viral diversity the probability of fixation of both a non-BnAb and BnAbs is the same. This argument was proposed to explain why BnAbs show up later in the infection when the viral diversity is large.

IX Population dynamics of pathogens and hosts

IX.1 Viral fitness models

The immune system of a population of hosts defines a fitness landscape $F$ in which viruses evolve. Models have been developed to describe how this immune pressure shapes the evolution of both flu Luksza2014 and HIV Ferguson2013a ; Shekhar2013 ; Barton2016 at the population level. Viruses are usually taken to evolve according to traditional population genetics models.

This formalism can be used to draw short term predictions about the future fate of existing strains, e.g. of circulating influenza over the course of a year. Viral populations are subject to stochastic Wright-Fisher dynamics, but once they reach large enough numbers they grow exponentially following the fitness function (see section XII.3.3). In that deterministic stage, the size of a given viral strain $\alpha$, $X_{\alpha}(t)$ grows (or decreases) according to the fitness it experiences:

[TABLE]

Alternatively, one can consider the longer-term prevalance of particular strains of very diverse viruses such as HIV across many hosts, and use that to infer the fitness landscape of the virus. To do this, we can assume that the population reaches an “equilibrium” as the result of many mutations and fixation events across time in the entire population of viruses and hosts. The probability of finding a given viral strain is then assumed to follow Boltzmann’s law:

[TABLE]

where $\beta$ plays the role of an inverse temperature setting the tension between entropy from mutations, or genetic drift, and the fitness advantage. One possibility is to make use of the equilibrium assumption and learn the fitness Hamiltonian directly from existing viral data. This has been done for the Gag envelope protein of HIV Ferguson2013a using maximum entropy approaches (see section XII.5.1), such that the probability of seeing a given HIV strain $\alpha$ defined by its amino acid sequence $\vec{\sigma}$, is $p_{\alpha}=p(\vec{\sigma})=Z^{-1}\exp(\sum h_{i}\sigma_{i}+\sum_{i,j}J_{i,j}\sigma_{i}\sigma_{j})$, where $h_{i}$ define single-site fields that constrain the probability of observing the wildtype amino acid at a given site $i$, and $J_{ij}$ are interaction terms between amino acids (assuming a binary representation of the virus with $\sigma_{i}=1$ denoting the consensus amino acid and $\sigma_{i}=0$ a mutation at time $i$). Such models have been inferred from multiple sequence alignment of the protein, and then used to evolve in silico viral proteins in a population of host individuals Ferguson2013a ; Barton2016 , using quasi-species equations introduced by Eigen Eigen1971 ; Leuthausser1998 , to make predictions about HIV evolution.

Instead of learning the fitness directly from the sequence data, another strategy is to derive the fitness of each strain by decomposing it into two components

[TABLE]

where the $F^{\rm ag}_{\alpha}(t)$ is due to the antigenic component of the immune pressure exerted by the on the virus, and $F^{\rm stability}_{\alpha}$ encodes the fact that mutations that change the ability of the protein to fold should be detrimental to fitness. The fitness due to folding depends on the free energy of folding $G$, assuming a two state thermodynamic model of folding Phillips_textbook :

[TABLE]

where $\bar{G}$ and $G_{0}$ set the energy scale and are learned from data. Similarly, the antigenic component of the fitness is taken to be a sigmoidal function of the binding affinity $H$:

[TABLE]

where again $\bar{G}$ and $G_{0}$ set the energy scale and are learned from data. These constants are learned from data assuming that non-epitope mutations in the virus decrease protein stability, whereas epitope mutations decrease binding affinity. This approach was applied to influenza by Łuksza and Lässig Luksza2014 . The antigenic component $F^{\rm ag}$ was inferred from ferret blood titers, in which sera of blood from ferrets are challenged with different viral strains and their antibody response measured. To gain precision, the growth prediction of Eq. 82 were considered at the level of clades rather than individual strains by summing over all strains in each clade. In the next section we will see how this model could be used to predict the upcoming dominant strain of influenza.

IX.2 Co-evolution between host and pathogen populations

The adaptive immune system of hosts also evolves under the selective pressure of the antigenic environment, by expanding immune receptors that led to successful recognition. The antigens and the immune systems of individuals in a population engage on an arms race, where one forces the evolution of the other. The interaction between the immune system and different viruses takes on different forms. While HIV strains undergo very dynamic in-host evolution Richman2003 experiencing many mutations and strong selection accompanied by clonal interference Zanini2015 , the influenza virus evolves more slowly, and at scale of the whole population of hosts, with few mutations in the same host. We have discussed HIV evolutionary models in section VIII. Here we turn to the effects of population level co-evolution on immune repertoires, characteristic of viruses such as influenza.

This is currently an emerging field. Most of the treatments of immune systems so far have been coarse-grained and reduced to the effective selective pressure they exert on pathogens. This pressure is seen on antigenic maps, which place pathogens in a common space according to the similarity of the immune response to them (measured in the sera of ferrets). In practice, antigenic maps are produced using dimensionality reduction algorithms to reduce the response to a two dimensional manifold. These methods have been extremely useful in tracking influenza evolution Smith2004 and show that single point mutations can result in a completely new response (new faraway cluster), while some multiple mutations do not change the type of response (same cluster).

However, the link between the molecular interaction between immune receptor (both BCR and TCR, as introduced in section II.4.2), and the these phenotypic maps remains to be explored. An outstanding question is what evolutionary constraints these molecular details impose. For example, is the order of mutations important, and are all the mutations independent? Deep mutational scanning experimental techniques exist now to both find the best binding proteins in multiple rounds of selection experiments, and to map out the spectrum of possible solutions.

Recently researchers have very successfully predicted short-term flu evolution Luksza2014 ; Neher2014 . There are two successful types of models in this class. The first uses the diversity of existing strains and extrapolates the tree branches that have recently expanded to predict the dominant strains in the near future Neher2014 . This approach makes use of the coalescence framework described in section XII.3.7 and since it relies mainly on evolutionary properties, can easily be extended to other globally evolving viruses. The main idea is to statistically infer the characteristics of the evolutionary process from the existing influenza trees, with no reference to the immune system. The second Luksza2014 is based on molecular information about the influenza antigen Smith2004 ; Bao2008 and identifying successful mutations. In practice these are stochastic models that assess the probability of future strains. This method incorporates an effective treatment of the immune system, based on the response of ferret serum to flu strains, as explained in the previous section. In practice this treatment ranks influenza strains by the strength of serum response, but has no information about the molecular basis of the response, potential overlap for similar strains and its evolution. It then builds a fitness model, where the future frequency of a given viral strain depends on its structural stability and ability to escape the immune system. The details of these interaction models are learned from data using advanced statistical inference techniques. The second method relies much more on biophysical details and is therefore in principle more adapted to influenza, although the approach has been extended to predict the evolution of tumorous clones Luksza2017 as a function of their immunogenicity score (see Sec. II.4.4). Despite their differences both approaches are similar in style: they rely on recent evolutionary traces (from the last couple of months) to predict the dominant strain up to a year in advance. From the evolutionary perspective both models encode the idea of clonal interference between viral strains (see section XII.3.9).

Another theoretical approach to the viral-immune co-evolution problem is based on Susceptible-Infected-Recovered (SIR) Models. These models have been used for a long time to look at epidemic spreading within populations Grenfell2008 ; Perelson2002 . Many of these approaches are based on simulating sets of nonlinear equations. Originally these models were used to study slowly evolving viruses, such as measles Grenfell2008 . However recently, they have been tied with viral fitness models that account for the evolution of the virus. We already discussed one type of these models, for in-host evolution of HIV in section VII. Influenza, unlike HIV that evolves within the host organism, evolves mainly at the level of a population. Recent models Rouzine2018 ; Yan2018 of influenza evolution combine traditional SIR approaches with data-derived knowledge about flu strain evolution. Specifically, in Yan et al. Yan2018 , the model is of the form given by Eqs. 170-172 for each antigenic strain $a$, with the Susceptible equation solved explicitly, $S_{\alpha}\approx\exp({-\sum_{\alpha^{\prime}}K_{\alpha\alpha^{\prime}}R_{\alpha^{\prime}}})$, where $S_{\alpha}$ and $R_{\alpha}$ are the fractions of individuals who are susceptible and recovered from (and thus immune to) strain $\alpha$, respectively. $K_{\alpha\alpha^{\prime}}=e^{-|{\alpha-\alpha^{\prime}}|/d}$ defines a cross-reactivity Kernel of range $d$, where distances correspond to the number of mutations in an infinite-genome model. This equation means that individual who are outside of the cross-reactivity range of the recovered individuals are susceptible to be infected by a given strain.

Within this setup, a mutation in the virus introduces a new strain $\alpha^{\prime}$ that increases its distance to all existing strains by one. In a perturbative limit, this mutation increases fitness by a fixed and constant amount $s\propto d^{-1}$. This situation can be mapped onto a fitness wave model where many beneficial mutations compete against each other Desai2007a (see section XII.3.9). In the fitness wave description, the viral population has a distribution of fitnesses and reproduces according to that fitness. The dynamics are dominated by events happening at the stochastic nose of the distribution, i.e. strains with the the largest fitness $f_{m}$. The beneficial mutations that govern the future fate of the population all occur at that stochastic nose, with rate $r\sim f_{m}/\ln(f_{m}/\mu)$, where $\mu$ is the mutation rate.

The cumulative effect of these stochastic beneficial mutations at the nose can be described by an effective Langevin equation:

[TABLE]

where $\xi(t)$ is a Gaussian white noise of amplitude $\langle{\xi(t)\xi(t^{\prime})}\rangle{}=r\delta(t-t^{\prime})$, and where the negative fitness term $-I_{\rm tot}$ comes from the population of immune systems catching up with the viral strains in the bulk of the distribution, decreasing the fitnesses of all existing strains. This Langevin description makes it possible to estimate the rate of “speciation”, whereby two mutant strains at the nose escape the influence of the current immune systems in different way, creating two independently evolving antigenic niches which separate once their genetic distance reaches at least $\sim d$, $r_{\rm sp}\sim re^{-d/q}$, where $q=f_{n}/s$ is the typical number of beneficial mutations that strains at the fitness nose accumulate relative to the rest of the population.

These arguments can be used to explain the deep splits in flu strain topologies that have been observed in data and simulated using similar types of equations Bedford2012 ; Marchi2019 . In other words, the traveling wave picture also explains how influenza can constantly escape the immune system without continuous accumulation of genetic diversity.

Rouzine and Rozhnova Rouzine2018 applied a similar mapping onto a fitness wave description, but in a one-dimensional antigenic space. They used the predictions of their continuous traveling wave framework to estimate an antigenic mutation rate of $\mu\sim 3\cdot 10^{-5}$ per transmission event and predict the cross-immunity distance of $d\sim 15$ nucleotide substitutions that agrees well with independent estimates.

While these models give general insights and scaling laws to understand viral-immune co-evolution, they ignore the specific molecular details and mechanisms of immunity. As more data becomes available about the specific interactions between immune cells and viral strains, it may be possible to combine these approaches with data-driven models of immune repertoires and viral genomes, and to use them to make specific predictions about the fate of particular viral strains and their relation to immune repertoires of hosts.

X Discussion

Despite the length of this review, there are many subjects we did not touch upon that definitely fall into the global topic of quantitative immunology. Some of these, like the whole area of transcriptomics Regev2017 applied to immunology, are new and while the methodology is quantitative, the current experiments are only just starting to give quantitative models of the immune system that can be linked with a physical understanding (although things are moving so fast that by the time this review is published, this sentence may be obsolete). We also did not go into the methodology of certain analysis or theoretical approaches that are widely used, such as machine learning Mehta2018a or stochastic gene expression and biochemical regulation Walczak2012 ; Tkacik2011 , because detailed reviews for the physics audience exist for these topics. We refer the curious reader to these papers for the necessary background. Lastly, we only briefly mention some amazing experimental advances, such as imaging Victora:2012gx since they are currently in the process of being used to verify quantitative models. We hope we have managed to give the idea of a vibrant and multi-direction field. We also note that we made presentation choices, which were not easy, because many of the presented topics are linked to other ones. These links will surely become better explored in the coming years as solid experimental quantification and validation of theoretical models becomes the norm.

XI Glossary

Immunology can be painful at time for physicists because of its “jargon”. To help with accessibility, we try not to focus on specific molecules, but sometimes we need to. So here we give a brief overview of the main players and rules of the game for the immune system, and also other molecules that are mentioned in the text. For an introduction to immunology, we strongly recommend reading “Your Amazing Immune System” courtesy of the European Federation Immunological Society http://www.oegai.org/oegai/2-PDF/AmazingImmuneSystem.pdf AmazingImm_book . After this first glance, we recommend the short but illuminating book “How the immune system works” by physicist Lauren Sompayrac Sompayrac1999 .

• •

Adaptive Immune Response – a response of the organism to specific pathogens that changes in the lifetime of each individual. It is based on lymphocytes (T and B-cells) recognizing pathogens, proliferating and then keeping a subset of cells called memory cells that are adapted to previously encountered pathogens. Recognition involves interactions with many components of the innate immune system (e. g. antigen presenting cells, cytokines) and well as with different cells of the adaptive immune system. Adaptive immune systems first appeared in jawed vertebrates.

• •

Innate Immune Response – the innate immune response involves non-specific types of defense, from physical and chemical barriers (e.g. skin, clotting, scratching) to cells and molecules that recognize non-specific pathogenic patterns, such as proteins of the bacterial cell envelope (LPS - lipopolysaccharides). The innate immune system is evolutionary older than the adaptive one and is found in plants, insects, funghi, as well as vertebrates.

• •

Cytokines - small protein secreted by leukocytes to enforce cell-to-cell communications in the immune system. Examples of specific cytokines mentioned in the text: IL-2 is a key anti-apoptotic cytokine for T cells as well as an activation cytokine for regulatory T cells; IL-7 is an hematopoietic growth (anti-apoptosis, pro-proliferation) factor for lymphocytes.

• •

Interleukins (IL): cytokines historically thought to be produced by leukocytes.

• •

Leukocytes: the white blood cells of the immune system. These include lymphocytes (B, NK or T cells), granulocytes, monocytes and macrophages (see Fig. 9).

• •

TCR: T cell receptor, the main receptor on the surface of T cells that recognizes pMHC as a ligand for T cell activation.

• •

BCR: B cell receptor, the main receptor on the surface of Bcells that recognizes membrane-bound molecules as ligands for B cell activation.

• •

Antigen: biomolecules that trigger an adaptive immune responses. Antigens can be proteins that are recognized by antibodies (B-cell mediated responses) or short peptides that are loaded onto MHC (T-cell mediated responses).

• •

Epitope: part or whole of the antigen that the immune receptor bind to. For B cells, part of the protein that the B cell receptor bind to. For T cells, this corresponds to the peptide loaded by the MHC.

• •

APC (antigen presenting cells) – dendritic cells and macrophages. Surveilling cells that internalize molecules and cells from tissues and present them on their MHC type II to T-cells.

• •

MHC: major histocompatibility complex; its main function is to ”present” short peptides, as antigen for T cells. It comes in three types, but two of them are more relevant for the purposes of this review: type I and type II. Type I MHC are expressed on most cells in the body, presenting random bits of protein fragments found in that cell – this informs the surveilling cells if the MHC presenting cell is healthy or not. Type II MHC are expressed only on specialized cells (antigen presenting cells) and carry information about the cellular environment (whether there is an infection in a given tissue).

• •

CD8(+) (killer) T-cell – a type of T-cell identified by its CD8 marker involved in interactions with MHC type I presenting cells. These T-cells trigger apoptosis (kill by forcing the cells to kill themselves) the infected cells.

• •

CD4(+) (helper) T-cell – a type of T-cell identified by its CD4 marker involved in interactions with MHC type II presenting cells that present peptides from antigens they engulfed. These T-cells produce cytokines and help orchestrate the response of other immune cells (e.g. B-cells, killer T-cells).

• •

pMHC: a complex of a short peptide (p) and MHC that constitutes a ligand for TCR. Depending on the nature of the embedded peptide, a pMHC can be agonistic (triggering an immune response), antagonistic (extinguishing an immune response) or null. Such hierarchy often lines up with self/non-self discrimination.

• •

Somatic (hyper)mutation: genetic alteration acquired by a cell that occurs in body cells (somaplasm) and that can be passed to the progeny of the mutated cell in the course of cell division; one of the driver for antibody maturation in B cells, and for repertoire generation in T cells.

• •

Germline mutation: genetic alteration that occur in the gamete-producing cells (sperm and eggs)

• •

VDJ recombination – a DNA editing process that creates T- and B-cell receptors (see section VI.2 for details).

• •

Ig class – B-cells can express different types of constant regions of their receptors coupled to the same variable region. Throughout their lifetime, they can change (albeit in a specific order) which gene they express – this is called class switching. The expressed constant region gene determines the so-called class of the B-cell and antibody. The different genes are aligned in the immunoglobulin locus, and after the B-cell moves on the next class, previous classes get deleted. For example the $\mu$ and $\delta$ genes are first in the heavy chain locus and they lead to the expression of IgM and IgD chains, both of which are expressed on naive cells. The order of the remaining classes for the heavy locus is IgG3, IgG1, IgA1, IgG2, IgG4, IgE, IgA2. Class switching is regulated by cytokines, through regulation of gene regulation (see sections IV.2 and V.1).

XII Methods

This section is devoted to introducing or expanding the more technical details of methods relevant to quantitative immunology. These methods and concepts may be of different kinds. Some pertain directly to the biophysics of immunology and the molecular details of binding between receptors and cognate ligands. Some may be classical tools from other fields, from statistical learning to population genetics, but which carry some conceptual similarity with approaches from statistical physics. The presented methods are almost always cited in the context of recent work in physical or quantitative immunology. While they are often not established concepts of immunology, we believe that they might play an increasing role in its future development, with more and more imports from other fields.

XII.1 Physical kinetics

XII.1.1 Diffusion-limited reaction rate

To compute the association rates of 2 of two molecules (Ligand and Receptor), let’s place ourselves in the reference frame of the receptor, so that the receptor is immobile, but the ligand diffuses with coefficient $D=D_{\rm Ligand}+D_{\rm Receptor}$. The spatial distribution of the ligand concentration $C_{\rm Ligand}(\vec{r},t)$ at steady state obeys Fick’s law:

[TABLE]

with boundary conditons $C(\vec{r},t)\xrightarrow{r\rightarrow\infty}C_{\rm Ligand}$ and $C(r=R)=0$, where $R=R_{\rm Receptor}+R_{\rm Ligand}$ is the sum of the receptor and ligand radii (each modeled by a sphere). Because of spherical symmetry, $C(\vec{r})=C(r)$, we can consider only the radial part of the Laplacian in spherical coordinates

[TABLE]

and the steady-state solution is:

[TABLE]

The flux of ligands $j_{\rm Ligand}$ onto the receptor sphere is

[TABLE]

and integrating over the sphere surface one can compute total flux:

[TABLE]

The total number of collisions per unit time is $4\pi\left(R_{\rm Ligand}+R_{\rm Receptor}\right)(D_{\rm Ligand}+D_{\rm Receptor})C_{\rm Ligand}$ and the rate of diffusion-limited interactions between receptors and ligands $k_{\rm collision}=\Phi_{\rm collision}/C_{\rm ligand}$ is:

[TABLE]

XII.1.2 The rates of dissociation between two biomolecules

The physical chemistry of dissociation of two biomolecules can be modelled in two steps: a slow one that corresponds the breakage of chemical bonds between molecules, and a fast one that corresponds to the rushing in of the solvent. Estimating the rate of dissociation $k_{\rm off}$ thus corresponds to estimating the slow step, which is truly driven by the quantum physics of bond breakage. The rate is then encapsulated by the frequency of bond vibration ($k_{B}T/\hbar\approx 10^{12}s^{-1}$) multiplied by the probability of of successfully dissociating the molecular pair ($\exp\left(-\Delta G_{\rm off}/(k_{B}T)\right)$). This bond-breakage free energy sums up the contribution of all the bonds (electrostatic interactions, hydrophobic bonds, Van der Waals forces, hydrogen-bonds etc.) holding a molecular complex together: it could be computed from purely quantum considerations, independently of the context in which the pair is considered (solvated in the intercellular medium, embedded in the plasma membrane, in vacuum etc.). Yet, this remains a tricky proposition as small errors in estimating the bond energies will yield to exponentially-inaccurate estimates of $k_{\rm off}$.

XII.1.3 The formation of ligand-receptor pairs: equilibrium and kinetics

Once association and dissociation rates are available, estimating the kinetics of formation of the ligand-receptor pairs is done by considering the reaction:

[TABLE]

and integrating the equation:

[TABLE]

We have three species whose kinetics we are estimating, and two conservation laws (for the ligand and for the receptor), hence the reaction coordinate for this reaction is one-dimensional and can be solved analytically. We call the complex concentration $x=[{\rm Complex}]$”

[TABLE]

where $L_{\rm total}$ is the total ligand concetration, free and complexed, and $R_{\rm total}$ is the total receptor concentration. In steady state (for $t\rightarrow\infty$) we have $dx/dt=0$ and

[TABLE]

which has two formal solutions

[TABLE]

of which only $x_{-}$ is physical and satisfies $x<L_{\rm total},R_{\rm total}$.

The kinetics of relaxation

[TABLE]

can also be integrated:

[TABLE]

where the characteristic rate for the reaction is: $k_{\rm reaction}=k_{\rm on}\sqrt{(L_{\rm total}+R_{\rm total}+K_{\rm D})^{2}-4L_{\rm total}R_{\rm total}}$, and $a$ is a constant depending on the initial condition.

In the limit of a few receptors $R_{\rm total}\ll K_{D},L_{\rm total}$, the expression simplifies to $x=L_{\rm total}R_{\rm total}/(L_{\rm total}+K_{\rm D})$ and $k_{\rm reaction}=k_{\rm on}L_{\rm tot}+k_{\rm off}$. Three regimes can be considered:

• •

When the ligand exists in high concentration, $L\gg K_{\rm D}$ and $x=R_{\rm total}$, the system relaxes and fluctuates with a characteristic rate $k_{\rm reaction}=k_{\rm on}L_{\rm total}$.

• •

When the ligand is sparse and in limited amount, $L_{\rm total}\ll K_{\rm D}$, few complexes form, $x\approx L_{\rm total}R_{\rm total}/K_{\rm D}$, and the system relaxes and fluctuates with a characteristic rate $k_{\rm off}$.

• •

When $[{\rm Ligand}]\approx K_{\rm D}$, half of the receptors are occupied and half are free of ligands, $x=R_{\rm total}/2$; this is the mid-point of the dose response when one assess receptor occupancy against increasing concentrations of ligands.

XII.2 Gene regulation

XII.2.1 Basic model

In the simplest quantitative model of gene regulation, we can write the deterministic dynamics of gene transcription and translation as:

[TABLE]

where $p_{i}$ are the production rates and $\gamma_{i}$ the degradation rates of the molecules, for $i={\rm protein},{\rm mRNA}$. Since transcription and degradation rates for mRNA are often larger than translation and degradation rates for the protein, timescale separation leads to

[TABLE]

At steady state:

[TABLE]

This straightforward expression can become arbitrarily complicated. Both the production and degradation rates of mRNA and proteins of interest can be complex functions of the signalling response, and post-translational modifications and regulated degradation can also happen.

In simplest case, these rates can be approximated as constants, when adiabatic conditions apply such that signalling responses and cytokine production and consumption occur on timescales much shorter or much longer then gene regulation. Including positive feedback loops produces more interesting dynamical behaviour, especially when the reinforcement is through mutltimerized transcription factors that introduce additional nonlinearities.

XII.2.2 Auto-amplification with a single transcription factor

A simple feedback loop corresponds to the case when the rate of transcription of a protein species present in the cell at concentration ${\rm[protein]}=X$ is a dose response of the $X$ itself (e.g. when $X$ is a transcription factor that binds to the promoter region of the gene for $X$). Then the dynamic equation for $X$ is:

[TABLE]

where $k$ is the maximum production rate, $K$ the concentration of protein $X$ at half-maximum expression, and $\gamma$ is a degradation rate (we have used the separation of time scales to eliminate the mRNA stage, as in the previous paragraph). This equation has one fixed point $X=0$ for ${k}/{\gamma}<0$, and $X_{\rm equilibrium}={k}/{\gamma}-K$, for ${k}/{\gamma}>K$. A simple stability analysis around that non trivial fixed point, taking $[X]=X_{\rm equilibrium}+\epsilon$,

[TABLE]

shows negative restoring force ($K-{k}/{\gamma}<0$), meaning that all fluctuations get quenched and the fixed point is stable.

XII.2.3 Auto-amplification with multiple transcription factors

The complexity and relevance of the auto-amplification gene regulatory circuit becomes more relevant when the expression of gene $X$ is regulated by two transcription factors $TF_{1}$ and $TF_{2}$, i.e. both transcription factors need to bind to the promoter region of gene $X$ to elicit its transcription. We can compute the state of the promoter using classical tools of statistical mechanics. The promoter for gene $X$ can exist in four possible states with the related probabilities: unoccupied promoter ($p_{0}$), $TF_{1}$ only ($p_{1}$), $TF_{2}$ only ($p_{2}$), and both ($p_{12}$), with

[TABLE]

with the partition function:

[TABLE]

where $[TF_{1}]$ and $[TF_{2}]$ denote concentrations of the two transcription factors and $K_{1}=c_{0}\exp(\Delta G_{1}/RT)$, $K_{2}=c_{0}\exp(\Delta G_{2}/RT)$, $K^{\prime}_{1}=c_{0}\exp(\Delta G^{\prime}_{1}/RT)$ and $K^{\prime}_{2}=\exp(\Delta G^{\prime}_{2}/RT)$ their equilibrium binding constants related to the free energies of binding individually (non-primed) and cooperatively (primed) to the promoter binding sites.

When the cooperative binding of transcription factors binding to the promoter region is usually thermodynamically favored, $K^{\prime}_{1}\ll=K_{1}$ and $K^{\prime}_{1}\ll K_{2}$, the partition function simplifies to $Z\approx 1+[TF_{1}][TF_{2}]/(K^{\prime}_{1}K^{\prime}_{2})$.

A common case if when $[TF_{1}]=[TF_{2}]=[TF]$ and gene $X$ is transcribed upon homodimerization of a transcription factor in its promoter region (then $K^{\prime}_{i}=K^{\prime}$). The probabilities of the promoter to be unoccupied (transcribed at a basal level for auto-activation), $p_{\rm off}=p_{0}$, and to be occupied (transcribed at an enhanced level for auto-activation), $p_{\rm on}=p_{12}$ are:

[TABLE]

This form of regulation, which follows a Hill function (with Hill coefficient $h=2$), introduces a nonlinearity in the production of $X$ that is of critical relevance in cell differentiation. For auto-activation by homodimers, the dynamic equation for the regulation of $X$ is:

[TABLE]

where $k$ and $\gamma$ are production and degradation rates for $X$ respectively.

As in the non cooperative case, this dynamical system has a single fixed point if $k/\gamma<K$: $[X]=0$. But for $k/\gamma>K$, it has three fixed points at [math], $X_{-}$ and $X_{+}$. The stability of these solutions can be deduced graphically (see section V.1 for details): $X_{-}$ is an unstable fixed point and [math] and $X_{+}$ are stable fixed points. Hence, this simple system of gene regulation will generate two types of cells: cells that do not express $X$, and cells that express a high level of $X$. We see many examples of such bimodal distributions in expression of transcription factors, cytokines and surface markers in cells of the immune system (see section V.1).

XII.3 Population dynamics, genetics

To study the somatic evolution of immune clones and cell types it is useful to summarize some basic results from population genetics. The material in this section can be found in many evolution and population textbooks Hartl_Genetics ; Gillespie_book . To describe the evolution of a mutating population, one typically considers a background of genetically identical individuals (or cells, in the context of lymphocyte population dynamics). Each one of these individuals can acquire a mutation with rate $\mu$, which starts a new mutant subpopulation. Each individual in both the ancestral and the mutant subpopulations can die or reproduce in each generation, and mutants carrying non-neutral amino acid substitutions have a different (higher for beneficial and lower for deleterious mutations) growth rates. Since it is more likely that a mutation will confer a disadvantage than an advantage, most mutants are deleterious. But as we shall see, even beneficial mutations are not necessarily destined to succeed and ultimately the fate of most mutants is to die. Out of the lucky few that establish subpopulations of significant frequency, some will completely take over the population — they will fix. In this section we calculate the probabilities for fates of mutant subpopulations in the simplified case, compared to somatic evolution of immune repertoires, when there are no sources of new clones other than mutations.

XII.3.1 Deterministic mutation-selection balance

Let us start with a situation where two subpopulations exist, the ancestral clone of size $n_{2}$ individuals and a mutant clone of size $n_{1}$ individuals. The subpopulations grow (which accounts for both reproduction and death) with rates $\gamma_{1}$ and $\gamma_{2}$, respectively, where $g_{i}=r_{i}-f(n_{1},n_{2})$ is the balance of growth rate $r_{i}$ and a death rate $f(n_{1},n_{2}$. Individuals also mutate from the ancestral population to the mutant one with rate $\mu_{1}$ and mutate back with rate $\mu_{2}$:

[TABLE]

We also assume, as is often the case in population genetics models that the population size is constant, $N=n_{1}+n_{2}=\rm{const}$, which sets $f(n_{1},n_{2})=(r_{1}n_{1}+r_{2}n_{2})/N$. The constant population size constraint means that knowing $n$ the size of both subpopulations is completely determined by the fraction of individuals in the ancestral population $x=n_{1}/n$, which follows:

[TABLE]

where the selection coefficient $s=r_{1}-r_{2}$ describes how much faster (or slower if $s<0$) the ancestral population grows compared to the mutant subpopulation.

At steady state, in the absence of mutations $\mu_{1}=\mu_{2}=0$, either the ancestral subpopulation fixes, $x*=1$ for $s>0$, or the mutant one $x*=0$ for $s<0$. In the absence of selection but with mutations, the relative ratio of the mutation rates determines the fraction of each subpopulation, $x^{*}=\mu_{1}/(\mu_{1}+\mu_{2})$. With both mutation and selection a mutation-selection balance is established (for $s<0$) at $x^{*}=(1/2)(1+{2\mu}/{|s|}-\sqrt{4({\mu}/{s})^{2}+1})$ for $\mu_{1}=\mu_{2}=\mu$. When mutations are rare compared to the fitness advantage, $|s|\gg\mu$, this balance simplifies to $x^{*}\approx\mu/|s|$.

XII.3.2 Genetic drift

In order to observe a situation described by the deterministic equations in Eq. 109, the mutant population needs to grow to a sizeable fraction of the population. However, every mutant appears first in only one individual and undergoes a subsequent random walk of reproduction and death, which means that the most likely fate is for it go immediately extinct. This makes the effect of small number noise coming from finite population sizes, called genetic drift in population genetics, relevant. We can explore this effect considering a population that produces only neutral mutants (meaning all mutants grow at the same rate as the ancestor), and keeping the population size fixed to $N$ individuals. If we focus on one individual at some initial time, and follow its offspring, at very long times only two outcome are possible. Either its offspring have taken over the whole population, or the lineage has gone completely extinct. Since we start with $N$ individuals, the probability of taking over is $1/N$. The argument generalizes to a subpopulation of size $n$: the probability that one of its members has taken over the population is $n/N$.

XII.3.3 Wright-Fisher model

While we cannot calculate deterministically the fate of any particular mutant individual, we can calculate the probability of the evolution of the ancestral and mutant fractions, assuming two subpopulations and a constant population size as we did in subsection XII.3.1. Two models, the Wright-Fisher and the Moran models, each with slightly different setups, describe the neutral evolution of populations. The Wright-Fisher model assumes discrete, non-overlapping generations: at each generation, individuals from the previous generation are cleared, a new sample of $N$ individuals is drawn, each new individual picking an ancestor from the previous generation with probability $1/N$. Following a subpopulation of size $n$ in the parent generation, the probability that their offspring comprise $m$ individuals is given by the transition probability:

[TABLE]

The probability distribution for the size of that subpopulation evolves according to the recursion: $P_{t+1}(m)=\sum_{n}P_{n\rightarrow m}P_{t}(n)$, where $t$ labels generations.

Within this model the mean frequency of the subpopulation does not change with the time $\sum_{n}nP_{n}(t)=\langle{n(t)}\rangle{}=\langle{n(t+1)}\rangle{}$. Since the binomial distribution for large N becomes Gaussian, the evolutionary trajectory of the number of individuals in each subpopulation size is well described by a random walk: $n(t+1)\approx n(t)+\eta(t)\sqrt{n(t)(1-n(t)/N)}$, with $\eta(t)$ normally distributed.

The Wright-Fisher model can be generalized to include selection. The probability to pick a given member $i$ of the parent generation depends on its selective advantage or disadvantage $s_{i}$, as $(1/N)(1+s_{i})/(1+\bar{s})$, with $\bar{s}=(1/N)\sum_{i}s_{i}$.

In the case of two subpopulations, wildtype (of size $n_{2}=N-n$) and mutant (of size $n_{1}=n$), with selective (dis-)advantage $s$ for the mutant, Eq.  is modified to: $P_{n\rightarrow m}={{N}\choose{m}}p^{m}(1-p)^{N-m}$, with $p=(n/N)(1+s)/(1+\bar{s})$ and $\bar{s}=1+sn/N$. For large populations, this leads to a biased random walk:

[TABLE]

whose determistic part reproduces Eq. 109 in the absence of mutations.

XII.3.4 Probability of extinction

We can calculate the probability of extinction of the mutant subpopulation in the setup of two populations we have described until now. This limit corresponds to a small mutation rate: each subpopulation can fix or go extinct before a new mutation appears. While many immunological situations (e.g. affinity maturation) may not be in that regime, it is an important result to know. The probability that a subpopulation of size $n$ goes extinct, $q(n)$, requires knowing the probability that all their offspring, calculated within the Wright-Fisher model, will go extinct:

[TABLE]

given a selective advantage $s$ for the population. For large $N$ the binomial distribution is peaked, and we can Taylor expand $q(i)$ around $n$, $q(i)\approx q(n)+(i-n)q^{\prime}(n)+\frac{1}{2}(i-n)^{2}q^{\prime\prime}(n)$, which results in:

[TABLE]

Solving for the extinction probability with boundary conditions $q(n=0)=1$ (an extinct subpopulation remains extinct) and $q(n=N)=0$ (a fixed population cannot go extinct) yields $q(n)=1-\frac{1-e^{-sn}}{1-e^{-sN}}$. The probability of fixation is then:

[TABLE]

If $s=0$ we recover the result in subsection XII.3.2.

In the limit of strong positive selection ($Ns\gg 1$, but $s\ll 1$) the probability that a founder mutant ($n=1$) fixes is proportional to the selection strength $P_{\rm fix}\approx s$, if the initial population is one individual, or more generally $P_{\rm fix}\approx ks$, if the initial population is small $k\ll 1/s$ individuals. On the other hand $P_{\rm fix}\approx 1$ for large founder populations $k\gg 1/s$. This defines a threshold of $1/s$ individuals to ensure the survival of the subpopulation. For weak selection pressures $|{Ns}|\ll 1$, the probability of fixation is equal to the frequency of individuals in the population and is independent of $s$, $P_{\rm fix}\approx n/N$, making this regime effectively neutral. Lastly, for strongly deleterious mutations, $sN\ll-1$, the probability of fixation is non-zero $P_{\rm fix}\approx e^{s(N-n)}$, but is exponentially small.

XII.3.5 Moran model, continuous limit, and time varying selection

The above results can also be obtained considering a diffusion process in the number of individuals in the subpopulation. To do this, we use the Moran model, in which at each time step we choose one individual to die and one to reproduce (it can be the same individual), ensuring the population size is kept fixed. This implies overlapping generations, and the typical generation time is of order $N$ time steps.

Assume a wildtype population of size $N-n$ and a mutant population of size $n$. The probability of having $n$ at time step $i$ is given by:

[TABLE]

where the first term describes killing one dominant allele and reproducing a mutant, the second one reproducing a dominant allele and killing a mutant, and the last term describes the two possibilities of killing and reproducing the same allele. In the limit of large population sizes, rescaling time by the typical generation $N$, $t=i/N$, so that $p_{n}(i)-p_{n}(i-1)=(1/N)\partial_{t}p_{n}(t)$, and Taylor expanding in $x=n/N$:

[TABLE]

Mutation and selection can be added following a similar derivation and contribute in a mean “drift”:

[TABLE]

This Fokker-Planck equation is the stochastic version of 109, with effective diffusion coefficient $x(1-x)/2N$. This expression is also consistent with the random walk approximation of the Wright-Fisher model (111), meaning that the two models are equivalent in that continuous limit.

Eq. 117 has a general solution in terms of Gegenbauer polynomials. At steady state, assuming $\mu_{1}=\mu_{2}=\mu$ for simplicity it takes the form $p(x)=Z^{-1}[x(1-x)]^{N\mu-1}e^{Nsx}$, where $Z$ is a normalization factor. We can reinterpret it as a Boltzmann distribution, $p(x)=Z^{-1}p_{0}(x)e^{Nsx}$, where each mutant individual gives an energy gain $-s$, and $p_{0}(x)=\left[x(1-x)\right]^{N\mu-1}$ is the neutral distribution in the absence of selection.

Since Eq. 117 is a one dimensional diffusion equation, the probability of fixation is calculated from the backward equation for reaching the absorbing barrier at $x=1$. This calculation gives the same result as Eq. 114.

To model the fluctuations of the environment, we can consider a time varying selection pressure $s(t)=s_{0}+\sigma(t)$ with mean $s_{0}$ and random white noise fluctuations $\langle{\sigma(t)\sigma(t^{\prime})}\rangle{}=\delta(t-t^{\prime})\Omega$ around this mean Takahata:1975ug . The probability $p(y,x,t)$ that the wildtype, starting at initial frequency $y$ at time $t=0$, reaches frequency at least $x$ by time $t$ is given by the backwards equation:

[TABLE]

where

[TABLE]

and where the diffusion term cumulates genetic drift and fitness fluctuations:

[TABLE]

With the boundary conditions $p(y\geq x,x,t)=1$ and $p(y=0,x,t)=0$, the probability of ever reaching frequency $x$ is:

[TABLE]

where $\alpha_{\pm}=\frac{1}{2}\left[1\pm\sqrt{1+2/N\Omega}\right]$ and $\lambda=s_{0}/(\Omega\sqrt{1+2/N\Omega})$. For $x=1$ we recover the fixation probability reported in Takahata et al Takahata:1975ug , and for $\Omega\to 0$ we recover the result for constant fitness (114), $P_{\rm fix}=(1-e^{-s_{0}Ny})/(1-e^{-s_{0}N})$.

For more than two genotypes Nourmohammad2013 , Eq. 117 generalizes to an equation for a vector $\vec{x}=(x^{1},...,x^{K-1})$ of linearly independent genotype frequencies ($x^{K}=1-\sum_{\alpha}^{K-1}x^{\alpha}$):

[TABLE]

where the covariance matrices are:

[TABLE]

the mutation coefficients $m^{a}(\vec{x})=\sum_{\beta}(\mu_{\beta\rightarrow\alpha}x^{\beta}-\mu_{\alpha\rightarrow\beta}x^{\alpha})$ and the selection coefficient is the relative growth rate of genotype $\beta$ compared to a reference genotype $s_{\alpha}=f^{\beta}-f^{{\rm ref}}$.

XII.3.6 Branching processes

Branching processes are useful for tracking the fate of the offspring of a individual through time. We will first introduce it in the context of the Moran model, and then present its more standard applications. Assuming birth (with rate $1+s$) and death (with rate $1$) as the only possibly processes, we can track the evolution of the probability of having 1 individual at time [math] and $n$ at time $t$. This probability satisfies a recursion that can be obtained by considering the possible events occuring in the first time step, between times [math] and $dt=1/N$ Allen2003 ; Desai2007a :

[TABLE]

The first term corresponds to the lineage going extinct between times [math] and $dt$ — and remaining extinct until $t$. The second term corresponds to a division of the initial individual between times $t=0$ and $dt$, and the lineage then reaching size $n$ from size $2$ in the remaining time, $t-dt$. The last term corresponds to no change at all.

The next step, which characterizes the branching process approach, is to assume that the outcome of a lineage can be deduced from the outcome of each of the descendants of the first division, taken independently: $p(2,n,t)\approx\sum_{m=0}^{n}p(1,m,t)p(1,n-m,t)$, which assumes that once one individual gives birth to two, the births in these lineages happen independently. This approximation is valid as long as $n\ll N$. Then the recursion becomes in the continuous time limit:

[TABLE]

Defining the generation function $G(z,t)=\sum_{n=0}^{\infty}p(1,n,t)z^{n}$ we obtain:

[TABLE]

which is solved with boundary conditions $G(z=1)=1$ and $G(z,t_{0})=z$:

[TABLE]

For $s=0$ (no selection), we obtain $G(z)=1-\frac{1}{1+(1-z)^{-1}}$, whose series expansion yields $p(1,n=0,t)=t/(1+t)$ and $p(1,n,t)=t^{n-1}/(1+t)^{n+1}$ for $n\geq 1$.

For arbitrary $s$ the mean number of individuals at time t is $\langle{n}\rangle{}=\partial_{z}G(z)|_{z=1}=e^{st}$.The probability of going extinct at time t is

[TABLE]

which goes to $1-s$ for $s>0$ at $t\to\infty$. We thus recover the result that the fixation probability goes to $s$, as obtained previously in the strong positive selection limit. In that limit, by the time the mutant escapes genetic drift, for $n\gg 1/s$, its population size is still small compared to the total population size $n\ll N$, making the branching process approximation appropriate.

While we have illustrated branching processes in the context of a Moran model of evolution, branching processes are ubiquitous and used in a variety of contexts. One of the simplest branching process is defined in discrete time, where at each step a individual can divide with probability $p$, or not divide with probability $1-p$ Allen2003 ; Yates2007 ; Zapperi1995 . A recursion relation for the generating function of the total number of individuals $n$ in a lineage at time $t$ (past and present), $G(z,t)=\sum_{n}p(n,t)z^{n}$, can be written following similar arguments as above:

[TABLE]

This recursion equation is solved for $t\to\infty$ by:

[TABLE]

which has a critical point as $p\rightarrow 1/2$ Zapperi1995 , when the average number of offspring equals 1. Rewriting the generating function as a series we recover the probability of the total number of individuals (which must be even)

[TABLE]

For $p<1/2$, $P_{k}$ is a decaying exponential as expected, and for $p=1/2$ it decays as a power-law $P_{k}\sim k^{-3/2}$. Note that this class of critical branching processes has been used to explain power laws in the distribution of activity in neural networks Beggs2003 .

XII.3.7 Coalescence process

As seen when calculating fixation probabilities, it is often useful to think about evolutionary processes backwards in time. This is the basic idea behind a coalescence approach, which can be formalised using branching processes in the context of the Wright-Fisher model. Here we give just some basic intuition about how thinking back in time about the history of coalescing sequences in lineages can be useful when studying affinity maturation processes, or tracing phenotypic lineages. W will present the coalescence process in a neutral evolutionary framework. As before the neutral framework provides us with a null model in the case affinity maturation where selection is important, but it may also be useful in immunological phenotyping. The coalescence approach does not concern itself with mutations, but simply tracks genealogies. Mutations can later be added to an existing genealogy (tree).

We consider two individuals and ask how long ago they shared a common ancestor. If one individual has a given parent, the probability that the second cell has the same parent, given there are $N$ cells, is $1/N$. The probability that they do not have the same parent is $1-1/N$. Following this reasoning, the probability that they have the same parent $t$ generations ago, but not during the $t-1$ generations is

[TABLE]

where we have expanded for large $N$ and $T_{2}$ stands for time to mean recent common ancestor (MRCA). The mean time for two cells to coalesce is simply the mean expectation time of this distribution

[TABLE]

The mean time for two cells to coalescence is equal to the population size, in units of generation time. More generally, the probability that $k$ cells do not share the same parent is $P_{\rm diff}=\left[1-1/N\right]\left[1-2/N\right]...\left[1-(k-1)/N\right]\approx 1-{{k}\choose{2}}N^{-1}$. The probability that at least two cells have the same parent (or coalescence in our backwards picture) is then $P_{c}=1-P_{\rm diff}={{k}\choose{2}}N^{-1}$. If $N\gg k$ then the probability that more than two cells share the same parent in a single generation can be neglected and we will assume that in each generation only two cells will share a parent. The distribution of times until the first coalescence is:

[TABLE]

After a coalescent event, there are $k-1$ individuals left, and the process can be repeated until the whole genealogical tree is reconstructed. The coalescence probability completely determines the statistics of the topology of the tree and the branch lengths in a neutral process. Once the branch lengths (times between each coalescence events) have been determined, neutral mutations are distributed randomly along the branch with some rate, so that their number on each branch follows a Poisson distribution. Specificially, the number of mutations $\pi$ between two individuals, also called pairwise heterozygocity, is given by:

[TABLE]

for coalescence time $t$ and per-generation mutation rate $\mu$ (the total time is $2t$ because it adds the two branches of length $t$ from the common ancestor). Of course, all of this breaks down in the presence of selection.

Thus, the coalescence probability determines the genetic diversity within a population. The distribution of $\pi$ in the population can be obtained by integrating over the coalescence time (132),

[TABLE]

where $\theta=2\mu N=\langle{\pi}\rangle{}$. Thus, in principle, in the regime of neutrality, characterizing the distribution of the mean mutational distance should allow us to read off the mutation rate.

Departures from the Poisson distribution is one of the signatures of selection. In the presence of selection, the distribution of branch lengths depends on the details of the type of selection we are studying. If we do think about affinity maturation (where selection plays an important role), we see that cells undergo bursts of selection in the germinal centers, followed by periods outside of the germinal centers. Two cells sampled at the same time from the blood may therefore have very different recent histories and using their mutation distance to infer a mean mutation rate would be misleading.

XII.3.8 Site frequency spectra and tree balancing

Another statistics that is easy to calculate within the neutral model, and can thus be used as a null model to compare to data, is the site frequency spectrum (SFS). The SFS is the number of individuals in the population (in our case cells) that have a mutation at a given position in the aligned sequence with respect to the dominant base pair in the most recent common ancestor, and presents it as a histogram: ${(n_{1},n_{2},....,n_{N-1})}$, where $n_{1}$ is the number of mutations present in a single cell, $n_{2}$ is the number number of mutations present in two cells, etc. Sites that are not polymorphic (all have the same base pair) do not contribute to the spectrum as they coincide with the most recent common ancestor of the whole population. However a mutation at a site that is close to fixation will contribute to the spectrum, although this mutation now dominates the population. In absence of information about the most recent common ancestor, in population genetics the mutation is often called with respect to the outgroup. In the case of B-cell receptors, a good estimate can be the infered from the best alignment to the V, D, and J germline sequences.

We can estimate the SFS for a neutrally evolving population of constant size $N$ within the Moran model. In the continuous time limit, the mean number of mutations shared by $k$ individuals evolves according to (for $k>1$):

[TABLE]

with

[TABLE]

the “current” of mutations across subpopulation sizes. The first term corresponds to death event occuring in the subpopulation of size $k+1$ carrying the mutation of interest, and the second term to birth events in the subpopulation of size $k-1$. New mutation always starts with one individual, so that

[TABLE]

where $\mu$ is the mutation rate. At steady state the current is constant and equal to $J_{k}=J_{1}=\mu/N$, because each new mutation has a probability $1/N$ to fix, resulting in a current $\mu/N$ of mutations traveling from from size $1$ to $N$. This implies:

[TABLE]

which is solved by

[TABLE]

The SFS can also be calculated approximately in models with selection, but still with fixed population sizes, so we do not recall these results here Neher2013 .

Another feature that can be used to identify selection through departure from the neutral model is how balanced lineage trees are. Intuitively, a neutral process does not favor adding a new mutation to any of the branches, hence the resulting trees should be symmetric and balanced. A process with selection will preferentially grow the favoured parts of the tree, resulting in some long branches, and other short ones. In the neutral model, the number of leaves $n$ in a sublineage at generation $t$, follows the distribution $p_{t}(n)$, which satisfies the recursion:

[TABLE]

The sublineage of interest can add one leave by reproducing, with probability proportional to its size $n$, (second term), or not (first term). Solving by recursion gives a uniform branch length distribution $p_{t}(n)={1}/({n-1})$.

XII.3.9 Clonal interference

Another concept from population genetics that has been shown to be relevant in BCR affinity maturation is clonal interference, which we discussed in more detail in section VIII. Here we recall the back-of-the envelope arguments of Desai and Fisher Desai2007a to show on a more classical example why clonal interference slows down the rate of adaptation.

Given a population size $N$ and mutation rate $\mu$, a mutation occurs in any individual with rate $N\mu$, so the time between mutations is $(N\mu)^{-1}$. Most of these mutations go extinct due to genetic drift. In the strong selection regime, a mutation survives genetic drift with probability $s$. Thus, when mutations are rare, the average time for a new mutation to occur and fix is $(N\mu s)^{-1}$, and the rate of adaptation is given by $v_{\rm rare}=N\mu s^{2}$.

However fixation itself may take time, and the above picture breaks down when mutations are no longer rare, i.e. when the typical time between succesful mutation becomes smaller or of the same magnitude as the time for it to fix. After the subpopulation carrying the mutation overcomes genetic drift which results in an initial population size of $1/s$, we can consider that the subpopulation grows deterministically with rate $s$, $n(t)\sim(1/s)e^{st}$. The mutation will fix when $n(t)=N$, which gives the fixation time $t_{\rm fix}=(1/s)\ln{Ns}$. Mutation are no longer rare when $(N\mu s)^{-1}\sim 1/s\ln{Ns}$, or $N\mu\sim(\ln{Ns})^{-1}$.

In that regime, many mutants can co-exist at the same time, and new mutations can appear in existing mutants before these have time to fix. Mutants existing at the same time will have different growth rates, depending on the number of mutations they have acquired, even if we assume for simplicity that all mutations give the same fitness advantage. The whole population can be described by a fitness distribution and we can notice that the fate of a mutation in the bulk of this distribution (a subpopulation that was created some time ago) is different that at its nose (a recently created subpopulation). The nose subpopulations are small and susceptible to genetic drift. The bulk subpopulations have more individuals and grow deterministically, but are subject to nonlinear effects of competition between individuals. We can consider these two subpopulations separately and then “stitch” the two solutions.

The individuals at the nose of the fitness distribution have a fitness advantage, lets call it $qs$, meaning that they have $q$ mutations, each confering an advtantage $s$, compared to the bulk of the distribution. The time needed for an even fitter individual (with fitness $(q+1)s$) to appear is given by $\tau$ that satifies $\int_{0}^{\tau}dt\mu n_{\rm nose}(t)qs=1$, where $n_{\rm nose}(t)=1/(qs)e^{qst}$ is the number of individuals at the nose of the fitness distribution, following the same arguments as for the rare mutation case. In the limit of $e^{qs\tau}\gg 1$, one gets:

[TABLE]

Every time that a new fitter class is added to the distribution, the mean of the distribution also increases. After $q\tau$, the old nose of the distribution will become the mean (note that the individuals that are less fit than the mean go extinct deterministically). The populations in fitness class $q$ grow as $qs$, so in class $q-1$ they grow as $(q-1)s$, etc. A lineage originating at the nose grows to dominate the population $q\tau$ later, while progressively losing its fitness advantage. By the time it dominates the population, its relative fitness advantage is zero. On average, it will have grown at rate $qs/2$. This gives the consistency equation $(1/qs)e^{qst/2}=N$, with $t=q\tau$. This gives a second estimate of the establishment time:

[TABLE]

Equating Eq. 143 and Eq. 144, we find the average rate at which the population grows:

[TABLE]

We can solve this implicit equation approximately assuming that $q$ is not too big and neglecting the $\ln{q}$ terms:

[TABLE]

The rate of adaptation in the regime of multiple competing mutations is then:

[TABLE]

This rate scales with the logarithm of the population size, $v_{\rm CI}\sim\ln{N}$, much slower than in the rare mutation regime $v_{\rm rare}\sim{N}$. The clonal interference regime is not only a regime in which there are many competing mutations, but a regime where new mutations arise on the background of still relatively low frequency mutations, forming competing lineages. The distinction between competition of different clones and clonal interference is especially important in affinity maturation, where the competition in germinal centers could be between very different clones (which is not clonal interference), or between clones with similar histories, in which case it is clonal interference.

XII.3.10 Quantitative traits

While models of population genetics are often defined in the space of genotypes or alleles, what is often measured in the resulting phenotype, which is a (possibly nonlinear) function of the genotype. For instance, for lymphocyte receptors, the relevant phenotype may be defined as the binding affinity to epitopes of interest.

A projection of a description in phenotypic space, where mutations occur independently at all loci $i$ along the genome, results in an effective description in phenotypic space, which is useful because selection occurs at the phenotypic level and we often do not have access to the genotype–phenotype map. Given a fraction $x^{\alpha}$ of the population carrying allele ${\alpha}$, the population mean and variance of a phenotype $E$ is:

[TABLE]

where $N$ is the population size, and $w(E)$ is the distribution of this trait in the population.

One can derive an effective equation for the joint evolution of the mean phenotype and its variance Nourmohammad2013 :

[TABLE]

where we calculate the drift and diffusion terms below.

First lets focus on the diffusion terms, which correspond to genetic drift are are independent of selection:

[TABLE]

where we have used the definitions of the covariance matrices in high dimensional space from Eq. 123. Similarly

[TABLE]

where in the last step we have assumed $w(E)$ is approximately Gaussian. A similar calculation shows that $C^{\Delta\Gamma}=\langle{(E-\Gamma)^{3}}\rangle{}\approx 0$ with the Gaussian assumption.

To compute the drift terms, which arise from selective effects, we need to specificy how fitness and phenotype are related. For simplicity, let us assume that fitness of each individual is a quadratic function of the phenotype (or expand fitness close to its peak at $E^{*}$), $f(E)=-c_{0}(E-E^{*})^{2}$, where $c_{0}$ is a prefactor that measures the width of the fitness peak. The mean fitness is then:

[TABLE]

The change in the mean phenotype $\Gamma$ due to selection is:

[TABLE]

Since (by analogy with Eq. 122) $\frac{d\Gamma^{\rm sel}}{dt}=C^{\Gamma\Gamma}s_{\Gamma}$, using Eq. XII.3.10 we can verify that the selection coefficient is $s_{\Gamma}=\frac{\partial}{\partial\Gamma}F(\Gamma,\Delta)=-2c_{0}(\Gamma-E^{*})$.

In principle we can calculate the change in the mean variance due to selection in the same way $\frac{d\Delta^{\rm sel}}{dt}=\frac{d}{dt}\int dE(E-\Gamma)^{2}w(E)$ assuming a peaked phenotype distribution, but in practice its easier to use the analogy with Eq. 122 and Eq. 151:

[TABLE]

To consider the changes of $\Gamma$ and $\Delta$ due to mutations we have to write down a more detailed genotype model since mutations act on base pairs. If each locus can take a WT value ($\sigma_{i}=0$) or mutant value ($\sigma_{i}=1$), each linearly additive phenotype (for simplicity) $E$ (for example binding energy) can be written as $E^{\alpha}=\sum_{i}^{L}E_{i}^{\alpha}\sigma_{i}$, where $E_{i}^{\alpha}$ is a given sites contribution to the overall phenotype. We have $\Gamma=(1/N)\sum_{i,\alpha}E_{i}^{\alpha}\sigma^{\alpha}_{i}x^{\alpha}$ and:

[TABLE]

where $\Gamma_{0}=(1/2)\sum_{i}E^{i}$ is a phenotype average.

Similarly, the variance defined in Eq. 148 is also a sum over sites, and the change in the variance due to mutations is:

[TABLE]

where $E^{2}_{0}=(1/4)\sum_{i}E^{2}_{i}$.

Now we have all the elements, the effective diffusion equation (XII.3.10) in phenotypic space reads:

[TABLE]

The stochastic equation for the evolution of the mean phenotype can then be read off as:

[TABLE]

where the first term describes mutations, the second term selection and $\partial_{\Gamma}F=-2{c_{0}}(\Gamma-E^{*}(t))$ and $\xi_{\Gamma}$ is a normalized Gaussian white noise. $E^{*}(t)$ can be a time dependent moving fitness maximum. An analogous equation holds for $\Delta$. Assuming that $\Delta$ changes on much faster timescales than $\Gamma$, it can be replaced by its mean, and the steady state solution has a Boltzmann form $Q_{\rm eq}(t)=Z^{-1}Q_{0}(\Gamma)e^{-2NF(\Gamma)}$, where as before Q_{0}(\Gamma)\sim\exp\Big{[}-\frac{1}{2}\frac{\left(\Gamma-\Gamma_{0}\right)^{2}}{\langle{\Delta}\rangle{}/(4N)}\Big{]}) is the distribution with no selection.

In a time dependent environment the peak of the distribution also changes with time according to a prescribed model. The population tries to track the fitness peak, without really ever reaching it. From the physical point of view, the system is maintained out of equilibrium. In a changing fitness landscape, where the population history is described by a sequence of phenotypic trait measurement $(\Gamma_{0},....\Gamma_{M})$ over time $(t_{0},...,t_{M})$ the fitness flux of a population history $\Phi$ Mustonen2010b describes a cumulative selective effect of phenotypic trait changes:

[TABLE]

The lack of equality holds in general because $F$ also has an explicit time dependency. One can show that this discrepancy is formally equivalent to dissipation in thermal systems, and can be related to the entropy production:

[TABLE]

where ${\cal P^{T}}$ is the probability of the forward trajectory and ${\cal P^{T}}$ of the backward trajectory and the Kullback-Leibler divergence is defined in Eq. 227. This dissipation has been called “fitness flux” in the context of population genetics.

These kinds of phenotypic trait models have been used as starting points for studying co-evolution as we describe in section VIII. The main difference is that the fitness objective $E^{*}$ is itself a function of the composition of the population with which the initial population evolves. That other population is also subject to selection and drift and evolves stochastically, giving rise to coupled equations Nourmohammad2016 .

XII.3.11 Lineage reconstruction

Lineage reconstruction is a necessary first step in both BCR evolutionary analysis and phenotypic tracking. While many software methods exist to reconstruct lineages for evolutionary problems, they are not always well adapted for immunological data. Nevertheless they are often used, and more adapted methods are usually built upon classical ones, so here we provide a general overview of these existing approaches.

The first problem in reconstructing lineages requires taking the sequence data and identifying clusters of sequences that share a common ancestor, and therefore belong to the same lineage. A classical strategy to cluster datapoints is single-linkage clustering, which builds hierarchical clusters by iteratively merging pairs together Everitt2011 . The Partis software Ralph2016a uses a likelihood ratio test to decide if a give set of sequences $\sigma_{1},\sigma_{2},...,\sigma_{N}$ can be grouped as descending from one ancestor or not. Specifically Partis uses a Hidden Markov Model (HMM) based method to annotate each nucleotide in a specific BCR sequence as coming from a given a set of hidden states corresponding to V, D or J genomic template, or from a set of exponentially distributed non-templated insertions. The sum over paths determines the probability $P(\sigma_{i})$ for each sequence. The same procedure repeated for a set of $N$ sequences, results in the total probability of a common scenario for these sequences $P(\sigma_{1},..,\sigma_{N})$. If the likelihood ratio $P(\sigma_{1},..,\sigma_{N})/\Pi_{i=1}^{N}P(\sigma_{i})>{\rm threshold}$, one concludes that the sequences originate from a common recombination. Partis starts with pairs of sequences, keeping the cluster with the largest likelihood ratio, and then adding new members to the cluster based on the same test. In practice, to speed up the algorithm Partis does not test all possible pairs but creates initial subsets of data with low Hamming distances. Joining a cluster is irreversible, which can lead to errors in clustering.

Once a lineage is defined, annotation softwares such as Partis Ralph2016 or IGoR Marcou2018 can be used to propose the naive root of the tree.

This problem alone is simpler than finding the whole tree genealogy, known as the topology of the tree or the branching pattern, because part of the sequence is templated by the V, D, and J genes (see Sec. VI.2). There are two main classes of methods for tree topology reconstruction: maximum parsimony and maximum likelihood. Let us first explain why it is practically impossible to exhaustively sample all tree topologies, and then explain the differences in the two approaches.

We will use the example of binary trees — trees in which each node gives rise to only two branches — since essentially all loopless trees can be cast into a binary form by adding branches of zero length. Starting from the simplest unrooted tree of two leaves connected by one branch, and adding new leaves one at a time, it is simple to convince oneself that in each step we add $1$ leave, $1$ internal node and $2$ edges, such that a tree with $N$ leaves has $N-2$ internal nodes and $2N-3$ edges (or branches). Adding the $N$th leave ($N-1\rightarrow N$) adds $2N-5$ topologies, such that the number of unrooted tree topologies with $N$ leaves $T_{N}$, is $T_{N}=(2N-5)T_{N-1}$, which can be recursively solved to give $T_{N}=(2N-5)!!$. For rooted trees, the number of rooted topologies with $N$ leaves is given by $T_{N+1}$ in terms of the number of non-rooted tree topologies, since it just requires adding a root to the existing unrooted tree topologies.

The goal of tree reconstruction is to find the topology that is consistent with the data. For a tree with $N=10$ leaves, that means exploring $T_{N}\sim 10^{6}$ trees and for $N=50$, $T_{N}\sim 10^{76}$ trees, which is prohibitively large. BCR clusters often have hundreds of leaves.

There are two ways of determining distance between two node sequences along an edge of a tree: one is to calculate the Hamming distance by simply counting the number of single nucleotide differences between two sequences $h_{j}$, the other is to use the estimated time between two mutations $t_{j}$. The number of mutations in a given time $t_{j}$ is distributed according to a Poisson distribution of mean $\mu t_{j}$: $P(h_{j})=\frac{\mu t_{j}}{h_{j}!}e^{-\mu t_{j}}$. The branch length refers to the time between two nodes. The tree length is defined as the total number of substitutions, $h_{\rm tot}=\sum_{\rm{branches}}h_{j}$, and is not equal to the sum of branch lengths. Maximum parsimony methods use Hamming distance to define distances on trees, which requires knowing the identity of all the internal nodes and the tree topology, whereas maximum likelihood sums over possibilities for the sequences on the internal nodes, but requires a mutational model.

The maximum parsimony approach is easy to formulate: for a given topology, it assigns internal nodes so as to minimize $h_{\rm tot}$. It then selects topologies with minimal $h_{\rm tot}$. The parsimony scores differ from [math] (when all the leaves are the same) to $NL$ (when all the leaves are different, which essentially means this position carries no information). There are $(2N-5)!!$ topologies to scan so many trees have the same $h_{\rm tot}$, including at its minimum. Therefore there is often a large family of most parsimonious trees. Apart from the practical problems, the assumptions behind the evolutionary model in this approach are not clear ZYang2006 . In practice, the reconstructed trees match the true structure if there is little variation in the internal nodes. The main advantage is that it embodies the basic intuition for reconstructing the phylogeny.

In the maximum likelihood method, given an overall lineage age $T$ and tree topology, the likelihood is given by a specific mutation model by varying the branch lengths $t_{j}$. Within the independent site assumption, the likelihood factorizes over sites and the log-likelihood $\ell=\sum_{r=1}^{L}\ln{\ell_{r}}$ independently traces the evolution of each site $r$. Thus the following equations apply to a single site, but generalize readily to many sites. The probability of mutating the base pair from $x_{i}$ into $x_{j}$ between nodes $i$ and $j$, is defined as $P_{x_{i},x_{j}}(t_{j})$. For instance, this probability can take the following form, although many others are possible:

[TABLE]

where $\pi_{x}$ (with $\sum_{x}\pi_{x}=1$) describes the probability that a mutation results in base pair $x$, and is equal to $1/4$ for the unbiased 4-base pair model, and $t_{j}$ is the branch length (time between nodes $i$ and $j$). The first term corresponds to no mutation, and the second to a mutation (including one into the same base pair). This formulation guarantees detailed balance: $\pi_{x}P_{xy}(t)=\pi_{y}P_{yx}(t)$ for all $t$.

We can now build the likelihood recursively, assuming a fixed root identity [math], with sequences denotes as $x_{0}$ that are not fixed. First consider a tree made of a root and two branches of length $t_{1}$ and $t_{2}$ leading to leaves $1$ and $2$, whose identity $x_{1}$ and $x_{2}$ is known. The likelihood of this tree is:

[TABLE]

We can now add an internal node $4$ to the tree, descending from [math] and ancestor of $1$ and $2$, and a leave $3$ descending from 0:

[TABLE]

where $L_{4}(x_{0})=P_{x_{4}x_{1}}(t_{1})P_{x_{4}x_{2}}$. This recursive form allows for efficient calculation of likelihood by successively joining roots of trees into a new root. Formally, joining root $j$ and $k$ into common parent $i$:

[TABLE]

Leaves are initialized to $L_{j}(x_{j})=1$, and the final results is given by:

[TABLE]

Given a fixed topology, we can now maximize the likelihood $L$ over the set of intermediate branch lengths $\{t_{j}\}$, which results in a score for each topology. The topology ${\cal T}$ with the best best ranking score, is the most likely topology (similarly to the maximum parsimony approach).

Existing maximum likelihood methods use the independent site assumption to decrease computational time. Many methods also assume homogeneous mutation rates across sites and the reversibility of the evolutionary process. The GTRGAMMA substitution model of the RAxML software Stamatakis2014 uses gamma-distributed mutation rates for different base pairs, and IgPhyML Hoehn2017 encodes a non-reversible mutation model that effectively accounts for the context dependence of mutations (although in a site-independent setup).

Davidsen and Matsen Davidsen2018 compared maximum likelihood and maximum parsimony methods for BCR lineage reconstruction. They concluded that an improved and informed maximum parsimony method outperforms classic maximum likelihood.

XII.3.12 Population growth rates

Another way to study the evolution of populations is to study dynamics that maximize the long term growth rate of the population given a fixed environment Kussell2005 ; Rivoire2011 :

[TABLE]

where $N_{t}$ is the population size at a given time, whose evolution is driven by the dynamics. One needs to specify classes of evolutionary dynamics to consider, and the class of interactions it has with the environment. The framework is clearly very general, and reduces to a constrained optimization problem. The main conceptual difference with population genetics models studied in this section is that populations cannot go extinct, and all possible variants are represented. Populations are in a regime where fast growth governs evolution. Alternative approaches could consider optimal strategies for avoiding extinction Bradde2019 , by minimizing the probability of extinction, $P_{\rm ext}$, whose form depends on the problem considered.

XII.4 Ecological models

Population of lymphocytes interact with other through signaling, competition for resources, cytokines, or antigens. For this reason, concepts and mathematical models from ecology are often useful to describe their dynamics.

XII.4.1 Generalized Lotka–Volterra models

Ecological generalized Lotka-Volterra models  Lotka20 ; volterra56 ; Keshet2005 describe the co-habitation of multiple species in the same environment and account for their interactions, both direct (one species needs another to reproduce) and indirect (through competition for external resources). In general these models have the form:

[TABLE]

where $x_{i}$ is the frequency of a given species, $c_{i}$ are external factors such as nutrients $\alpha(\vec{x},\vec{c})$ is the growth rate and $\beta(\vec{x},\vec{c})$ the death rate with $\vec{x}=(x_{1},x_{2},\ldots)$ and $\vec{c}=(c_{1},c_{2},\ldots)$. The form of the dependence of the growth rate on the different species and external factors determines the non-linearity of the problem. This is usually encoded by an interaction matrix between species. Models where the growth and death rates are the same for all species are called neutral. Often one introduces a carrying capacity which is a form of non-linearity that describes competition for resources.

This formulation of the model is very general and the solution depends on the specific assumptions. In general these models are solved numerically, and the results depend on the numerical values of the parameters that are often not known with great certainty. These models are often used to ask questions about the coexistence of different species, as well as speciation itself — why are ecological environments with different species stable? Since we often do not have experimental information to parametrize the interaction matrix, random matrix models have been succesfully used to effectively describe the interactions May2001 . In this approach one uses the fact that within a family of random matrices (i.e. matrices whose elements are chosen from the same distribution) the eigenvalues of these matrices are the same. Near the fixed point, the stability of the system is explored by linearizing the system of equations in Eq. 167, and the eigenvalues of this linearized matrix determine the stability fo the ecosystem:

[TABLE]

where $\epsilon$ is the interaction strength and $\bf{K}$ is a random interaction matrix. As a result the stability of the system does not depend on the realization of the interaction matrix, just on its statistics. For the matrix to be stable, all the eigenvalues of $\bf{K}$, $\lambda_{i}$ have to satisfy $\epsilon\lambda_{i}\leq 1$, which is fulfilled in the largest eigenvalue satisfies $\lambda_{\rm max}<1/\epsilon$. In the case of Gaussian random matrices, the properties are determined by their first two moments and a strong transition occurs for $N\rightarrow\infty$ where the system is stable for $\epsilon<1/\sqrt{2}$ and otherwise unstable. This result holds only if the connectivity scales with the size of the system.

Incidentlly, random matrix theory also has been used when looking at covariances in sequence variation (see section XII.5.6). The general idea is that beyond the first couple of eigenvalues of the covariance matrix of amino acid variability (which often stems from phylogenetic bias), eigenvalues are well approximated by a random matrix. For $R\times C$ Gaussian random matrices, where the elements $K_{ij}$ are chosen from a Gaussian distribution with mean [math] and variance $\sigma$, the density of eigenvalues of the covariance matrix $Y=\langle{KK^{T}}\rangle{}/R$ is given by the Marchenko-Pastur distribution MarchenkoPastur1967 :

[TABLE]

where $\theta=C/R$ and $\theta_{\pm}=\sigma^{2}(1{\pm}\sqrt{\theta})$.

XII.4.2 Susceptible-Infected-Recovered (SIR) Models

Susceptible-Infected-Recovered (SIR) models are used to study the spread of epidemics at the population level. An SIR model considers three possible kinds of individuals: ones that are susceptible to the infection (S), infected (I) and recovered (R) and therefore usually immune. The simplest SIR model assumes that the birth-death of the host individuals happens on slower timescales than the spread of the epidemic itself. This model is a good description of measles and other infections by slowly evolving pathogens.

[TABLE]

where $\beta$ is an effective infection rate, $\nu$ is the rate of recovery, and $\nu(\beta S/\nu-1)=x$ defines an effective growth rate, or fitness of the host population. Within this formulation no one can die. Since $\partial_{t}S+\partial_{t}I+\partial_{t}R=0$, the total population size is constant, $N=S+I+R=$const. The typical recovery time is $T_{\rm r}=\nu^{-1}$ and the typical timescale for a infected-susceptible individual interaction is $T_{\rm inf}=\beta^{-1}$, meaning that an infected individual manages to infect $R_{0}=T_{\rm r}/T_{\rm inf}=\beta/\nu$ individuals before recovering. $R_{0}$ is also called the infection radius or basic reproduction number.

Using Eqs. 170-172, we find $S=S(t=0)e^{-R_{0}(R-R(t=0)}$. At infinite time ($t\rightarrow\infty$, $I=0$), the number of recovered individuals reads:

[TABLE]

which means that at an end of an epidemic there are still susceptible individuals. From Eq. 171 we find that $S(t=0)R_{0}>1$, leads to $\partial_{t}I>0$ and an exponentially growing infection (i.e. epidemic outbreak), whereas $S(t=0)R_{0}<1$ leads to $\partial_{t}I<0$ and quenching of the infection. Since these equations are deterministic, the initial condition determines the infection spread for all times. The effective number of susceptible individuals that can propagate the infection determines its future.

SIR equations can be extended to account for host birth (with rate $\alpha$) and death (with rate $\mu$), that lead to equations of the form:

[TABLE]

with two stable steady state solution: pathogen free, $(S,I,R)=(\alpha/\mu,0,0)$; and permanent infection, $(S,I,R)=\left((\nu+\mu)/\beta,(\mu/\beta)(R_{0}-1),(\gamma/\beta)(R_{0}-1)\right)$, where $R_{0}=\beta\alpha/\left(\mu(\mu+\nu)\right)$ is the infection radius. For $R_{0}\leq 1$, the system evolves towards the pathogen-free solution, while for $R_{0}\geq 1$ the permanently infected solution is reached.

The basic SIR model can be made more complicated in many different ways. One other variation called the SIS (Susceptible – Infected – Susceptible) model describes a situation where an infected individual does not have long lasting immunity. This model is inspired by influenza, where the virus mutates fast enough that a past infection does not necessarily guarantee immunity (although current models of influenza account for more detailed descriptions of cross-reactivity – see section IX). The model just has two states:

[TABLE]

Using $S=N-I$, we obtain a single differential equation:

[TABLE]

which has two fixed points $I=0$ for $R_{0}=(\beta/\nu)N<1$ and $I=(\beta N-\nu)/\beta$ for $R_{0}>1$. The second fixed results from the solution:

[TABLE]

We note that the SIR equations (Eqs. 170-172) differ in its assumptions from the two species Lotka-Volterra equations (see section XII.4.1) between and a host and a pathogen. Specifically, SIR equations do not assume that everyone can get infected, while Lotka-Volterra equations do. However, generalized Lotka-Volterra equations can easily be modified to account for different subclasses of hosts. Lotka-Volterra equations on the other hand do explicitly model the viral population, and allow for host-pathogen oscillations that SIR models do not, and they do do not assume a constant population size. Notably, the SIS equations (Eq. 177) are a two-species realization of Lotka-Volterra dynamics, with the additional assumption for constant population size.

In section VII we describe equations that model the in-host evolution during HIV infections. They are very similar in spirit to SIS equations which explicitly consider viral dynamics.

XII.4.3 Solution of stochastic population dynamics with a source

Here we present the calculations leading the results presented in Sec. VII, some of which are original to this review. To solve Eq. 62, we define the generating functions:

[TABLE]

The time evolution of the population is then governed by:

[TABLE]

At steady state, we have:

[TABLE]

where we have used the balance between birth and death of clones $dG(1)/dt=\Theta(1)-\nu N_{1}=0$. Using $(\Theta(1)-\Theta(x))/(1-x)=\sum_{C>0}\theta_{c}\sum_{i=0}^{C-1}x^{i}$, and expaning $1/(\nu-\mu x)=(1/\nu)\sum_{j=0}^{\infty}(\mu/\nu)^{j}x^{j}$, we obtain:

[TABLE]

which reduces to

[TABLE]

in the simple case $\theta_{C}=\theta\delta_{C,1}$. The total number of clones reads:

[TABLE]

and

[TABLE]

for $\theta_{C}=\theta\delta_{C,1}$. In the limit of rare divisions $\nu\ll\mu$, only clone sizes below the initial clone size, $C\leq k$, contribute, so that

[TABLE]

where $H_{k}=\sum_{i=1}^{k}1/i$ are harmonic numbers. In the opposite limit of division balancing death, $\mu-\nu\ll\mu$, the long tail of the second sum dominates:

[TABLE]

with $\epsilon=(\mu-\nu)/\mu$.

Eq. 183 can be solved out of steady state with the method of characteristics. Defining $F(x)=G(x)-G(1)$ and making the change of variable $y=\ln[(x-\mu/\nu)/(x-1)]/(\mu-\nu)$, this equation reduces to:

[TABLE]

This equation can be solved in absence of a source term ($\Theta=0$): $F_{0}(x,t)=A(t+y)$, with $A(y)=F_{0}(x,0)$, yielding:

[TABLE]

Starting with a single clone of size $s$ at $t=0$, $G_{0}(x,0)=x^{s}$ and $F_{0}(x,0)=x^{s}-1$, the coefficients of the solution $G_{0}(x,t)$ to the homogeneous equation (Eq. 191 with $\Theta=0$) give Green’s function, defined as the probability that the clone starting at size $s$ at $t=0$ has size $C$ at a later time $t$, $G_{0}(x,t)=\sum_{C}P(C,t|s,0)x^{C}$. Expanding (192) in $x$, one obtains:

[TABLE]

with the shorthands $z=e^{-t(\mu-\nu)}$, $r=\nu/\mu$. The combinatorial factor in the sum comes from counting the number of ways there are to choose $s-n$ factors of order [math] in $x$, $\binom{s}{n}$, and $n$ factors of order $(j_{1},\ldots,j_{n})$ so that $\sum_{i=1}^{n}j_{i}=C$, $\binom{C-1}{n-1}$, in a series expansion of the form $[a+b((cx)+(cx)^{2}+(cx)^{3}+\ldots)]^{s}=[a+bcx/(1-cx)]^{s}$ to which (192) can be reduced.

We can now use Green’s function to calculate the dynamics of the distribution with an abitrary initial condition:

[TABLE]

In the case of an empty immune system at $t=0$, $N_{C}(0)=0$, and constant thymic output $\theta_{k}$, we get:

[TABLE]

where $F_{1}$ is the Appell hypergeometric function of two variables, and where we made the change of variable $z=e^{-(t-t^{\prime})(\mu-\nu)}$ to carry out the integral in (194). This equation simplifies greatly for $\theta_{k}=\theta\delta_{k,1}$:

[TABLE]

XII.4.4 Solution to foward jump process with opposing drift and source

Consider the process described by (68) generalized to an arbitrary distribution of jumps:

[TABLE]

where absorbing boundary condition at $x=0$. where $J(y)dy$ is the rate (per unit time) of jumps of size between $y$ and $y+dy$. This equation can be solved exactly by considering the Laplace transform of $\rho$:

[TABLE]

which satisfies:

[TABLE]

where $\hat{J}(k)$ and $\hat{\rho}(k)$ are the Laplace transforms for $J$ and $\tilde{\theta}$. The steady state solution is:

[TABLE]

To get a physical solution, a pole at $k=0$ cannot exist, so that the total number of clones $\hat{\rho}(0)$ is well defined. This is satisfied for $\mu\rho(0)=\hat{\theta}(0)$, which is simply the condition that the rate of birth of new clones, $\hat{\theta}(0)$ be balanced by the rate of death, $\mu\rho(0)$. Then solution becomes

[TABLE]

To understand the behaviour of $\rho(x)$ at large $x$, we have to examine poles for negative $k$, $k=-\alpha$, which satisfy the condition:

[TABLE]

The left-hand side has derivative $\bar{\nu}\equiv\int dx\,J(x)x$ with respect to $k_{0}$, which is the average growth rate of a clone. To guarantee that all clones eventually go extinct, that number must be smaller than $\mu$. Therefore (202) has only one solution $\alpha>0$. At large $x$ the behaviour is dominated by that pole, $\hat{\rho}(k)\sim 1/(k+\alpha)$, yielding $\rho(x)\sim e^{-\alpha x}$ and:

[TABLE]

The total numbers of clones and cells read:

[TABLE]

In the specific case of fixed-size jump, $J(x)=sp\delta(x-\Delta x)$, we get back (70) and (71). Eq. 202 simplifies in the limit of many very small jumps, corresponding to $J_{\rm tot}\equiv\int_{0}^{\infty}dx\,J(x)\to\infty$ while keeping the average effect of jumps, $\bar{\nu}<\mu$, as well as its second moment, $\Gamma=\int_{0}^{\infty}dx\,J(x)x^{2}$, finite, meaning that the average jump size $\langle\Delta x\rangle=\nu/J_{\rm tot}$ goes to zero. Expanding (202) at small $x$, we get:

[TABLE]

which gives back (72) in the case of fixed-size jumps.

XII.4.5 The Yule process

The Yule-Simon Yule1925 ; Simon1955 process is another type of “neutral” process that is characterized by a distribution with power long tails, which in its tail is reminiscent of Zipf’s law (see section VII.2). The model was initially formulated to describe the distribution of biological genera. It was later adapted in network science and is known under the name of preferential attachement, or the rich-get-richer model. In general it is a realization of a birth-death process, where the birth coefficient is proportional to the number of individuals. In the model, new individuals are added at each time step. With probability $\alpha$, the new individual gets a mutation (or undergoes a speciation event if we think about species), and with probability $1-\alpha$ the new individual’s type is chosen among all the possible ancestors, i.e. proportionaly to the abundance of each type.

The probability that the rank of an individual is $k$ results in a power law distribution at large $C$:

[TABLE]

where $\Gamma(x)$ denotes a gamma function, and $\rho=1/(1-\alpha)$.

A precise derivation of (207) can be found in Simon1955 , but the power-law behaviour can be understood intuitively in a continous approximation. Assume that the rate of division is $\mu$, so that there are $C_{\rm tot}=e^{\mu t}$ individual at time $t$, starting from a single individual at time $t=0$. The rate of new emergence of new mutants is $\alpha\mu e^{\mu t}$, so that at time $t$ there are $n(t)=\int_{0}^{t}dt^{\prime}\alpha\mu e^{\mu t^{\prime}}=\alpha(e^{\mu t}-1)$ new types. The key point is that the abundance of a mutant is an exponentially increasing function of its age. A mutant that arose at $t^{\prime}$ has size $C=e^{\mu(1-\alpha)(t-t^{\prime})}$ at time $t$. $n(t^{\prime})$ can thus be viewed as the rank of the mutant, ordered by increasing frequency. The abundance of the mutant arising at time $t^{\prime}$, of rank $n=n(t^{\prime})$ is, at time $t$:

[TABLE]

where the approximation is only valid at large $n$. We recognize Zipf’s law (56), which implies a power law in the clone size distribution according to (55):

[TABLE]

XII.5 Inference

Here we introduce basic concepts and methods for constructing models and inferring their parameters from data. These approaches are becoming increasingly important as large datasets are being produced by high-throughput methods, from imaging to sequencing. Some of these methods have been used in only a few applications in computational immunology, but we anticipate that their usage will become widespread in future studies.

XII.5.1 Probabilistic inference, maximum likelihood and Bayesian statistics

Often, empirical data originates from a stochastic process. The first source of variability, which is present in almost any measurement, is experimental noise, which is random. In many cases related to biological systems, from gene expression to cell signaling and random recombination of DNA, the underlying processes are intrinsically stochastic, and so are the models that must describe them.

Let us assume that a model can be encapsulated into a probability distribution for observations $P(\bm{x}|\theta)$, where $\bm{x}$ is the empirical data, and $\theta$ the set of parameters defining the model. In many cases, the data can be decomposed into a set of independent datapoints, $\bm{x}=(x_{1},\ldots,x_{N})$, so that the probability distribution factorizes over these observations:

[TABLE]

This probability is called the likelihood (of the data given the model). Because of its often multiplicative form, it is common to consider instead its logarithm, the log-likelihood $\mathcal{L}(x|\theta)$, which is additive, or its normalized variant $\ell(x|\theta)=(1/N)\sum_{i=1}^{N}\ln P(x_{i}|\theta)$. One can think of $x_{i}$ as a sequence, for instance of an immune receptor or of a pathogen; a fluorescence signal coming from flow cytometry or single-cell microscopy imaging; an abundance of RNA transcripts for gene expression, or of sequencing reads of immune receptors when estimating the abundance of lymphocyte clones. The models may be statistical models for the occurence of a particular data point, or stochastic dynamical models of the process.

A popular way to infer the parameters $\theta=(\theta_{1},\ldots,\theta_{K})$ of a model is to find those that maximize the likehood of the data:

[TABLE]

This estimate of the parameters is called the maximum likelihood estimator, and it can be shown to be unbiased and optimal in the limit of large numbers of observation, in the sense that it saturates the Cramer-Rao bound. The Cramer-Rao bound gives a lower bound on the variance of unbiased estimators of a fixed and deterministic parameter. In practice, this implies that the fluctuations of the estimator around its true value, $\theta^{*}(\bm{x})-\theta$ are given by a multivariate Gaussian distribution of mean 0 and co-variance $I^{-1}/N$, with:

[TABLE]

is called the Fisher information matrix.

When one has information about the process and its parameters independently of the data, through a prior belief in the distribution of $\theta$, $P_{\rm prior}(\theta)$, this information can be combined with the data using Bayes rules:

[TABLE]

This expression defines the posterior distribution of parameters, given the data. For instance, some parameters may have been measured with some uncertainty, which can be modelled using a Gaussian distribution $P_{\rm prior}(\theta_{i})=(1/\sqrt{2\pi\sigma_{i}^{2}})\exp[-(\theta_{i}-\theta_{i,\rm meas})^{2}/\sigma_{i}^{2}]$, or there exists natural physiological range for them. A uniform prior, $P_{\rm prior}=$const, amounts to considering the likelihood alone. Also note that in the limit of large datasets, $N\gg 1$, the log-likelihood dominates over the logarithm of the prior, which becomes negligible. In that case, the posterior over possible parameters, $P(\theta|\bm{x})$, is given by a Gaussian distribution over $\theta$ with mean $\theta^{*}$ and covariance given by the same Fisher information matrix $I^{-1}/N$. This means that the Bayesian fluctuations, which correspond to our uncertainty about the true values of the parameters, match the fluctuations of the error we make when picking the maximum-likelihood estimate.

Using Eq. 213, one can either consider the full range of acceptable parameter values from the posterior, by e.g. sampling from it. Alternatively, one can take its maximum, $\theta^{*}=\operatorname*{arg\,max}_{\theta}P(\theta|\bm{x})$, called the maximum a posteriori estimator. In the limit of large $N$, that estimator is equivalent to the maximum likelihood estimator.

We now turn to two simple examples of inference problems for which the maximum likelihood estimator recovers an intuitive answer. First consider linear regression. Suppose we have data points $(x_{1},\ldots,x_{N})$ taken at times $(t_{1},\ldots,t_{N})$, which we want to fit with a linear model, $x=\theta_{1}+\theta_{2}t$. Implicitly, we must assume that the data is noisy, otherwise a fit would not be necessary. Assuming Gaussian noise, the stochastic model reads:

[TABLE]

where $\sigma$ is the amplitude of the noise. Examining this expression, one can see maximizing the likelihood corresponds to minimizing the mean squared error, $\sum_{i=1}^{N}(x_{i}-\theta_{1}-\theta_{2}t_{i})^{2}$, which is the standard way to do linear fit. This argument generalizes to any fit of a parametrized function, assuming Gaussian distributed noise.

The second example is estimating frequencies from a finite set of outcomes. For concreteness, let us assume that we want to describe the probability of finding a particular nucleotide, $x=$A,C,G,T, at a particular location in a genomic sequence. The parameters of the model are the frequencies of A, C, and G, denoted $\theta_{1},\theta_{2},\theta_{3}$, while the frequency of T is $1-\theta_{1}-\theta_{2}-\theta_{3}$. The probability of observing a certain set of observations $\bm{x}=(x_{1},\ldots,x_{N})$ is given by:

[TABLE]

where $n_{A}=\sum_{i=1}\delta_{x_{i},A}$ and similarly for the other nucleotides. Maximizing the likelihood with respect to the parameters yields the intuitive counting estimator: $\theta^{*}_{1}=n_{A}/(n_{A}+n_{C}+n_{G}+n_{T})$, and similarly for the other nucleotides. Note that this argument holds for sequences of nucleotides rather than single nucleotides, provided that the underlying probabilistic model assumes that the choices of nucletoides at each position along the sequence are independent of each other.

XII.5.2 Model selection

The maximum likelihood and Bayesian rules outlined above assume that the class of models, if not its parameters, are known, but this is not always the case. For instance, we may have competing hypotheses between different models, or we may consider models of increasing complexity. For instance, in the example of linear regression, we may want to know whether the data is not better explained by a quadratic model $x=\theta_{1}+\theta_{2}t+\theta_{3}t^{2}$. Comparing the likelihood directly always favors models with more parameters, which can in turn lead to overly complex models with the risk of overfitting the data. We need a way to compare models with different structures and numbers of parameters.

The general approach to compare two model A and B, described by probabilities $P(\bm{x}|\theta_{A},A)$ and $P(\bm{\sigma}|\theta_{B},B)$, is to evaluate the overall probability of the data given each model, by integrating the over the parameters:

[TABLE]

The relative probabilities of each model is then evaluated by Bayes’s rule:

[TABLE]

In the limit of large $N$, $P_{A}(\bm{x}|\theta_{A})$ (and the same for B) is very peaked around its maximum likelihood value, with Gaussian fluctuations around that peak given by the Fisher Information matrix $I_{A}$:

[TABLE]

By the saddle-point approximation we obtain:

[TABLE]

If we replace this expression into Eq. 216, and focus on two terms that do not depend on our prior assumptions, we see two terms: the maximized likelihood for each model $P_{A}(\bm{x}|\theta^{*}_{A})$, which quantifies the quality of the fit, and “parsimony” term $\sqrt{(2\pi/N)^{K_{A}}/\mathrm{det}(I_{A})}$, which quantifies the volume of parameters that are consistent with each model. As evident for this expression of the parsimony, the more parameters, the smaller that volume will be.

The ratio of probabilities in (216) between A and B is sometimes called Bayes’s factor, while the ratio of the parsimony terms is called Occam’s factor. If we expand the logarithm of (216) in $o(\ln(N))$ using (220), we obtain, up to a $-2$ factor, the so called Bayesian Information Criterion (BIC):

[TABLE]

The BIC is a popular score for comparing models, and can be intuitively understood as a correction to the likelihood by a term that penalizes large numbers of parameters. In the large $N$ limit, Eq. 218 tells us that the model with the smallest BIC is more likely to be correct. Another popular score is the Akaike information criterion (AIC), defined as $2K-2N\ell(\bm{x}|\theta^{*})$, which puts a smaller penalty on the number of parameters.

Applications of model selection criterion have been numerous in computational immunology, see e.g. Mccoy2015 ; Marcou2018c ; Cotari2013 ; Laydon2015 ; Kepler14 .

XII.5.3 Expectation-Maximization

Some models are better defined in terms of variables that are not accessible to the observer. If we call those hidden variables $\bm{h}$, and $\bm{x}$ the visible variables, then the likelihood of the data can be expressed as the sum over the hidden variables:

[TABLE]

Expectation-Maximization (EM) consists of maximizing $P(\bm{x}|\theta)$ iteratively. Starting with a guess $\theta^{t}$, one can write a pseudo log-likelihood:

[TABLE]

which is essentially an average of the log-likelihood of the full model (i.e. including the hidden variables), weighted by the posterior distribution of the hidden variable $h_{i}$ under the current model $\theta_{t}$. In an iteration step, one sets $\theta_{t+1}=\operatorname*{arg\,max}_{\theta}\mathcal{\tilde{L}}(\bm{x}|\theta)$. As this scheme is iterated, $\theta_{t}$ converges to a maximum of the true likelihood $P(\bm{x}|\theta)$.

Expectation-maximization has notably been used in the inference of generation model for antigen-specific lymphocyte receptors (TCR and BCR), in Refs. Murugan2012 ; Elhanati2014 ; Elhanati2015 ; Elhanati2016 .

XII.5.4 Hidden Markov models

Hidden Markov models (HMM) are models with a hidden variable $h$ which is the trajectory of a Markov process, $h=(h_{x})_{x=1,\ldots,L}$:

[TABLE]

The visible variable $x$, is also a vector of size $L$, whose elements $x_{j}$ are called ’emissions’ and are drawn according to the value of the hidden variable according to the emission probabilities $P(x|h)=\prod_{j=1}^{L}e_{j}(x_{j}|h_{j})$. A key feature of HMM is that high-dimensional sums over $h=(h_{1},\ldots,h_{L})$ can be performed recursively using technique which is equivalent to the transfer matrix technique in statistical physics. For instance, the forward algorithm allows one to calculate $P(x)$, which requires to sum over all possible trajectories of the hidden variable. This is done by defining $z_{j}(h_{j})=P(x_{1},\ldots,x_{j},h_{j})$ and using the recursive equation:

[TABLE]

The probability $P(x)=P(x_{1},\ldots,x_{L})$ is then obtained as $\sum_{h_{L}}z_{L}(h_{L})$. A similar trick allows one to calculate any marginal of the model, e.g. $P(x_{j})$ or correlation functions. Combining these recursions with EM defines a powerful algorithm for learning the parameters of the model, $\theta=\{w_{j}(\cdot|\cdot),e_{j}(\cdot|\cdot)\}$, called the Baum-Welch algorithm.

HMM have been applied to model BCR and TCR annotation to the germline in various software, e.g. SODA Munshaw2010 , Partis Ralph2016 , and repgenhmm Elhanati2016b .

XII.5.5 Information theory

The entropy of a probability distribution $P(x)$ is defined as:

[TABLE]

It quantifies the randomness of the distribution, and is maximum for a uniform distribution (and zero for perfectly peaked one). For this reason is is often used as a measure of diversity. Unlike all other diversity measures, it is additive, meaning that the entropy of a joint distribution of two independent variables $(P(x),Q(y))$ is given by $S[P]+S[Q]$.

The entropy is the building block for several information-theoretic measures. The Kullback-Leibler divergence or relative entropy, defined as

[TABLE]

is used a distance measure between probability distributions (although not a metric in the mathematical sense because of its asymmetry). Another popular measure is the mutual information between two variables $x$ and $y$:

[TABLE]

It quantifies how much the knowledge of $y$ reduces the entropy of $x$, and vice-versa since it is symmetric in $x$ and $y$. It is often used as a non-parametric measure of correlations between two variables.

XII.5.6 Maximum entropy models

When the underlying mechanisms that give rise to the data are not known, it can be useful to define phenomenological models based on the observables that are deemed important. A convenient way to do this is to infer maximum entropy models Presse2013a , which are probabilistic models $P(x)$ of maximum entropy $S[P]$ subject to the constraint that they agree with the data on a choice of key mean observables $\mathcal{O}_{a}(x)$: $\sum_{x}P(x)\mathcal{O}_{a}(x)=(1/N)\sum_{i=1}^{N}\mathcal{O}_{a}(x_{i})$. It can be shown that the distribution takes an exponential form:

[TABLE]

where $\theta_{a}$ are Lagrange multipliers that must be adjusted to satisfy the constraints, and $Z$ is a normalization constant. It turns out that maximizing the likelihood of a dataset $(x_{1},\ldots,x_{N})$ under Eq. 229 over the model parameters $(\theta_{a})$ is equivalent to satisfying the constraint over the mean observables. The inference procedure is computationally hard, and is usually performed using a combination of Monte Carlo sampling methods and gradient descent, or mean-field techniques Barton2016a . Often, the observables are taken to be marginals of variables of interest, as well as pairwise of higher-order correlations between them. The resulting models then fall into classes of inverse statistical physics models, such as disordered Ising or Potts models Nguyen2017 .

XII.5.7 Machine learning and Neural networks

Modern machine learning techniques are increasingly popular in computational immunology, and are anticipated to become even more popular in the near future. We refer the reader to Ref. Mehta2018 for a review of machine learning methods aimed at physicists.

The best established use of deep neural networks in immunology is for predicting peptide-MHC binding through a set of tools call netMHC Buus2003 ; Andreatta2015 . More traditional machine learning approaches have been used to compare and characterize immune repertoires Thomas2014b ; Cinelli2017 .

Current efforts aim at predicting TCR-antigen binding using machine learning techniques Jurtz2018 ; Sidhom2018 ; Jokinen2019 , although the amount of data necessary to obtain truly predictive models is probably still insufficient at this stage.

XII.6 High-throughput repertoire sequencing

The variable region of an antigen receptor chain is about 400 bp long. The sequencing challenge is to capture the whole region in one read. Modern methods of high-throughput genomic or metagenomic sequencing use shotgun sequencing which breaks the genome of interest into short fragments and after sequencing pastes the reads together by putting together overlapping reads with the help of a reference guiding template. Since Repertoire Sequencing (RepSeq) focusses on highly variable regions which includes non-templated insertions and deletions, the sequence must be acquired in one read. Read lengths required for TCR sequences are typically shorter ($\sim 100-150$bp) since short fragments of V and J genes are enough to distinguish different V genes with the help of known genomic templates. However, since B-cells carry a large number of hypermutations in the V gene, sequencing of longer fragments ($\sim 600$bp) that encompass explicitly most of the V gene is necessary.

Different methods have been developed for sequencing DNA and mRNA. In general, the protocols start with isolating the mRNA or DNA of the cells of interest (BCR, TCR, in subsets of cells of interest sorted using FACS). The mRNA product is then reverse transcribed onto cDNA and barcoded, sometimes with unique molecular identifier (UMI), and amplified by Polymerase Chain Reaction (PCR) before sequencing. Different techniques exist for PCR amplification that either add primer specific sites or multiple primers. Rapid amplification of cDNA ends (RACE) is one of the most common techniques, which is based on adding a linker with a primer binding site to a conserved region, on the 5’ end of the V gene where no constant template is available (unlike on the 3’ end). In DNA sequencing multiple primers that target specific regions of the DNA, introducing a primer specific amplification bias that can be controlled for by spiking in known sequences.

Protocols based on sequencing directly DNA have the advantage that the experiments is free of mRNA expression bias. However, currently mRNA based technologies are able to reliably report sequence counts thanks to the use of UMI barcoding techniques that are still being developed for DNA sequencing. The difficulty with using UMIs in the DNA protocol lies in introducing the barcode before the original sequence is amplified. In the mRNA protocole the barcode is introduced during the initial reverse transcription of mRNA into cDNA for every RNA molecule. The barcoded sequence is then purified and the whole product is amplified by PCR. For DNA an additional ligation step is needed, in which the DNA is cut close to the region of interest (keeping in mind that the sequenced lengths is short). Ligation is not an efficient reaction. This approach is simple when many gene copies are available in a sample, but this is not the case in immune repertoire sequencing. Alternatively, the barcode could be added in the first PCR cycle. The product then needs to be well purified making show no barcode carrying primers are left, at the same time making sure that the barcoded sequences is not lost due to its low fraction compared to non-barcoded sequences, which remains technically tricky.

The above procedures sequence the chains of TCR or BCR repertoire in bulk, leaving no possibility to figure out the pairing of alpha and beta, or light and heavy chains in the same cell. Naturally, this limits the discussions of repertoire diversity, and more importantly attempts to link sequence to a functional phenotype. Recently these limitations have been overcome by the development of single cell repseq sequencing, where the PCR reaction is performed either in wells Dekosky2014 or in droplets Klein2015 ; Mcdaniel2016 ; Grigaityte2017 . Alternatively, a computational approach that exploits the statistics of co-occurence of the rare events partitioning the chains separately into wells is used to identify pairs of TCR sequences that come from the same cell from bulk DNA sequencing experiments Howie2015 . The results of these analysis show that the pairing of alpha-beta chains is largely independent for TCR Dupic2018 ; Grigaityte2017 , so the results of bulk analysis for diversity hold. However the affinity and functional properties of the BCR and TCR need to be estimated at the level of the whole receptor.

Lastly, novel high-throughput single-cell barcoded mRNA technologies Zheng2017 are promising to change the field of immune repertoire sequencing, providing paired chain reads in large numbers. These technologies are very similar to those used in 5’RACE with UMI. Droplet based platforms allow in principle for large numbers of cells (up to 80 000 per chip) to be analyzed at the same time, although these numbers are still much smaller than the $10^{6}$ cells in bulk experiments. Given that clones have very low frequencies in bulk experiments (possibly down to a single cell out of $10^{11}$), sampling of the order of $10^{4}$ does not guarantee reproducible experiments. Additionally the cost remains high, at about $1$ USD per cell. At the moment of writing, the first analyses based on this technology are still underway.

XIII Acknowledgments

The authors thank all past and current collaborators for the many useful discussions. TM and AMW thanks Meriem Bensouda Koraichi, Barbara Bravi, Victor Chardès, Thomas Dupic, Cosimo Lupo, Carlos Olivares, Jacopo Marchi, Maria Ruiz Ortega and Natanael Spisak for their comments on the manuscript and Mikhail Pogorelyy and Anastasia Minervina for discussions. GAB thanks Sooraj Achar, Emanuel Salazar Cavazos and Van Truong for their comments on the manuscript. This work was partially supported by the European Research Council Consolidator Grant n. 724208 [TM & AMW] and the intramural research program of the National Cancer Institute [GAB]. We also thank the KITP, Les Houches EdP and IES Cargèse institutes for hospitality during workshops on the topics reported in this review.