Fluctuating fitness shapes the clone size distribution of immune repertoires

The adaptive immune system relies on the diversity of receptors expressed on the surface of B and T-cells to protect the organism from a vast amount of pathogenic threats. The proliferation and degradation dynamics of different cell types (B cells, T cells, naive, memory) is governed by a variety of antigenic and environmental signals, yet the observed clone sizes follow a universal power law distribution. Guided by this reproducibility we propose effective models of somatic evolution where cell fate depends on an effective fitness. This fitness is determined by growth factors acting either on clones of cells with the same receptor responding to specific antigens, or directly on single cells with no regards for clones. We identify fluctuations in the fitness acting specifically on clones as the essential ingredient leading to the observed distributions. Combining our models with experiments we characterize the scale of fluctuations in antigenic environments and we provide tools to identify the relevant growth signals in different tissues and organisms. Our results generalize to any evolving population in a fluctuating environment.

The adaptive immune system relies on the diversity of receptors expressed on the surface of B and T-cells to protect the organism from a vast amount of pathogenic threats. The proliferation and degradation dynamics of different cell types (B cells, T cells, naive, memory) is governed by a variety of antigenic and environmental signals, yet the observed clone sizes follow a universal power law distribution. Guided by this reproducibility we propose effective models of somatic evolution where cell fate depends on an effective fitness. This fitness is determined by growth factors acting either on clones of cells with the same receptor responding to specific antigens, or directly on single cells with no regards for clones. We identify fluctuations in the fitness acting specifically on clones as the essential ingredient leading to the observed distributions. Combining our models with experiments we characterize the scale of fluctuations in antigenic environments and we provide tools to identify the relevant growth signals in different tissues and organisms. Our results generalize to any evolving population in a fluctuating environment.

I. INTRODUCTION
Antigen-specific receptors expressed on the membrane of B and T cells (BCRs and TCRs ) recognize pathogens and initiate the adaptive immune response [1]. An efficient response relies on the large diversity of receptors that protect us against the many pathogens in the environment. The adaptive repertoire diversity is maintained from a source of newly generated cells, each expressing a unique receptor. These progenitor cells later divide or die, and their offspring make up clones of cells that share a common receptor. The sizes of clones vary, as they depend on the particular history of cell divisions and deaths in the clone. The clone size distribution thus bears signatures of the challenges faced by the adaptive system. Understanding the form of the clone size distribution in healthy individuals is an important step in the characterization of the antigenic recognition process and the functioning of the adaptive immune system. It also presents an important starting point for describing statistical deviations seen in individuals with compromised immune responses.
High throughput sequencing experiments in different cell types and species have allowed for the quantification of clone sizes and their distributions [2][3][4][5][6][7][8][9]. Generically, clone size distributions exhibit heavy tails [2,[9][10][11]-the number of clones of a given size is well described by a power law over a wide range of clone sizes (see Fig. 1A-B). Such universal behavior for a variety of cells types in different species (B cells, T cells, naive cells, effector cells in fish, mice, and humans) suggests common underlying rules for growth, death and homeostatic control of adaptive repertoires. Previous work has described the specifics of immune dynamics for a certain cell type [12,13], or a certain signaling pathway, using detailed mechanistic models [14][15][16]. It remains unclear, however, what essential features of these dynamics may lead to the observed power-law distributions, and what are the key biological parameters of the repertoire dynamics that govern its behavior.
The wide range and types of interactions that influence a B or T cell fate happen in a complex, dynamical environment with inhomogeneous spatial distributions. They are difficult to measure in vivo, making their quantitative characterization elusive. To overcome this, we describe the effective interaction between the immune cells and their environment as a stochastic process governed by only a few relevant parameters. This effective description is consistent with the idea, common in physics, that the apparent universality of the observed distribution is underlied by broad model universality classes. Since the heavy tail behavior exists in both B and T cells, we consider general properties that are common to both cell types, and ignore hypermutations for simplicity.
All cells proliferate and die depending on the strength of antigenic and cytokine signals they receive from the environment, which together determine their net growth rate. This effective fitness that fluctuates in time is central to our description. We find that its general properties determine the form of the clone size distribution. Two broad classes of models are distinguished, according to whether these fitness fluctuations are clone specific (mediated by their specific BCR or TCR) or cell specific (mediated by phenotypic fluctuations such as the number of cytokine receptors). We identify the models that are compatible with the experimentally observed distributions of clone sizes. These distributions do not depend on the detailed mechanisms of cell signaling and growth, but rather emerge as a result of self-organisation, with no need for fine-tuned interactions.

A. Clone dynamics in a fluctuating antigenic landscape
The fate of the cells of the adaptive immune system depends on a variety of clone-specific stimulations. The recognition of pathogens triggers large events of fast clone proliferation followed by a relative decay, with some cells being stored as memory cells to fend off future infections. Naive cells, which have not yet recognized an antigen, do not usually undergo such extreme events of proliferation and death, but their survival relies on short binding events (called "tickling") to antigens that are natural to the organism (self-proteins) [17,18]. Because receptors are conserved throughout the whole clone (ignoring B cell hypermutations), clones that are better at recognizing self-antigens and pathogens will on average grow to larger populations than bad binders. By analogy to Darwinian evolution, they are "fitter" in their local, timevarying environment [12,19,20].
We denote by a j (t) the overall concentration of an antigen j as a function of time. We assume that after its introduction at a random time t j , this concentration decays exponentially with a characteristic lifetime of antigens λ −1 , a j (t) = a j,0 e −λ(t−tj ) as pathogens are cleared out of the organism, either passively or through the action of the immune response. Lymphocyte receptors are specific to certain antigens, but this specificity is degenerate, a phenomenon refered to as cross-reactivity or polyspecificity. The extend to which a lymphocyte expressing receptor i interacts with antigen j is encoded in the cross-reactivity function K ij , which is zero if i and j do not interact, or a positive number drawn from a distribution to be specified, if they do. In general, interactions between lymphocytes and antigens effectively promote growth and suppress cell death, but for simplicity we can assume that the effect is restricted to the division rate. In a linear approximation, this influence is proportional to j K ij a j (t), i.e. the combined effect of all antigens j for which clone i is specific. This leads to the following dynamics for the evolution of the size C i of clone i: where ν and µ are the basal division and death rates, and where Bξ i (t) is a birth-death noise of intensity B 2 = (ν + j K ij a j (t) + µ)C i , with ξ i (t) a unit Gaussian white noise (see Appendix A for details about birth death noise). New clones, with a small typical initial size C 0 , are constantly produced and released into the periphery with rate s C . For example, a number of the order of s C = 10 8 new T-cells are output by the thymus daily in humans [21]. Since the total number of T cells is of the order of 10 11 , this means that cells die with an average rate of 10 −3 days −1 in homeostatic conditions [21]. Because the probability of rearranging the exact same receptor independently is very low (< 10 −10 ) [22], we assume that each new clone is unique and comes with its own set of cross-reactivity coefficients K ij . Assuming a rate s A of new antigens, the average net growth rate in Eq. 1 is f 0 = ν + a j,0 K s A λ −1 − µ < 0, and the stationary number of clones should fluctuate around N C ≈ s C |f 0 | −1 clones. This is just an average however, and treating each clone independently may lead to large variations in the total number of cells (i.e. the sum of sizes of all clones). To maintain a constant population size, clones compete with each other for specific resources (pathogens or selfantigens) and homeostatic control can be maintained by a global resource such as Interleukin 7 or Interleukin 2.
Here we do not model this homeostatic control explicitly, but instead assume that the division and death rates ν, µ are tuned to achieve a given repertoire size. We verified that adding an explicit homeostatic control did not affect our results (see Fig. S2 and Appendix B).
We simulated the dynamics of a population of clones interacting with a large population of antigens. Each antigen interacts with each present clone with probability p = 10 −7 , and with strength K ij drawn from a Gaussian distribution of mean 1 and variance 1 (truncated to positive values). A typical trajectory of the antigenic stimulation undergone by a given clone, j K ij a j , is shown in Fig. 1E (green curve), and shows how clone growth tracks the variations of the antigenic environment. When the stimulation is particularly strong, the model recapitulates the typical behaviour experimentally observed at the population level following a pathogenic invasion [23,24], as illustrated in Fig. 1D: the population of a clone explodes (red curve), driving the growth of the total population (blue curve), while taking over a large fraction of the carrying capacity of the system, and then decays back as the infection is cleared.
On average, the effects of division and death almost balance each other, with a slight bias towards death because of the turnover imposed by thymic or bone marrow output. However, at a given time, a clone that has high affinity for several present antigens will undergo a transient but rapid growth, while most other clones will decay slowly towards exctinction. In other words, locally in time, the antigenic environment creates a unique "fitness" for each clone. Since growth is exponential in time, these differential fitnesses can lead to very large differences in clone sizes, even if variability in antigen concentrations or affinities are nominally small. We thus expect to observe large tails in the distrubution of clone size. Fig. 2A shows the cumulative probability distribution function (CDF) of clone sizes obtained at steady state (blue curve) showing a clear power-law behaviour for large clones, spanning several decades.
The exponent of the power-law is independent of the introduction size of clones (see inset of Fig. 2A), and the specifics of the randomness in the environment (exponential decay, random number of partners, random interaction strength) as long as its first and second moment are kept fixed (See Fig. S3 and Appendix C).

B. Simplified models and the origin of the power law
To understand the power-law behavior observed in the simulations, and its robustness to various parameters and A. B cell zebrafish experimental cumulative clone size distribution for fourteen fish as a function of the fraction of the population occupied by that clone from data in Weinstein et al. [2]. B. Clone size distribution for murine T-cells from Zarnitsyna et al. [11] (data plotted as presented in original paper). C. The dynamics of adaptive immune cells include specific interactions with antigens that promote division and prevent cell death. New cells are introduced from the thymus or bone marrow with novel, unique receptors. Division, death and thymic or bone marrow output on average balance each other to create a steady state population. D-E. Example trajectories from simulations of the immune cell population dynamics in Eq. 1. The total number of cells (D) shows large variations after an exceptional event of a large pathogenic invasion. One or a few cells that react to that specific antigen grow up to a macroscopic portion of the total population, and then decrease back to normal sizes after the invasion. A typical clone size trajectory along with its pathogenic stimulation j Kijaj(t) shows the coupling between clone growth and variations of the antigenic environment (E). Parameters used: sC = 2000 day −1 , C0 = 2, sA = 1.96 · 10 7 day −1 , aj,0 = a0 = 1, λ = 2 day −1 , p = 10 −7 , ν = 0.98 day −1 , µ = 1.18 day −1 .
sources of stochasticity, we decompose the overall fitness of a clone at a given time (its instantaneous growth rate) into a constant, clone-independent part equal to its average f 0 < 0, and a clone-specific fluctuating part of zero mean, denoted by f i (t). This leads to rewriting Eq. 1 as: with B 2 ≈ (|f 0 | + 2µ)C i . The function f i (t) encodes the fluctuations of the environment as experienced by clone i. Because antigens can be recognized by several receptors, these fluctuations may be correlated between clones. Assuming that these correlations are weak, f i (t)f j (t ) ≈ 0, amounts to treating each clone independently of each other, and thus to reducing the problem to the single clone level. The stochastic process giving rise to f i (t) is a sum of Poissondistributed exponentially decaying spikes. This process is not easily amenable to analytical treatment, but we can replace it with a simpler stochastic process with the same temporal autocorrelation function. This autocorrelation is given by f i (t)f i (t ) = A 2 e −λ|t−t | , with the antigenic noise strength A 2 = s A pa 2 0 K 2 λ −1 , and where we recall that λ −1 is the characteristic lifetime of antigens. The simplest process with the same autocorrelation function is given by an overdamped spring in a thermal bath, or Ornstein-Uhlenbeck process, with η i (t) a Gaussian white noise of intensity 1 and γ = A √ λ quantifies the strength of variability of the antigenic environment (see Appendix D). This is also the process of maximum entropy or caliber [25] with that autocorrelation function (see Appendix E and [26]).
The effect of the birth death noise Bξ i (t) is negligible when compared to the fitness variations for large clones and it has no effect on the tail (see Fig. S5 and Appendix F). It can thus be ignored when looking at the tail of the distribution and its power law exponent, but it will play an important role for defining the range over which the power law is satisfied.
The population dynamics described by Eqs. 2 and 3 can be reformulated in terms of a Fokker-Planck equation for the joint abundance ρ of clones of a given log-size x = log C and a given fitness f : where the source term s(x, f ) describes new clones arriving at rate s C with size C 0 and normally distributed fitnesses of variance f 2 = γ 2 /λ. This Fokker-Planck equation can be solved numerically with finite element methods with an absorbing boundary condition at x = 0 to account for clone extinction. The solution, represented by the black curve in Fig. 2A, matches closely that of the full simulated population dynamics (in blue). The power-law behaviour is apparent above a transition point that depends on the distribution of introduction sizes of new clones and the parameters of the model (see below). Intuitively, the microscopic details of the noise are not expected to matter when considering long time scales, as a consequence of the central limit theorem. However, the long tails of the distribution of clone sizes involve rare events and belong to the regime of large deviations, for which these microscopic details may be important. Therefore, the agreement between the process described by the overdamped spring and the exponentially decaying, Poisson distributed antigens is not guaranteed, and in fact does not hold in all parameter regimes (see Fig. S8).
We can further simplify the properties of the noise by assuming that its autocorrelation time is small compared to other timescales. This leads to taking the limit γ, λ → ∞ while keeping their ratio constant σ = γ/λ constant, so that f i (t) is just a Gaussian white noise with Fig. S4). The corresponding Fokker-Planck equation now reads ). This equation can be solved analytically at steady state, and the resulting clone size distribution is, for C > C 0 : with The full solution, represented in Fig. 2A in red, captures well the long-tail behaviour of the clone size distribution despite ignoring the temporal correlations of the noise, and approaches the solution of the colored-noise model (Eq. 3) as λ, γ → ∞, as expected (see Fig. 2A). The power law behaviour and its exponent depend on the noise intensity, but are otherwise insensitive to the precise details of the microscopic noise, including its temporal properties. Fat tails (small α) are expected when the average cell lifetime is long (small |f 0 |) and when the antigenic noise is high (large σ or A). The explicit expression for the exponent of the power law 1 + α as a function of the biological parameters can be used to infer the antigenic noise strength A 2 directly from data. An inverse lifetime of |f 0 | ≈ 10 −3 days −1 of a typical T-cell can be estimated as the ratio of thymic output to the total population of lymphocytes in the body [27][28][29]. The characteristic lifetime of antigens λ −1 is harder to estimate, as it corresponds to the turnover time of the antigens that the body is exposed to, but is probably of the order of days or a few weeks, λ ≈ 0.1 day −1 . We estimated α = 1 ± 0.2 from the zebrafish data of Fig. 1A [2, 10] using canonical methods of power-law exponent extraction [30] (see Appendix G for details), and also found a similar value in human T cells [31]. The resulting estimate, A = 10 −2 day −1 , is rather striking, as it implies that fluctuations in the net clone growth rate, A, are much larger than its average f 0 .
While the distribution always exhibits a power law for large clones, this behavior does not extend to clones of arbitrarily small sizes, where the details of the noise  and how new clones are introduced matter. We define a power-law cut-off C * as the smallest clone size for which the cumulative distribution function (CDF) differs from its best power-law fit by less than 10%. Using numerical solutions to the Fokker-Planck equation associated to the colored-noise model, we can draw a map of C * as a function of the parameters of the system. In Fig. 2B-C we show how C * varies as a function of the introduction size for different values of the dimensionless parameter related to the effective strength of antigen fluctuations relative to their characteristic lifetime at fixed power law exponents. In principle one can use this dependency to infer effective parameters from data. In practice, when dealing with data it is more convenient to consider the value of the cumulative distribution at C * , rather than C * itself. For example, fixing C 0 = 4 and fitting the curve of Fig. 1A with our simplified model using λ as an adjustable parameter, we obtain λ ≈ 0.14 day −1 (see Appendix G), which corresponds to a characteristic lifetime of antigens of around a week. Although this estimate must be taken with care, because of possible PCR amplification biases plaguing the small clone size end of the distribution, the procedure described here can be applied generally to any future repertoire sequencing dataset for which reliable sequence counts are available.

C. A model of fluctuating phenotypic fitness
So far, we have assumed that fitness fluctuations are identical for all members of a same clone. However, the division and death of lymphocytes do not only depend on signaling through their TCR or BCR. For example, cytokines are also growth inducers and homeostatic agents [32,33], and the ability to bind to cytokines depends on single-cell properties such as the number of cytokine receptors on the membrane of a given cell, independent of their BCR or TCR receptor. Other stochastic single-cell factors may affect cell division and death. These signals and factors are cell specific, as opposed to the clone specific properties related to BCR or TCR binding. Together, they define a global phenotypic state of the cell that determines its time-varying "fitness," independent of the clone and its T-cell or B-cell receptor. This does not mean that these phenotypic fitness fluctuations are independent across the cells belonging to the same clone. Cells within a clone share a common ancestry, and may have inherited some phenotypic properties of their common ancestors, making their fitnesses effectively correlated with each other. However, this phenotypic memory gets lost over time, unlike fitness effects mediated by antigen-specific receptors.
We account for these phenotypic fitness fluctuations by a function f c (t) quantifying how much the fitness of an individual cell c differs from the average fitness f 0 . This fitness difference is assumed to be partially heritable, which we model by: where λ −1 is the heritability, or the typical time over which the fitness-determining trait is inherited, γ c quantifies the variability of the fitness trait, and η c (t) is a cell-specific Gaussian white noise of power 1. Despite its formal equivalence with Eq. 3, it is important to note that here the fitness dynamics occurs at the level of the single cell (and its offspring) instead of the entire clone. The dynamics of the fitness f i (t) of a given clone i can be approximated from Eq. 7 by averaging the fitnesses f c (t) of cells in that clone, yielding: where η i (t) and ξ i (t) are clone-specific white noise of intensity 1, and ν and µ are the average birth and death rates, respectively, so that f 0 = ν − µ (details in Appendix I). The difference with Eq. 3 is the 1/ C i (t) prefactor in the fitness noise η i (t), which stems from the averaging of that noise over all cells in the clone, by virtue of the law of large numbers. Because of this prefactor, the fitness noise is now of the same order of magnitude as the birth-death noise, which must now be fully taken into account. Taking Eq. 8 and Eq. 9 at the population level gives a Fokker-Planck equation with a source term accounting for the import of new clones. We verify the numerical steady state Fokker-Planck solution against Gillespie simulations (Fig S6, see Appendix H for details). Fig. 3A-B show the distribution of clone sizes for different values of the phenotypic relaxation rate λ and environment amplitude γ c . These distributions vary from a sharp exponential drop in the case of low heritability (large λ) to heavier tails in the case of long conserved cell states (small λ). To quantify the extend to which these distributions can be described as heavytailed, we fit them to a power law with exponential cutoff, ρ(C) ∝ C −1−α e −C/Cm , where C m is the value below which the distribution could be interpreted as an (imperfect) power law. Fig. 3C shows a strong dependency of this cut-off with the phenotypic memory λ −1 . The longer the phenotypic memory λ −1 , the more clone-specific the fitness looks like, and the more the distribution can be mistaken for a power law in a finite-size experimental distribution. Larger birth-death noise also extends the range of validity of the power-law. As a result, and despite the absence of a true power-law behaviour, these models of fluctuating phenotypic fitnesses cannot be discarded based on current experimental data.
The model can be solved exactly at the two extremes of the heritability parameter λ. In the limit of infinite heritability (λ → 0) the system is governed by selective sweeps. The clone with the largest fitness completely dominates the population, until it is replaced by a better one, giving rise to a trivial clone-size distribution. In the opposite limit, when heritability goes to 0 (λ → +∞), the Fokker-Planck equation can be solved analytically (see Appendices I and J), yielding an exact power-law with exponential cutoff, ρ(C) ∝ The numerical solution of Fig. 3B is close to this limit. Note that even with a negligible exponential cutoff, the predicted α < 0 contradicts experimental observations.

III. DISCUSSION
Clone size distributions from high throughput data have been available for a few years and hold great promise for understanding the processes that shape the composition of the repertoire. Yet their interpretation requires models. The fate of lymphocytes making up immune repertoires is governed by a wide variety of mechanisms and pathways that vary across cell types and species. Previous population dynamics approaches to repertoire evolution have taken great care in precisely modeling these processes for each compartment of the population, through the various mechanisms by which cells grow, die and change phenotype [12,13,34]. However, one of the most striking properties of repertoire statistics revealed by high-throughput sequencing-the observation of power laws in clone size distribution-hold true for various species (human, mice, zebrafish), cell type (B and T cells) and subsets (naive and memory, CD4 and CD8), and seems to be insensitive to these context-dependent details. These observations call for a stochastic description of the various properties affecting cell fate, and their ubiquity points to some universal features of the dynamics.
The model introduced in this paper describes the stochastic nature of the immune dynamics with a minimal number of parameters, easing the numerical and analytical treatment and helping interpret the different regimes. These parameters are effective in the sense that they integrate different levels of signaling, pathways, and mechanisms. We assumed that they are general enough that different cell types (B and T cells) or subsets can be described by the same dynamical equations. Most of these parameters have natural interpretations: the typical lifespan of an infection or antigen stimulation is given by λ −1 , and the effect of these stimulations by A. For memory cells, λ −1 corresponds to the duration of an infection. For naive cells, the interpretation is less clear, as non-pathogenic antigens such as self-antigens may not be as transient as pathogens; in this case λ −1 could correspond to the characteristic timescale of changes in a clone's micro-environment. Likewise, in the context of cell-specific noise, the parameters λ −1 and γ c correspond to the correlation time and variability of the cell phenotype. All these parameters may vary greatly across cell types and species.
We showed that the relevant sources of stochasticity for the shape of the clone-size distributions fall into two main categories, depending on how cell fate is affected by the environment. Either the stochastic elements of clone growth act in a clone-specific way, through their receptor (BCR or TCR), leading to power-law distributions with exponent ≥ 1, or in a cell-specific way, e.g. through their variable level of sensitivity to cytokines (and more generally through any phenotypic trait affecting cell fitness), leading to exponentially decaying distributions with a power-law prefactor. These two types of signals (clone specific and cell specific) are important for the somatic evolution of the immune system [20,32,33,[35][36][37] and our analysis shows that the shape of the clone size distribution is informative of their relative importance to the repertoire dynamics. It provides a first theoretical setting and an initial systematic classification for modeling immune repertoire dynamics. Our method applied to high-throughput sequencing data can be used to quantify how much each type of signal contributes to the overall dynamics, and what is the driving force for the different cell subsets. For example, although it is reasonable to speculate that clone-specific signals should dominate for memory cells (through antigen recognition), and cell-specific selection for naive cells (through cytokinemediated homeostatic division), the relative importance of these signals for both cell types is yet to be precisely quantified, and may vary across species. A clear power law over several decades would strongly hint at dynamics dominated by interactions with antigens, while a faster decaying distribution would favor a scenario where individual cell fitness fluctuations dominate. Applying these methods to data from memory cells can give orders of magnitude for the division and half-life of memory lymphocytes, as well as the typical number of cells C 0 from a clone that are stored as memory following an infection.
The application of our method to data from the first immune repertoire survey (B cell receptors in zebrafish [2]) suggests that clone-specific noise dominates in that case, allowing us to infer a relation between the dynamical parameters of the model from the observed powerlaw exponent ≈ 2. However, there are a few issues with applying our method directly to data in the current state of the experiments. First, the counts (i.e. how many cells have the same receptor sequence and belong to the same clone) from many high-throughput repertoire sequencing experiments are imperfect because of PCR bias and sampling problems. New methods using single-molecule barcoding have been developed for RNA sequencing [8,38,39], but they do not solve the problem entirely, as the number of expressed mRNA molecules may not faithfully represent the cell numbers because of possible expression bias. In addition, most studies (with the exception of [40]) have been sequencing only one of the two chains of lymphocyte receptors, which is insufficient to determine clone identity unambiguously. As methods improve, however, our model can be applied to future data to distinguish different sources of fitness stochasticitiy and to put reliable constraints on biological parameters.
Thanks to its generality, our model is also relevant beyond its immunological context, and follows previous attempts to explain power laws in other fields [41][42][43]. The dynamics described here corresponds to a generalization of the neutral model of population genetics [44] where thymic or bone marrow outputs are now reintepreted as new mutations or speciations, and where we have added a genotypic or phenotypic fitness noise (receptor or cell-specific noise, respectively). It was recently shown that such genotypic fitness noise strongly affects the fixation probability and time in a population of two alleles [45,46]. Our main result (Eq. 6) shows how fitness noise can cause the clone-size distribution (called frequency spectrum in the context of population genetics) to follow a power law with an arbitrary exponent > 1 in a population of fixed size, while the classical neutral model gives a power law of exponent 1 with an exponential cut off (as shown in our exact solution with γ c = 0). Our results can be used to explain complex allele frequency spectra using fluctuating fitness landscapes.
Acknowledgements. This work was supported in part by grant ERCStG n. 306312. We compare results from a full Gillespie simulation (blue crosses) of a system with only birth-death dynamics with analytical prediction for a discrete system (black crosses, Eq. A3) and a continuous system (red curve, Eq. A12). The prediction with discrete variables is more accurate for small clones but the behaviour of all systems is the same for large populations. The parameters are µ = 1.45 day −1 , ν = 1.5 day −1 , C0 = 2 and we introduce 2000 new clones per day.
Appendix A: Simple birth-death process with no fitness fluctuations, and its continuous limit In this Appendix we derive the steady-state clone size distribution for a system that does not experience any environmental stimulation or noise, but is governed by a birth death process. We will show that the small number fluctuations arising from the discrete nature of birth and death are not sufficient to explain the observed distributions. We also show that our choice of a continuous birth death process is equivalent to its discrete version.
The multiplicative birth-death process corresponds to the following discrete dynamics: where µ is the division rate, ν the death rate. We assume that the population of cells of size n is maintained out of equilibrium by a source of new cells. The steady state solution for cell numbers above the value of the source satisfies detailed balance P (n)µn = P (n + 1)ν(n + 1) (A2) and, assuming the death rate is larger than the birth rate, takes the form The continuous counterpart of this discrete stochastic process corresponds to the following linear-noise approximation: where ξ i (t)ξ i (t ) = δ(t − t ) and f 0 = µ − ν < 0 (and we use theÎto convention ). In terms of x = log C the Langevin equation is and the corresponding Fokker-Planck equation reads where s(x) is the distribution of sizes of newly arriving clones. At steady state, we find where K is an integration constant. Defining To ensure convergence we set K = s C /(|f 0 |C m ) and the steady solution of the Fokker-Planck equation is This result is exactly equivalent to that of Eq. A3 when ν−µ = |f 0 | µ, ν. The accuracy of the approximation is verified in Fig. S1. Even for very large exponential cutoff values, C m , the apparent exponent is α = 0, corresponding to a flat cumulative distribution. This distribution is inconsistent with experiments, regardless of sequencing depth and we conclude that pure birth-death noise is not sufficient to explain the observed distributions.

Appendix B: Effects of explicit global homeostasis
In the simulations of clone dynamics in a fluctuating environment presented in the "Clone dynamics in a fluctuating antigenic landscape" Results section of the main text, we did not explicitly include a homeostatic control term, but tuned the division and death rates to achieve a given repertoire size. Here we add an explicit homeostatic term to the growth and degradation terms in the Langevin simulations described by Eq. 1 of the main text where N is a carrying capacity, h is the homeostatic constant multiplicator and r is the exponent of homeostatic response that described the sharpness of the response when approaching then carrying capacity limit. Comparing in Fig. S2 the resulting clone size distribution obtained with the explicit homeostatic term to the distribution from the simulations in the main text, we see that the explicit homeostatic term does not have an effect on the form of the distribution. It does have an effect on the trajectory of certain clones, and in particular on the response of the system to a very large invasion, making it an important feature of the dynamics of the immune system. However, as shown by the results in Fig. S2 its net effect on the clone size distribution can be taken into account by tuning division and death. When considering specific trajectories in the mean field approximation homeostatic control will add a systematic negative drift to the clonal population and can be accounted for by an additional contribution to f 0 .

Appendix C: Details of noise partition do not influence the clone size distribution function
In the simulation of the dynamics of receptors experiencing a clone-specific fitness presented in the "Clone dynamics in a fluctuating antigenic landscape" Results section of the main text we distributed the noise between the different random distributions: the poisson distributed number of new antigens (s A ), the variance of the initial concentrations (a j,0 ) and the variance of the binding probability (the values of K ij ). We made specific choices for this reparation by picking specific parameters of the random processes. Here we show that these specific choices of repartitioning the contributions to the noise do not influence the clone size distributions. Fig. S3 compares clone size distributions obtained with different values of the poisson distributed number of newly arriving antigen N a and the variance of the Gaussian distributed binding probabilities K ij , reproducing the same distributions in both cases. In the "Simplified models and the origin of the power law" Results section of the main text we make a series of approximations to effectively describe the dynamics of immune cells: we first approximate the antigenic environment by a random process with time correlated (colored) noise and we later neglect these temporal correlations. In this section and Appendix F we give the details that lead to the specific forms of the effective equations. In this Appendix we derive the Fokker-Planck equations for the time correlated noise model. In Appendix F we will consider the limit of an infinitely quickly changing environment.
The Langevin equations describing the dynamics of cells experiencing clone specific fitness fluctuations with a finite correlation time are where ξ i (t)ξ i (t ) = δ(t − t ) represents birth death noise in the linear-noise approximation (with theÎto convention) and η i (t)η i (t ) = δ(t − t ) is the noise of anti- The parameters were taken to be (as in Fig. 1) sC = 2000 day −1 , C0 = 2, day −1 , aj,0 = a0 = 1, λ = 2 day −1 , p = 10 −7 , ν = 0.98 day −1 , µ = 1.18 day −1 . For the red curve the variance of the entries of Kij is 1, so that K 2 = 2 and sA = 1.96 · 10 7 while for the black curve the variance of the entries of Kij is 3, so that K 2 = 4, and sA = 0.98 · 10 7 . genic environment. The autocorrelation function of this Ornstein-Uhlenbeck process is (D3) We pick the steady-state value of the initial fitness distribution to cancel the first in Eq. D3, f i (0) 2 = γ 2 /λ and obtain (conditioned on the integral of the net growth rate f + f 0 being positive so that the clone does not go extinct). Setting x = log C, we obtain a new set of Langevin equations where the birth-death noise is now treated in theÎto convention. The corresponding Fokker-Planck equation for the distribution of fitness and clone size at time t, ρ(x, f, t), verifies where s(x, f ) is the source of new clones. We solve this equation numerically using finite element methods to obtain clone size distributions for the clone-specific fitness model.

Appendix E: The Ornstein Uhlenbeck process and maximum entropy
In this Appendix we show that the maximum entropy or maximum caliber process with autocorrelation function x(t)x(t + s) = A 2 e −λ|s| corresponds to the Ornstein-Uhlenbeck process. We consider this continuous maximum entropy process as the continuous limit of a simpler maximum entropy system in discrete time. Burg's maximum entropy theorem [47] states that the maximum entropy process in discrete time that constrains X n (t) 2 = A 2 and X n (t)X n+1 (t) = A 2 e −λτ corresponds to the following Markovian dynamics: where η is Gaussian white noise. In the limit of τ → 0 we recover the constrained autocorrelation function in the vicinity of s = 0 + : x(t) 2 = A 2 , (d/ds) x(t)x(t + s) | s=0 + = −λA 2 , and Eq. E1 converges to an Ornstein-Uhlenbeck process.
Appendix F: Model solution for white-noise clone-specific fitness fluctuations In the limit of infinitely quickly fluctuating environments, γ → +∞ and λ → +∞ while keeping their ratio σ = γ/λ constant, the autocorrelation of the fitness noise approaches a Dirac delta function, and the fluctuating part of the growth rate f i (t) converges to Gaussian white noise, f i (t)f i (t ) = 2σ 2 δ(t − t ). Effectively the immune cell dynamics are now described by a one dimensional Langevin equation for the clone size where η i (t)η i (t ) = δ(t − t ) follows the Stratanovich convention and ξ i is as before. The equation for the logarithm of the clone size x = log C is We explicitly checked that the numerical solution to the clone specific fitness model in Eqs. D1 and D2 converged to the dynamics described by Eq. F1, as demonstrated in Fig. S4. We now solve this equation analytically, starting with the case of no birth-death noise: Eq. F1 simplifies to The equation for x = log C (using the Stratanovich convention) is with the corresponding Fokker Planck equation where s(x) is the source term describing the size of newly introduced clones. Assuming a constant initial clone size, s(x) = s C δ(x − x 0 ), the steady state solution is where we have defined and K is an integration constant. Imposing that ρ vanishes at infinity sets K = s C σ 2 and the final form of the steady state clone size distribution is In all simulations and solutions we find that for large clones, the model of temporally correlated fitness fluctuations behaves as the its white noise limit. This behaviour can be explained by the fact that large clones need a long time to become large. At these long timescales, the characteristic time of noise correlation is negligible and the noise may be approximated as white. For this reason, the exponent α of the power law computed assuming a white noise for the fitness fluctuations is still valid even when that noise is actually correlated in time.
Next, we re-introduce the birth-death noise and solve the general equation. The Langevin equation for x = log C, results in the Fokker-Planck equation for the distribution of clone sizes Assuming that the initial size is constant, the steady state solution is given by the solution of the inhomogeneous linear equation: The full solution is the sum ρ = ρ 0 + ρ 1 of the particular solution, and the solution ρ 1 to the homogeneous equation of solution: we set K = s for convergence and obtain the steady state clone size distribution for large x or in terms of the clone size We see that the white noise solution with birth-death noise has the same large clone power law behaviour as without birth-death noise. Fig. S5 illustrates how birth death noise in the clone-specific fitness models with time correlated noise also does not affect the power law exponent but only the cut off of the power law.

Appendix G: Data analysis
In the main text we report values of the power law exponents and power law cut off values obtained from the high throughput sequencing repertoire study of clone size distributions of zebrafish B-cell heavy chain receptors of Weinstein et al. [2]. We extracted the power law exponent and the best fit for the starting point of the power law, defined as its lower bound cutoff, from the discrete clone size distributions plotted in Fig. 1 of the main text using the methods discussed by Clauset and Newman [30]. Specifically, for each point of the cumulative clone size distribution we compute an estimate of the power law exponent with that point as cutoff (i.e the best fit of the power law including only the values of the distribution above that point) using where C min is the cut off and n is the number of points with y-axis values above C min . For each of these cut-off values we compute the Kolmogorov-Smirnov distance between the data and the estimated power law distribution: where the maximum is taken over all values above the cut off C min , F d is the cumulative distribution function (CDF) of the data and F e (C; C min ) is the CDF of the estimated power law distribution with C min as a cutoff, using Eq. G1. The the cut off is taken to be the minimum of this distance over all possible cut off values and the exponent is the exponent found for this value. The obtained power law parameters are presented in Table I. The power law exponent gives reproducible values for different individuals and agrees with values of the same exponent obtained from human data [31]. We note that the power law exponent of the cumulative distribution function is α for a power law distribution with exponent 1 + α. As discussed in detail in the main text, the reliability of the cutoff estimate C * is sensitive to experimental precision of capturing the rare clones. In the presented dataset the reads were not barcoded and the counts had to be renormalized by a known PCR amplification factor. Therefore, these normalized counts could not to used as normal counts, making the definition of a cut-off clone size problematic. To overcome this problem, we estimate the power law cut-off from the value of the cumulative distribution function at the cut-off clone size (instead of the cut-off clone size itself). That value is invariant under rescaling of absolute clone size values, unlike C * .
We notice that the steady state solution is invariant under a full rescaling of time in the equations of the dynamics. This means that the system can be described by two dimensionless parameters, α = f 0 λ 2 /γ 2 and λ 3 /γ 2 , and the introduction size C 0 . Fitting α to data and assuming value for C 0 , we can compare the value of the power law cut-off in data and in simulations to fit the remaining dimensionless parameter, λ 3 /γ 2 . Estimating f 0 based on thymic output we can predict the order of magnitude of λ and γ.

Appendix H: Cell specific simulations
In the "A model of fluctuating phenotypic fitness" Results section of the main text, we present results of Fokker-Planck simulations for the cells dynamics. Here we verify that the stochastic dynamics of cells subject to  a fluctuating cell-specific fitness are well approximated at the population level by a Fokker-Planck equation with a source term accounting for the import of new clones by comparing its numerical steady-state solution obtained by a finite elements method to explicit Gillespie simulations. We simulated the dynamics of clones using a Gillespie algorithm where cell division and death are accounted for explicitly and depend linearly on a fitness f c (t) fluctuating according to Eq. 7. The death rate is kept constant (above the average birth rate) and the fluctuations of the fitness only affect the birth rate (with the constraint that the birth rate is always positive). The agreement between the results of this detailed simulation and the Fokker-Planck solution, shown in Fig. S6, validates the linear-noise approximation for the birth-death noise as well as the averaging argument leading to Eq. 8 and 9. This allows us to rely on the Fokker-Planck solution to explore parameter space.
Appendix I: Model of cell-specific fitness fluctuations, and its limit of no heritability The cell specific fitness model described in the "A model of fluctuating phenotypic fitness" Results section of the main text arises as a description of a population where each cell experiences its own growth fluctuations but cells deriving from the same lineage remain correlated. In this Appendix we derive the equations that describe the dynamics of clones in this system.
Each cell c experiences a time-correlated multiplicative noise from environmental growth factors. For cells j in a given cell lineage (or clone) i, each individual cell's fitness follows the stochastic dynamics: . Averaging over all cells in the clone, we obtain where f i is the average fitness in clone i and where we have added a birth-death noise term (µ + ν)C i ξ i . We use theÎto convention for the birthdeath noise, ξ i (t)ξ i (t ) = δ(t − t ) and the Stratanovich one for the environmental noise η i (t)η i (t ) = δ(t − t ). The equivalent equations for x = log C are and the Fokker-Planck equation is where s(x, f ) is the joint distribution of size and fitness or newly arriving clones (from thymic or bone marrow output). This is the full Fokker-Planck equation that is solved numerically in the main text using the finite elements method.
Because of the 1/ √ C i prefactor in front of the noise term, we could expect fitness fluctuations to behave like a birth-death noise in the limit of low heritability (λ → ∞). In the remainder of this Appendix we show that this is not the case, and we show how to take the limit of no heritability properly.
Consider the limit of λ → ∞ and γ c → ∞, keeping the ratio γ c /λ constant, so that f does not become infinitesimally small. The equation for the environmental stimulation f in x = log C space is given by (in Stratanovich convention) With infinite precision, for any value of t, we set the integral of η to be bounded and obtain the first integral is with probability 1 − smaller in norm than √ 2γ c √ tK( )e −k , where K(t) is a constant to control the variations of the integral of ξ with probability (where K(t, ) = Φ −1 (1 − In Appendix I we have shown that in the limit of λ → ∞ and γ c → ∞, the system reduces to the one dimensional equation with theÎto rule for the white noise η i . The corresponding Fokker-Planck equation is Assuming a deterministic introduction size s(x) = s C δ(x − x 0 ), at steady-sate we get which for x > x 0 is solved by where K is an integration constant, Ei is the exponential integral function and The divergence of Ei at infinity sets K = s C λ 2 /(γ 2 c ) and the clone size distribution is ρ(x) = Ei(e x /C m ) − Ei(C −1 m ) e −e x /Cm+x for x < x 0 Ei(e x0 /C m ) − Ei(C −1 m ) e −e x Cm+x for x > x 0 (J10) or in terms of x = log C ρ(C) = e −C/Cm Ei(C/C m ) − Ei(C −1 m ) for C < C 0 e −C/Cm Ei(e x0 /C m ) − Ei(C −1 m ) for C > C 0 (J11) The validity of this solution is checked in Fig. S9 and the convergence of the full solution of Eq. I6 (with no birthdeath noise) to the analytical solution in the limit of no heritability (λ → ∞) is show in Fig. S10.
For comparison, in a pure birth-death process (no fitness fluctuations) the clone-size distribution is, for C large enough, ρ(C) ∼ e −C/Cm /C where C m = (µ + ν)/(2(µ − ν)), as shown in Appendix A. These two solutions both have an exponential cutoff, but have very different power-law exponents, corresponding to α = 0 and α = −1, respectively.
We now add the birth-death noise, i.e. consider both types of noise, still in the limit of no heritability. The corresponding Langevin equation reads: and gives the solution: with which is a power-law with an exponent 0 ≤ 1 + α ≤ 1 and an exponential cutoff The convergence of the solution of the full system, Eq. I6, to this solution is checked in Fig. S11.