A branching stochastic evolutionary model of the B-cell repertoire

We propose a stochastic framework to describe the evolution of the B-cell repertoire during germinal center (GC) reactions. Our model is formulated as a multitype age-dependent branching process with time-varying immigration. The immigration process captures the mechanism by which founder B cells initiate clones by gradually seeding GC over time, while the branching process describes the temporal evolution of the composition of these clones. The model assigns a type to each cell to represent attributes of interest. Examples of attributes include the binding affinity class of the B cells, their clonal family, or the nucleotide sequence of the heavy and light chains of their receptors. The process is generally non-Markovian. We present its properties, including as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t\rightarrow \infty $$\end{document}t→∞ when the process is supercritical, the most relevant case to study expansion of GC B cells. We introduce temporal alpha and beta diversity indices for multitype branching processes. We focus on the dynamics of clonal dominance, highlighting its non-stationarity, and the accumulation of somatic hypermutations in the context of sequential immunization. We evaluate the impact of the ongoing seeding of GC by founder B cells on the dynamics of the B-cell repertoire, and quantify the effect of precursor frequency and antigen availability on the timing of GC entry. An application of the model illustrates how it may help with interpretation of BCR sequencing data.


Introduction
During an infection, the immune system initiates a response to protect the body from the invading pathogen and establish immunity against future reinfections.This response often relies on B cells, a class of lymphocytes specialized in the production of antibodies that neutralize invaders by binding to foreign target molecules on the pathogen called antigens.Since each antibody binds to specific molecules, a highly diverse B-cell repertoire is essential for achieving robust immune protection.While naive B cells express antibodies with considerable molecular diversity due to somatic recombination (Dreyer and Bennett 1965;Tonegawa 1983), the B-cell repertoire must undergo additional diversification to effectively respond to a vast array of evolving pathogens (Janeway et al. 2005).
The B-cell repertoire further diversifies during affinity maturation, the Darwinian evolutionary process that occurs in germinal centers (GC) (Janeway et al. 2005;MacLennan 1994;De Silva and Klein 2015).GC are temporary structures that develop in secondary lymphoid tissues such as lymph nodes, spleen, tonsils, and Peyer's patches in the gut.They provide the microenvironment that supports and regulates adaptation of the B-cell repertoire to the invading pathogen.They are continually seeded by B cells selected for their ability to bind to the antigen (Schwickert et al. 2007).They include two distinct anatomical compartments, the dark and light zones, between which B cells traffic back and forth to undergo successive rounds of proliferation and somatic hypermutation followed by antigen-mediated selection (Eisen and Siskind 1964;Weigert et al. 1970;Jacob et al. 1991;Muramatsu et al. 2000).Somatic hypermutation occurs at an extraordinarily high rate estimated at ∼ 2 × 10 −4 − 10 −3 per base pair per generation which is about a million times greater than the mutation rate observed in other parts of the genome (Berek and Milstein 1987;McKean et al. 1984).
Understanding the rules of affinity maturation has major clinical applications.For example, several clinical trials currently in progress are testing novel HIV vaccines that seek to harness this process to elicit B cells able to secrete broadly neutralizing antibodies (bnAbs) similar to those isolated in people living with HIV (Leggat et al. 2022).By deciphering the rules governing antibody maturation, immunologists could manipulate GC to devise effective immunogens (i.e., molecules capable of inducing an immune response) and vaccination strategies aimed at preventing HIV acquisition and infection against other pathogens.
Mathematical models have played a crucial role in advancing our understanding of the dynamics of GC B cells, shedding light on this complex evolutionary system.Notably, pioneering work in the field can be found in (Agur et al. 1991;Kepler andPerelson 1993a, b, 1995;Oprea and Perelson 1997;Kleinstein and Singh 2001;Meyer-Hermann et al. 2001;Iber and Maini 2002) while a recent review is provided in Buchauser and Wadermann (2019).Many of these models are formulated as agentbased models that evaluate properties of GC via simulations.While this approach considers the stochastic nature of GC reactions, it may be computer-intensive and only yield conclusions for chosen parameter values.In this paper, we propose a comprehensive stochastic framework that allows for a more systematic exploration of properties across a range of plausible parameter values.We demonstrate the utility of this framework by showing that it predicts established properties of GC, while offering explanations for these properties and enabling the identification of new ones.
Somatic hypermutation and the competition between GC B cells are two central factors that structure evolution and adaptation of the B-cell repertoire.We propose a model that captures the most salient features of the process by assuming that the dynamics of GC B cells is controlled by three mechanisms: (1) the continual recruitment of founder B cells; (2) clonal expansion; and (3) somatic hypermutation of the BCR-encoding immunoglobulin (Ig) gene loci.
The continual recruitment of founder B cells induces an influx of B cells into GC and initiates the formation of clones (Schwickert et al. 2007(Schwickert et al. , 2009;;Turner et al. 2017).We model this mechanism using a counting process, potentially time-inhomogeneous.Evolution of B-cell clones is driven by cellular division, death, and differentiation, which we describe using an age-dependent branching process.Finally, accumulation of somatic hypermutations, the mechanism which randomly alters the binding properties of BCR, is captured by assigning a type to each cell and defining a network specifying the set of admissible transitions between types as cells divide.It is important to note that types need not represent nucleotide or amino acid sequences of BCR, but could instead reflect associated properties such as binding affinity classes or clonal families.Differences in the probability of death between types capture differential fitness (competition) between different types of B cells.The resulting model is a multitype age-dependent branching process with immigration, potentially time-varying, extending previously considered models (e.g., Sevastyanov 1957;Jagers 1968;Hyrien et al. 2005Hyrien et al. , 2017Hyrien et al. , 2015;;Pakes and Kaplan 1974;Mitov et al. 2018;Slavtchova-Bojkova et al. 2023a, b;Yakovlev and Yanev 2006).The model allows tailored network architectures to study evolution of a variety of features of the BCR repertoire.Numerical methods have been proposed to implement these models (e.g., saddlepoint approximations to moments Hyrien et al. 2010).
We demonstrate the flexibility of the framework by studying three model structures.The first one classifies B cells according to their binding affinity for a given antigen (see Sect. 3.5.1).We show that, as a class, lower binding affinity B cells have a competitive advantage over their higher-affinity counterparts because of their propensity to engage in GC reactions earlier, attributed to their higher abundance.We also find that clonal families should contain B cells with diverse binding affinity levels that grow in size at the same rate within their respective clones, a finding aligning with previous observations (Kuraoka et al. 2016;Bannard and Cyster 2017).No definitive mechanisms have been put forth to explain this phenomenon; our model shows that it could be, in part, attributed to the combined effects of clonal expansion and somatic hypermutation (Sect.3.7).Mathematically, this property emerges from the communication between cell types (or, equivalently, the irreducibility of the mean offspring matrix).From a biological standpoint, it results from the fact that somatic hypermutation may alter affinity for the antigen, and that these alterations can be reversed.
Several studies have shown that early GC harbor B cells from multiple clonal families, potentially hundreds, and gradually loose this clonal diversity to become oligoclonal (i.e., dominated by a few clones) (Kroese et al. 1987;Janeway et al. 2005;Faro and Or-Guil 2013;Tas et al. 2016).To study evolution of clonal dominance, we consider a second model in which each type identifies a clonal family.Analysis of the Fig. 1 Principle of a successful sequential immunization strategy.The vaccine regimen begins with a priming immunogen seeking to induce B cells with genetic and molecular signatures expected to allow HIV bnAb production.B cells induced by the prime consist of both desired and non-desired (on-and off-target) B cells that express BCR with limited somatic hypermutations and lack the ability to neutralize the virus.To induce B cells with neutralizing activities, the immune system is boosted by a series of immunogens designed to support further maturation of B cells from clonal families induced by previous vaccine doses and select somatic hypermutations relevant for broad neutralization against HIV.After each vaccination, on-target B cells compete with off-target B cells for antigen-mediated selection model supports that oligoclonality arises from differential binding affinity between clones.It also suggests that clonal dominance may be a temporary feature repeatedly lost to newer clones until the GC reaction stops (see Sects. 3.5.2,5).
We introduce a third model designed to examine the divergence of BCR sequences from their unmutated (germline) ancestors, driven by the accumulation of somatic hypermutations.This model is relevant for characterizing the evolution of the B-cell repertoire in the context of immunization.Most vaccines achieve protection through induction of antibodies (Plotkin 2008).Protective vaccines against the human immunodeficiency virus (HIV) have remained elusive to date.However, about 20-30% of people with HIV infection develop bnAbs which neutralize many of the globally circulating HIV strains through binding conserved regions of the virus (Doria-Rose et al. 2009).VRC01 is one bnAb that specifically targets the CD4 binding site on the HIV virus.The efficacy of passive immunization with VRC01 against HIV acquisition was recently tested in two clinical trials which showed that protection against HIV infection by VRC01 was confined to viruses highly susceptible to neutralization (Corey et al. 2021).Novel vaccines that seek to induce HIV bnAbs via sequential immunization are currently being evaluated (Haynes and Mascola 2017; Leggat et al. 2022) (see Fig. 1 for a description of their principle).One challenging goal of this strategy is selection of somatic hypermutations similar to those of HIV bnAbs, including rare substitutions, capable of neutralizing the virus.One metric for assessing the vaccine-induced immune response is the distance between the sequences of vaccine-elicited BCR to Each cell is assigned a type at birth; when it divides it may produce offspring of different types (e.g., reflecting the impact of somatic hypermutations).In the example depicted here, the population includes two types of B cells, but more complex models with arbitrary number of types may be constructed that of a designated target bnAb.The third model describes the progression of these distances over time (Sect. 6).
The remainder of this paper is organized as follows.Section 2 provides a brief review of how antibodies are encoded and their structure.The general model is defined in Sect.3; see Fig. 2 for a graphical overview.Section 3.2 discusses the model that describes the seeding of GC.Section 3.3 focuses on the particular case where this model is a time-inhomogeneous Poisson process.Section 3.4 presents the multitype branching process that describes clonal expansion and the impact of somatic hypermutation.The specific models introduced earlier are defined in Sect.3.5.Model properties are presented in Sects.3.6, 5 and 6 to explore patterns of the dynamics of GC reactions.Our study of clonal dominance relies on measuring the alpha and beta diversity of B-cell repertoires; these indices are presented in Sect. 4 in the general context of multitype branching processes.Throughout, results are discussed and illustrated in the context of sequential immunization, introduced in Sect.6.1.Technical details (proofs) are provided in the Appendix.

Encoding and structure of antibodies
Antibodies consist of four polypeptide chains arranged in two identical pairs, with each pair consisting of a heavy and a light chain connected by disulfide bonds (Fig. 3).The stem of the (Y-shaped) antibody molecule, called the constant region, determines the class of the antibody (IgG, IgA, IgM, IgD, and IgE).The tips of the antibody molecule are known as the antigen-binding regions or variable regions.The variable regions of the antibody molecule contains regions called complementarity-determining regions (CDR).There are three CDRs in both the variable heavy and variable light chains, denoted as CDR1, CDR2, and CDR3.These loops are the most variable parts of the antibody, and their unique sequences are responsible for antigen recognition.The unique sequence of amino acids in the CDRs allows each antibody to bind to a specific target antigen or a closely related group of antigens.Between the CDRs, there are four regions known as framework regions (FR1-FR4) which provide the structural scaffold for the variable region.While they are less variable than the CDRs, variations in the nucleotide sequence of the FRs can still affect the antigen-binding properties of the antibodies.The combination of the CDRs and the adjacent FRs forms the antigenbinding site of the antibody.
The genetic information encoding the variable regions is divided into multiple variable (V), diversity (D, only for heavy chains), and joining (J) gene segments.During B cell development, the DNA of the B cell randomly rearranges one gene segment from each V(D)J category per chain.The selected gene segments are joined together, and any excess DNA between them is removed.This process which takes place in the bone marrow is called somatic recombination.It results in a functional gene encoding the variable region of the antibody.Variable regions are approximately 110-130 amino acid long in heavy chains, and slightly shorter in light chains.Most positions are encoded by the V gene.
After the B cell encounters an antigen, it may undergo somatic hypermutation in the variable regions of its Ig genes.This process introduces random mutations into these genes, primarily in the CDRs.B cells with mutated antibodies that bind more effectively to the antigen receive a survival advantage and are selected for further proliferation, eventually leading to the production of high-affinity antibodies.Somatic hypermutation further improves these antibodies for optimal antigen recognition.This process produces more effective antibodies over time during an immune response.See Janeway et al. (2005) for further details.

Classification of B cells into types
All BCR expressed on the surface of a B cell (around 120,000; Alt et al. 2015, p. 154) share the same nucleotide sequence that determines their antigen specificity (Janeway et al. 2005;Alt et al. 2015).During a GC reaction, proliferating B cells accumulate somatic hypermutations in their Ig gene loci, causing BCR sequences to drift away from the original (naive) sequence assembled during somatic recombination.This evolutionary process, which enables adaptation of the B-cell repertoire, may be studied from many angles.Aside from the sequence of the B-cell receptors, relevant properties (attributes) include their binding affinity class and clonal family.We propose a unified modeling framework in which a generic type is assigned to each B cell.These types are indexed by a set K, and we let K = |K| denote the total number of types, possibly infinite, in the model.Types represent arbitrary attributes, but it may be convenient  Schwickert et al. (2007) experimentally showed that GC are open structures that are continually visited and seeded by naive follicular B cells (Schwickert et al. 2007).We describe the ongoing seeding of GC using K type-specific immigration processes, one per cell type.These processes are formulated by defining K sequences of positive and increasing random variables (r.v.) {T k , = 1, 2 . ..}where 0 ≤ T k1 ≤ T k2 • • • a.s., k ∈ K.Each of these sequences represents the successive time points at which type-k founder B cells join the GC reaction.The first type-k founder B cell enters the GC at time T k1 , the second one at time T k2 , and so on.The sequence {T k , = 1, 2 . ..} generates a counting process k (t) = ∞ =1 1 {T k ≤t} , t ≥ 0, which represents the number of type-k founder B cells that have entered the GC by time t.Let {T , = 1, 2 . ..} = k∈K {T k , = 1, 2 . ..} denote the collection of all time points at which founder B cells enter the GC, and define the overall immigration process

Seeding of germinal centers as a counting process
which represents the total number of founder B cells that have joined the GC by time t, regardless of type.The specification of the counting processes { k (•)} k∈K is contextdependent.Given its importance in applications, the Poisson process is treated in detail in the next section.An example where they are not defined as Poisson processes is given in Sects.3.5.2and 5.

A time-inhomogeneous Poisson race
Throughout this section, we assume that: The rate at which founder B cells enter GC has not been fully elucidated, and each k (•) could be formulated as a time-homogeneous Poisson process for simplicity.The ability of naive antigen-specific B cells to enter pre-existing GC is impacted by multiple factors, including antigen availability.Turner et al. (2017) showed that B cells that acquire a small amount of antigen can enter GC at any stage of the reaction; however, naive antigen-specific B cells may be recruited during a more restricted time window (within 6-10 days) (Turner et al. 2017).Taken together, these findings support specifying k (•) as a time-inhomogeneous Poisson process with rate r (•) functionally related to antigen availability.
We refer to Resnick (1992) and Durrett (2004) for discussions on Poisson processes, and recall the following useful five facts about them: To identify the type of the -th founder, define a r.v.I = (I k , k ∈ K) where I k denotes the number of type-k founder B cells that entered the GC at time T , = 1, 2 . ... We assume that a single B cell immigrates at each time T ; hence, I k ∈ {0, 1} and k∈K I k = 1 with probability one.We also assume that {I } ∞ =1 are independent and identically distributed (i.i.d.) r.v. with p.g.f.g(s) = E s I = k∈K g k s k where s α = k∈K s α k k .On occasion, we will make the following assumption: (A2) There exists positive constants {g k } k∈K such that k∈K g k = 1 and r k (t) = g k r (t) for every t ≥ 0 and k ∈ K.
Under Assumption (A2), temporal fluctuations in the influx of founder B cells remain identical across types, up to the multiplicative constants g k .Together with Fact 4, this implies that g k = r k (t)/ i∈K r i (t) is the probability that a founder B cell that joins the GC at time t is of type k.Assumption (A2) implies that these probabilities remain unchanged throughout time, even when r (•) is time-dependent.When each process k (•), k ∈ K, is a time-homogeneous Poisson process with rate g k r 0 for some constant r 0 ∈ (0, ∞), then (•) is a time-homogeneous Poisson process with rate r (•) ≡ r 0 .
For every k ∈ K, the r.v.T k1 denotes the waiting time until the first type-k founder enters the GC.According to Fact 5, T k1 has a distribution with hazard rate r k (•).If immigration is time-homogeneous, with r k (t) = g k r 0 , t ≥ 0, then T k1 ∼ E x p(g k r 0 ).If r k (t) = g k r 0 t γ for some constants r 0 > 0 and γ > −1, then T k1 has a Weibull distribution with shape and scale parameters 1 + γ and 1+γ g k r 0 1 1+γ , heavy tailed when the immigration rate decreases over time at a polynomial rate γ ∈ (−1, 0).

Its median and mean are 1+γ
, and the ratio of the medians or means of the waiting time until first type-k and first type-k founders join the GC are both given by (g k /g k ) 1 1+γ .We note that the assumption of a decreasing immigration rate may partly reflect current vaccination strategies where doses are typically administered as a bolus.
Figure 4 shows the median, mean, and ratio of medians or means of T k1 as a function of γ , for various values of the precursor frequency g k and of g k /g k .As g k decreases, both the median and mean increase.The median is neither a decreasing nor monotone function of γ ∈ (0, ∞): it increases for γ ∈ (−1, (g k r 0 e 1 − ln 2)/ln 2) and decreases for γ ∈ (g k r 0 e 1 − ln 2)/ln 2, ∞).The ratio comparing the medians for two cell types (less versus more abundant) decreases for γ ∈ (−1, ∞).When the influx of B cells recruited by the GC slows down (γ ∈ (−1, 0)), the median waiting time is disproportionately longer for rarer cell types, and shorter when the influx accelerates.For instance, with an immigration rate decreasing at a rate of γ = −1/2 (i.e., r k (t) = g k r 0 t −1/2 ), every halving of the frequency of the founder B cells (i.e., g k /g k = 2) causes the median or mean waiting time to increase by a factor of 4 Remark 1 (Biological relevance) During an immune response, the recruitment of new B cells by GC may slow down when the amount of antigen that fuels the reaction decreases (Turner et al. 2017).The time until rare B cells (which tend to join the reaction later than more common ones, as discussed below and in Remark 2) join the GC may be longer than the waiting time for more common B cells.The fueling of GC by a continual supply of antigen, perhaps characteristic of chronic infections or GC in the gut, could increase the likelihood of rare B cells entering GC.This observation has clinical applications, particularly for the development of vaccines that seek to elicit responses from B cells with rare genetic and molecular signatures which could benefit from sustained immunogen delivery.A related phenomenon has been recently reported (Lee et al. 2022).
The recruitment of founder B cells into GC can be cast as a competition between classes of B cells that race against each other to place some of their members into the GC reaction.Let denote the number of type-k 1 founder B cells, k 1 ∈ K 1 , that have joined a given GC by the time the first type-k 2 founder B cell, k 2 ∈ K 2 , joins the same GC; this time is Suppose that K 1 and K 2 are two non-overlapping subsets of K, that Assumption (A2) holds, and Moreover, . Thus, when Assumption (A2) holds, the mean number of type-k 1 founder B cells, Remark 2 (Biological relevance) The proportions {g k , k ∈ K} represent the frequencies of antigen-specific B cells among each type of B-cell precursors recruited by GC.
It has been estimated that, in humans, the frequency of VRC01-class CD4 binding site antibody precursors is about 1 in 2.4 million B cells.Precursors for other HIV-like bnAbs (e.g., PGT121) may be even more rare (Jardine et al. 2016;Steichen et al. 2016).The practical implication of the above results is that immunogens used by germline-targeting vaccines must substantially enrich the pool of on-target B-cells so they will not be outcompeted by off-target B cells selected for entry in GC induced by boost vaccination.
The continual seeding of GC has also implications on the diversification of the B-cell repertoire.Our results suggest that GC tend to be initially colonized by commonly found precursors (e.g., B cells founding clones with 'average' binding affinity potential, sufficient to enter GC) before they recruit the rarest ones (e.g., B cells able to produce high binding affinity B cells).By the time the first high affinity B cell joins the reaction, the GC will have been visited by many B cells with lower binding affinity potential.Each of these early founder B cells may initiate a clone at a time where competition for the antigen may not yet be at its peak.This mechanism may therefore support diversification of GC B cells during the earliest stage of the reaction.The waiting time until the highest binding affinity B cells join the reaction allows B cells with lower binding affinity to undergo somatic hypermutations, further diversifying the BCR repertoire.Thus, as a class, lower binding affinity B cells have a temporary competitive advantage over high binding affinity B cells due to the timing of their arrival which allows them to clonally expand.Although (rare) high binding affinity founder B cells tend to join GC later, their clones may still outnumber lower affinity clones if they have higher fitness.Another consequence of the continual seeding is clonal instability in the GC, a phenomenon studied in Sect. 5 which also diversifies the B-cell repertoire through the turnover of dominant clones.

Branching mechanism and intra-clonal evolution
For every k ∈ K, the lifespan of any type-k cell is described by a non-lattice r.v.
At the end of its lifespan, every type-k cell produces a random number of offspring described by a r.v.ξ k = (ξ k j , j ∈ K) where ξ k j denotes the number of typej daughter cells.We assume that cells either die or divide; hence, j∈K ξ k j ∈ {0, 2}.Let , s = (s 1 , . . ., s K ), denote the probability generating function (p.g.f.) of ξ k .Lastly, within a clone, every cell evolves independently of every other cell.
Let C(k) denote any type-k cell.Define p k 0 = P{C(k) → ∅} and p k i j = P{C(k) → (C(i), C( j))} where p k 0 is the probability that any type-k cell completes its lifespan without producing any offspring either because it dies, or differentiates into a plasma or memory B cell, or exits the GC; and p k i j denotes the probability that any typek cell divides into two cells of types i and j, respectively, with i, j, k denote the conditional probability that any type-k cell divides into one type-i and one typej cell.The conditional probabilities {q k i j } i, j∈K capture the impact of somatic hypermutation on type-k cells that manifests at division.Notice that jn denote the expected number of typej cells produced by any type-i cell at the end of its lifespan.Define Fig. 5 Four examples of network architectures describing the interconnection between cell types.In each example, the nodes represent the types and the edges identify the set of admissible transitions between types.For simplicity, self-renewing division (without change in type) and cell death are not shown in the graphs.Architecture 1: The process includes five types, with each type representing a binding affinity class (from lowest (1) to highest ( 5)) positioned along two mutational pathways; type-1 cells may produce type-2 and type-3 cells, and vice versa, allowing offspring to reversibly enter either one of the pathways: type-2 cells are less likely than type-3 cells to produce offspring of highest affinity (type-5); this first example is introduced in Sect.3.5.1.Architecture 2: The model includes countably many types, each type representing one of the clonal families seeding a GC, ranked in the order in which they join the GC; in this second example, introduced in Sect.3.5.2, the types do not connect.Architecture 3: Each type in the model represents the nucleotide sequence of the immunoglobulin gene loci that encode the heavy and/or light chains of the BCR; in this third example, briefly introduced in Sect.3.5.3,transitions between types arise through somatic hypermutation.Architecture 4: The model offers a simplified description of the accumulation of somatic mutations in the Ig gene loci of B cells; the state space is a 5-dimensional vector: its first two entries represent the number of mutated and unmutated positions among those at which the germline and a target bnAb match (ν 11 ∪ ν 12 ); its last three entries represent the number of mutated and unmutated positions among those at which the germline and bnAb do not match (ν 21 ∪ ν 22 ∪ ν 23 ); specifics about the model, in particular why three entries are needed in the latter case, are explained in Sect.6 the associated matrix of mean offspring and write m i = j∈K m i j for the total mean number of offspring produced by any type-i cell.Then, is the number of type-k cells at time t, k ∈ K, generated by the above-defined multitype branching process.Put for the size of the clonal family at time t.

Examples of model structures
We consider a few illustrative examples of model structures, shown in Fig. 5.

Example 1: modeling evolution of binding affinity classes
The role of GC is to generate high binding affinity B cells through somatic hypermutation, antigen-mediated selection, and clonal expansion.Our first model classifies the B cells that populate a GC into types that each represents a binding affinity class for one of the antigens fueling the reaction.The affinity level of an antibody refers to its strength in binding to a specific antigen and the stability of the resulting complex.A higher affinity indicates a stronger interaction between the antibody and antigen.Many clonal families originate from a founding B cell with a relatively low affinity, but nevertheless sufficient to secure entry into the GC.Their descendants are subjected to somatic hypermutation which may either leave binding affinity virtually unchanged or altered.When their affinity improves, B cells are more likely to persist in the GC, captured by a decrease in the probability of death ( p k 0 ).The branching mechanism and topology of the transition network are formulated to describe propagation of the members of a clonal family over the various affinity classes under a given scenario.The example shown in Fig. 5 allows two mutational pathways.Both of them enable enhancement of affinity, but one is a mutational dead-end, prematurely restricting improvement of affinity.It extends the catenary model studied in Kleinstein and Singh (2001).

Example 2: an infinite-type model of clonal dominance
In the second example, GC B cells are classified according to the founder B cell from which they descend: descendants of the first founder are all of type 1; descendants of the second founder are all of type 2; and so on.Thus, here, each type identifies a particular clonal family, types do not communicate with each other, K = {1, 2 . ..}, and K = ∞.The offspring mean matrix is and the lifespan of any type-k cell follows a distribution with c.d.f.G k (•).We use this model to study the dynamics of clonal dominance within GC.The model predicts that GC are initially clonally highly diverse, before gradually loosing their clonal diversity and converging toward dominance by a few clonal families.These predictions corroborate previous experimental observations.Additionally, clones that dominate a GC reaction are eventually replaced by new ones, until the reaction stops.

Example 3: mutational gains and losses relative to a germline gene and a target antibody
Evolution of the BCR repertoire may be described at the sequence level by building a model in which each type represents the nucleotide sequence of a particular BCR and postulating how somatic hypermutations induce transitions between types at division (Fig. 5, architecture 3).However, the complexity of the rules of somatic hypermutation and affinity-based selection makes the formulation of such a model difficult (requiring specification of ≥ 4 L (4 L − 1) transition probabilities for nucleotide sequences of length L).In this paper, we propose a simpler evolutionary model that focuses on gains and losses of mutations relative to a germline and a target antibody sequence.
The model partitions positions of the heavy and/or light chain sequence into 5 distinct subsets to capture mutations at positions at which the bnAb and germline match or differ (Fig. 5, architecture 4; see also Sect.6).

Properties of the general model
Let X , = 1, 2 . .., be i.i.d.r.v.representing the type of the -th founder B cell.When assumption (A1) holds, it follows from Fact 4 that denotes the number of type-k cells in the -th clone t units of time after its initiation, and with initial condition Z ( ) (0) = e X where e k is a K -dimensional vector with all entries equal to 0 except for a one in the k-th place.We assume that: (t) |X 1 = k denotes the conditional p.g.f. of Z (1) (t), given the process begins with one type-k cell.These p.g.f. are the unique solutions of the nonlinear integral equations: Mode (1971), Athreya and Ney (1972)].
For every t ≥ 0, define A k j (t) = E(Z (1) where, for every i, j, k ∈ K, A k j (t) denotes the average number of typej cells at time t in a clone started from a single type-k founder B cell of age 0 at time 0, B k i j (t) is the second-order factorial moment of the number of type-i and typej cells at time t in a clone started from a single type-k cell at time 0, and and If M * (s) is irreducible, then its Perron-Frobenius root ρ * (s) always exists, and we can define [see Mode (1971)]: Definition 1 (Malthusian parameter) The Malthusian parameter α is the solution, assuming it exists, of the equation ρ * (α) = 1 where ρ * (α) denotes the Perron-Frobenius root of M * (α) (i.e., ρ * (α) is the largest eigenvalue of M * (α)).The Malthusian parameter α always exists in the super-and critical cases; its existence must be verified in the subcritical case.When α exists, we have α < 0, = 0, or > 0 depending on whether the process is subcritical, critical, or supercritical, respectively.
Suppose that M is positive regular (i.e., ∃n ∈ {1, 2 . ..} such that the entries of M n are all strictly positive), and the p.d.f.
Remark 4 Intra-clonal dynamics is affected by whether the clonal family ultimately survives.Assuming (A4) holds, and letting q = P{W = 0} denote the probability of extinction, we have conditional on non-extinction that Otherwise, the process converges to 0 as t → ∞ (extinction).
To describe the dynamics of GC B cells, define the vector denotes the number of type-k B cells at time t in the GC.We assume that the GC does not contain any B cell at time t = 0, and set Y(0) = 0 where 0 denotes the K -dimensional null vector.Define Y (t) := K k=1 Y k (t).For every t ≥ 0, Y(t) can be expressed as where k (u) k∈K represents the composition of the -th clone u units of time after it started, setting Z ( ) (u) = 0, u < 0. We refer to Y(•) as a K -type age-dependent branching process with non-homogeneous immigration.
Let (t; s) = E{s Y(t) |Y(0) = 0} denote its p.g.f., denote the expectation and covariance of Y(t).Adapting a result from Mitov et al. (2018), we deduce when Assumptions (A1,A2) hold that with the boundary condition (0; s) = 1.Moreover, for every j, k ∈ K, where A j (x) = i∈K g i A i j (x) and C jk (x) = i∈K g i B i jk (x).

Classes of B cells that communicate do not vanish
Under Assumption (A4), the proportion of type-k cells in a clonal family not going extinct tends a.s. to v k / K k =1 v k as t → ∞.Since v has strictly positive entries, all types are predicted to grow at the same exponential rate α, regardless of the probability at which they die or self-renew.This property arises from the irreducibility of the mean offspring matrix M (i.e., the fact that all types communicate).It is illustrated in Fig. 6 using simulations based on a two-type process.From a biological standpoint, communication between types is driven by somatic hypermutation.Its biological interpretation in the context of Model 1 (Sect.3.5.1) is that clonal families in which high binding affinity B cells can mutate into lower binding affinity B cells (and vice versa) may include non-negligible numbers of B cells of varying binding affinity levels, even those less likely to survive.Kuraoka et al. (2016) made a similar observation in mouse experiments, reporting clonally related B cells competing for a same epitope but with affinity orders of magnitude different (Kuraoka et al. 2016).They observed cases where low-affinity B cells accounted for 22-27% of GC B cells.This observation, which has been interpreted as an indication of the permissiveness of GC to retain low or moderate affinity B cells (Kuraoka et al. 2016;Bannard and Cyster 2017), is not fully understood.While several factors could explain it, including the possibility of an affinity threshold (Kuraoka et al. 2016), the asymptotic behavior of Z(t) indicates that clonal expansion and somatic hypermutation could also jointly contribute.

On selection and competition between B cells
Several processes contribute to the regulation of B cells in GC.Homeostasis maintains a balance between B-cell amplification and excessive immune activation.Additionally, antigen-mediated selection favors B cells that have undergone mutations that result in higher affinity for the antigen.These selected B cells receive stronger survival signals compared to those with lower affinity, which are more prone to undergoing apoptosis.
The model accounts for some, but not all, of these processes.For instance, it captures, at least partially, variations in antigen-mediated selection among B cells by introducing type-specific probabilities of cell death ( p 0 k ).These probabilities can differ, resulting in some cell types in the model having lower fitness and being more susceptible to cell death than others.
However, the model does not consider the homeostatic mechanisms that regulate the growth of GC B-cell populations.Specifically, in cases where a clonal family is driven by a supercritical process that does not go extinct, its size will grow exponentially.Therefore, when describing the composition of a population based on cell counts, the model is best suited for depicting the organization of GC or individual clonal families during their early phase.Over extended time periods, the model can be used to track the temporal evolution of the relative frequencies of B cell classes using Y(t)/Y (t) or Z(t)/Z (t).
Another critical aspect of GC dynamics involves the competition among B cells, both between and within clonal families and types.When a new B cell with higher affinity than those already present in the GC joins the population, competition could increase the likelihood of death of other, less competitive cell types, resulting in a reduction of their Malthusian parameter.While the model does not explicitly detail this phenomenon, it partially captures it when comparing the growth of two distinct cell types.To illustrate, consider two different, non-communicating cell types ing the odds of selecting a cell of type k versus type k , provides insights into how the frequencies of cell types k and k evolve over time.We use it in Sect. 5 to compare the size of clonal families.For every t ≥ max(T k , T k ), conditional on non-extinction for both clonal families (W k W k > 0), this ratio satisfies: Thus, the odds would remain similar to those in a model where the Malthusian parameter in the k-th clone decreases upon the initiation of the more competitive k -th clonal family, provided that the difference between the Malthusian parameters is α k − α k .Overall, despite some limitations, we have observed that the model can effectively replicate certain experimental observations.

Temporal diversity of the B-cell repertoire
A diverse BCR repertoire is essential for effective immune protection against infection and a potential benefit for responses to antigens structurally related to those from past infections or vaccinations (Kuraoka et al. 2016).However, it can also give rise to significant amounts of B cells expressing binding antibodies unlikely to progress toward pathogen neutralization; instead, these B cells may compete with and restrict those capable of generating neutralizing antibodies.We propose to measure evolution of diversity over time in terms of the K types of the model to gain further insights into the organization of the B-cell repertoire.In this regard, we introduce indices of temporal α and β diversities for multitype branching processes.We applied them in Sect. 5 to study the turnover of clonal families within GC.

Temporal alpha diversity
Alpha diversity is a fundamental concept in ecology to quantify species diversity within a given habitat.Popular alpha diversity indices include the Hill numbers.When applied to the proposed multitype branching process with immigration to quantify the diversity of types in a GC at time t, the Hill number of order q, defined on the event where q ∈ (0, 1) ∪ (1, ∞) determines the sensitivity of the index to the relative frequency of each type, with the influence of the most abundant types increasing with q (Hill 1973).The higher D q Y (t), the higher diversity.Letting q → 0, (species richness), representing the total number of types in the GC at time t.Setting q = 2 yields the inverse Simpson diversity index (Simpson 1949) used in Sect. 5 to assess evolution of clonal diversity, whereas letting q → 1 gives the exponential of the Shannon index: Clone-specific Hill numbers are defined on . Since the summands in D q Y (t) and D q Z (t) are positive, both D q Y (t) and D q Z (t) are well defined when K = ∞.
= K as q → 0, and The proof follows easily from Assumption (A4) and Proposition 1 in Hill (1973).

Temporal beta diversity
The variation in species composition between two habitats or time points may be measured using beta diversity indices (Whittaker 1960(Whittaker , 1972)).We consider these indices to quantify the degree of temporal dissimilarity (turnover) of the K types in a clonal family and a GC.These indices are defined using either presence-absence, or abundances (counts), or relative abundances (proportions) of types.The Jaccard distance (Jaccard 1912) between the collection of types present in a clonal family at times t 1 and t 2 is defined using absence-presence of types as where S(t) = {k ∈ K : Z k (t) > 0} is the set of types in the clone at time t.A value of 0 indicates identical type representation at times t 1 and t 2 , and 1 no overlap in types.This index does not take frequency of types into account in measuring diversity.The Bray-Curtis index, which quantifies the degree of compositional dissimilarity within a clonal family at times t 1 and t 2 using is defined based on the vector of counts Z(t) and normalized by Z (t 1 ) + Z (t 2 ) to range between 0 (no compositional variation) and 1 (full turnover of types).A popular beta diversity index based on relative abundances to measure the overlap between the composition of a clonal family at times t 1 and t 2 is the Hellinger distance is the Bhattacharyya coefficient, representing the cosine of the angle between the vectors √ Z(t 1 )/Z (t 1 ) and √ Z(t 2 )/Z (t 2 ) (Bhattacharyya 1946).Another similar index also based on the relative abundances Z(t)/Z (t) is the proportion of similarity (Whittaker 1972): Both B Z (t 1 , t 2 ) and P Z (t 1 , t 2 ) range between 0 (no compositional similarity between t 1 and t 2 ) and 1 (perfect similarity).The Hellinger distance H Z (t 1 , t 2 ) also ranges between 0 (absence of change in the proportion of types from t 1 to t 2 ) and 1 (complete turnover of types).These indices are all defined conditional on Z (max(t 1 , t 2 )) > 0.
See Sect.7.1 for a proof.Proposition 2 shows that the Bray-Curtis index of the compositional beta diversity of a clonal family not becoming extinct converges over time to a finite constant in the admissible range [0, 1], depending on δ.Convergence to 1 as δ → ∞ occurs exponentially fast with δ and reflects the continual expansion of the clonal family as δ increases.Interestingly, asides from δ, the limit in t depends solely on α.Thus, as t → ∞, the Bray-Curtis index becomes insensitive to any aspect of the offspring and lifespan distributions beyond those captured by the Malthusian parameter, and to the number of types, K .It partitions the family of supercritical multitype branching processes into equivalence classes indexed by α and including processes with identical Bray-Curtis diversity.Hence, two non-extinct multitype branching populations with identical Bray-Curtis indices have similar Malthusian parameter α, whereas when their Bray-Curtis indices differ, the one with the highest index is expected to have a larger Malthusian parameter because bc Z (δ, α) increases with α.It follows easily from Proposition 2 that for every δ ≥ 0 if Assumption (A4) holds and conditional on non-extinction.
The indices H Z (t, t + δ), B Z (t, t + δ), and P Z (t, t + δ) tend to 0 or 1 as t → 0 regardless of δ.The behavior of these indices reflects intra-clonal stabilization of the frequency of the types as t increases.See (Jagers and Nerman 1996) for further discussion on the asymptotic composition of supercritical multitype branching processes.
Indices of beta diversity may also be defined at the GC level; for example, and The asymptotic behavior of GC-specific indices must be treated on a case-by-case basis.For the model studied in Sect.5, H Y (t, t + δ) converges a.s. to 0 as δ → ∞ when the branching mechanism is the same across clonal families.

Background and notation
Germinal center reactions begin polyclonal and tend to finish oligoclonal (Kroese et al. 1987;Küppers et al. 1993;Faro and Or-Guil 2013;Tas et al. 2016).The cross-sectional experiments that established these findings did not capture the temporal trajectory of individual clones.Thus, for example, we do not know whether clones that dominate a GC remain dominant over time.
To gain further insights into the dynamics of clonal dominance, we consider the infinite-type process introduced in Sect.3.5.2and formulate the overall immigration process (•) as a time-inhomogeneous Poisson process with rate r (•).The order of arrival induced by this process determines the types in the model.Since each type represents a specific clone, the sequence {T k } ∞ k=1 a.s.satisfies T k1 = T k and T k = ∞, = 2, 3 . ... The first identity indicates that the time at which the first member of the k-th clone joins the GC (T k ) is the time at which the first type-k founder B cell joins the GC.The second identity ensures that each clone has a unique arrival time in the GC.Consequently, the type-specific immigration processes { k (•)} ∞ k=1 are neither Poisson processes nor independent.The composition of the k-th clone at time t is For every k ∈ K = {1, 2 . ..}, the lifespan of every type-k B cell in the GC is a r.v. with c.d.f.G k (•), and their probability of division is p k 2 .Since binding affinity may differ between clones, we assume that these probabilities are independent and identically distributed (i.i.d.) r.v. with c.d.f.
Likewise, the c.d.f.G k (•) could be assumed to be randomly sampled from a family of distribution functions.Here, none of the types communicate and, conditional on { p k 2 , G k (•)}, we define the Malthusian parameter of the k-th clone as the root α k of the equation assuming such a root exists a.s.; it always does when p k 2 ≥ 0.5.When each G k (•) is an exponential distribution with parameter λ, the Malthusian parameters are given by α k = (2 p k 2 − 1)λ and form a sequence of i.i.d.r.v. with c.d.f.
The conditional probability of extinction of the k-th clone, given { p k 2 , G k (•)}, defined as the smallest non-negative root of the equation and Ney 1972).
Definition 2 (Clonal dominance) We say that the k-th clone dominates the GC reaction at time t if its size is equal to or exceeds that of any other clone at time t; that is, if We study temporal evolution of clonal dominance using , the fraction of B cells in the GC at time t belonging to the k-th clone.We set F k (t) = 0 when Y (t) = 0. We study evolution of the alpha and beta diversities of clonal families within GC using indices introduced in Sect. 4. Following initiation of a GC, species (clonal) richness is given by Assuming independence of (t) and {Z (k) 1 (t) = 0}.Species richness does not account for the abundance of non-extinct clonal families and grows indefinitely with t when Q > 0. In particular, when r (t) = r 0 , t ≥ 0, it can be shown that E(D 0 Y (t)) ∼ r 0 (1 − Q)t, such that increasing in the probability of extinction of clonal families (and of cell death) will proportionally reduce the number of clonal families in the GC.The Simpson index (Simpson 1949) which accounts for the abundance of clonal families and specializes to is the probability that two cells sampled with replacement at time t from a given GC belong to a same clonal family.The inverse of the Simpson index D 2 Y (t) = S Y (t) −1 provides a measure of the alpha clonal diversity: the higher S(t) −1 , the larger the clonal diversity, accounting for evenness (Rempala and Seweryn 2013).The Bray-Curtis index and proportion of similarity comparing clonal composition in a GC at times t 1 and t Here, we ask what these indices reveal about clonal dynamics in the absence and presence of differential fitness between clones, as previously considered (Tas et al. 2016).We note that the clone-specific diversity indices are D 2 Z (t) = 0, S Z (t) = 1, and B Z (t, t + δ) = 0 for every t, δ ≥ 0, because each clone includes a single type.

Clonal dynamics in the absence of differential fitness
In this section, we assume absence of differential fitness between clones.
Proposition 3 Suppose that all clones have identical offspring and lifespan distributions.Let α denote their common Malthusian parameter.Suppose that α > 0, Assumption (A4) holds, and (A5) there exists a positive, increasing function f See Sect.7.2 for a proof.Proposition 3 shows that, in the absence of differential fitness between clones, clonal dominance is a stable property unlikely to be challenged by other clones once established.Since the relative abundance F k (t) is asymptotically proportional to e −αT k a.s.
−→ 0 (k → ∞) and {W k } ∞ k=1 is an i.i.d.sequence, clones that arrive early in the GC are the most likely to become dominant.The relative size of the dominant clone tends to increase with α and decrease with the immigration rate (causing T k − T k−1 to stochastically decrease).Although the smaller α, the longer it may take for clonal dominance to establish, the replacement of a dominant clone by another clone joining the reaction later is unlikely to happen because this would require a sequence of independent r.v. with a decreasing trend to break a record value that increases with k: suppose the k * -th clone is dominant and W k * > 0; then, informally, in order for the k-th clone (k > k * ) to become dominant, we must have log ; in particular, if the immigration process is time-homogeneous, the increment between log W k * and log W k must grow linearly with k − k * , requiring the distribution of log W k (given W k > 0) to be heavy tailed, which is not possible since W k has finite moments (Harris 1963).Convergence of Simpson's index to a strictly positive r.v.reflects stabilization of the composition of GC in terms of clonal families and that the established order and composition are GC-specific.Also, lim t→∞ S Y (t) < 1 shows that GC do not converge toward monoclonality.The higher α, the higher S Y (t) (i.e., the more concentrated the distribution of clone size).Convergence of the Bray-Curtis index to the same limit as for Z(•) established in Proposition 2 arises from the fact that all clones obey the same offspring and lifetime distributions.Convergence of other beta diversity indices (to 0 or 1) captures absence of detectable temporal variation in clonal dominance.

Differential fitness induces unstable clonal dominance
We now assume that {α k } ∞ k=1 is a sequence of non-degenerate r.v., and show that any dominant clone is eventually outcompeted.
Proposition 4 Suppose that Assumption (A4) holds, See Sect.7.3 for a proof.Proposition 4 states that (i) the relative size of every clone within a GC becomes eventually negligible and (ii) every dominant clone is eventually outcompeted with probability one and regularly replaced as long as the GC reaction continues: an average of at least 3 weeks in lymph nodes and spleens, and much longer in Peyer's patches.The replacement of dominant clones contributes to the diversification of the BCR repertoire by allowing new clones to expand over time.Informally, the replacement of a dominant clone occurs primarily when a new clone that joins the GC is not on a path to extinction and has a Malthusian parameter higher than that of the dominant clone.The phenomenon is also driven by a record process; however, unlike the one describes in Sect.5.2, here the record process associates with the Malthusian sequence: α(t) = max{α k , k ∈ {1 . . ., (t)} : W k > 0}, t ≥ 0. Let L(t) denote the number of such record values between times 0 and t.Since the number of clones joining the GC between the L(t)-th and (L(t)+1)-th records tends to increase exponentially with L(t) (Neuts 1967), the rate of clonal turnover should rapidly slow down with L(t), unless, for instance, the Malthusian parameters of founder B cells visiting the GC increase stochastically over time.This scenario could potentially occur if the GC was to recruit B cells released by other concurrent GC, exhibiting a high affinity for any of the antigens driving the reaction.

Simulations: impact of the shape of immigration
We assumed that every G k (•) is an exponential distribution with mean 1 so each unit of time represented the average duration of one mitotic cycle.For reference, the mean mitotic cycle of activated B cells has been estimated to last 8-10 h (Chan et al. 2021).We also assumed that p k 0 was uniformly distributed in the interval [0.35, 1].Thus, each process Z (k) k (•) was super-critical with probability 0.23 (= 0.15/0.65).We were interested in studying the impact of the shape of the immigration rate on the dynamics of clonal families, considering three scenarios with different immigration rates: (1) r (t) = r 0 t −1/2 , where immigration intensity decreases at a rate of −1/2 as discussed in Sect.3.3; (2) r (t) = r 0 , with a constant rate over time; and (3) r (t) = r 0 t 1/2 , where immigration intensity increases at a rate of γ = 1/2.The immigration r (t) has not yet been estimated based on experimental data.However, in the stationary k (t) and F k (t)) within a GC (each line represents one clonal family).Bottom: Evolution of clonal diversity within the GC, as measured by the Inverse Simpson index, in three scenarios corresponding to r (t) = r 0 t −1/2 (left), r (t) = r 0 (center), and r (t) = r 0 t 1/2 (right), with r 0 = 10 in all three cases.See Sect.5.2 for detail case, it seems reasonable to assume that it is at least one order of magnitude faster than the mean mitotic cycle duration.Consequently, we set r 0 = 10 in all three cases.For computational reasons, we simulated the process over 20 units of time.Assuming that mitotic cycles and GC reactions last on average 10 h and 3 weeks, respectively, a duration of 20 units of time would cover almost half of the lifespan of a GC.We simulated the process 10 times in each scenario to identify patterns in the temporal trajectories of the Simpson index.
Results from the simulations are presented in Fig. 7.The top two panels show the evolution of the absolute and relative size of individual clonal families.Over the 20 time units, three clonal families reached clonal dominance (right panel).The bottom three panels indicate that clonal diversity, as measured by the Simpson diversity index, increases initially before leveling off between 5 and 10 units of time after initiation of the reaction.This pattern holds in the three immigration scenarios.In scenario (1) where the immigration rate decreases, the Simpson index tends to gradually decrease thereafter, suggesting that clonal diversity consistently reduces as the influx of founder B cells slows down.In the other two scenarios, diversity appears to also decrease but at a slower rate or less consistently among simulations.These results illustrate the potential benefit in sustaining the influx of B cells into GC in order to increase clonal diversity.

Biological relevance
Recent studies estimated that around a few hundreds of founder B cells colonize GC, clearing the field from a commonly held belief that GC are seeded by 2-10 B cells (Faro and Or-Guil 2013;Tas et al. 2016;Kroese et al. 1987).Our model leads to a similar conclusion and supports the hypothesis that oligoclonality results from differential binding affinity between clones.We note that past studies were based on cross-sectional data.They could not track the number of clones accumulated in GC over their entire lifespan and may not have accounted for those that were outcompeted and eliminated before the GC were observed and for clones that had not yet joined the GC.Therefore, reported estimates may represent lower bounds for the number of founder B cells that join GC, and their actual number could be larger.Proposition 4 also suggests that clones with the highest binding affinity will tend to arrive late in GC, such that clonal dominance at a particular time point need not reflect optimal binding affinity.

Steering the repertoire via sequential immunization
Antibodies play a key role in protective immunity against various pathogens.Although only a handful of people are known to have been cured from HIV thus far, numerous studies have demonstrated that the bnAbs that emerge in a significant proportion of individuals with HIV infection can neutralize a wide array of circulating viral strains by targeting conserved HIV epitopes, the specific sites on the antigen that these bnAbs bind to, marking a promising avenue for potential therapeutic interventions.Particularly, a successful vaccine regimen might be designed by recapitulating the process giving rise to bnAbs via a sequence of strategically designed prime-boost immunogens that can trigger the activation of B cells with the molecular and genetic traits of targeted bnAbs and guide a subset of the B-cell repertoire toward neutralization against HIV by selecting appropriate somatic hypermutations.Figure 1 provides an illustration of the principle of this strategy assuming two boosts; a successful vaccine might require a larger number of boosts.
Two obstacles to the elicitation of bnAbs through vaccination are the exceptionally high mutation frequency characteristic of HIV bnAb which vaccines may have to induce and the difficulty of the germline precursors of HIV bnAbs to bind with the envelopes of circulating viral strains.In addition, some mutations necessary for neutralization may be rare, and the process of somatic hypermutation can either lead to mutations at positions where the bnAbs align with the germline or erase previously acquired mutations conducive to bnAb development.Neither of these outcomes may be desirable as they can impede or prevent maturation towards neutralization.
To quantify the process of steering BCR towards bnAbs, we introduce a model that describes the accumulation of somatic hypermutations in BCR relative to their germline and the targeted bnAb.To fix ideas, our model focuses on the region encoded by the V gene, also known as the V segment.

Foreword to the modeling approach
In HIV vaccine research, evolution of the B-cell repertoire may be assessed by sequencing the heavy and light chains of thousands of antigen-specific B cells and comparing the mutations accumulated in their Ig gene loci with those found in relevant bnAbs.Multiple metrics may be designed to quantify this evolution in individual BCR sequences, depending in part on the class of B cells that the vaccine is expected to induce.Here, we opt for a simple, generic metric that enumerates mutations in five distinct regions of the heavy chain V segment.These regions are defined by whether mutations occur at positions where the aligned bnAb and germline gene sequences match or do not match.The proposed metric can be restricted to specific positions of the V segment without altering the model structure and interpretation.
Acknowledging the intricate nature of somatic hypermutation and antigen-mediated selection, we adopt a high-level modeling approach: instead of detailing BCR maturation on a per-position and per-nucleotide basis, we focus on counting the number of amino acid substitutions within the five identified regions.On one hand, aggregating multiple positions and nucleotide substitutions obscures the ability to specify how the accumulation of mutations alters antigen-mediated selection or the fitness of B cells.On the other hand, aggregation is expected to average out the impact of mutations on selection of B cells by the vaccine, at the level at which assumptions are postulated.
Our model offers a useful framework to understand some aspects of the evolution of the BCR repertoire.We use it in Sect.6.7 to interpret BCR sequencing data.The analysis underscores the potential challenges in precisely replicating a particular antibody.Regardless, achieving neutralization through vaccination might not necessarily require reproducing a specific bnAb, as other unidentified antibodies could also induce broad neutralization.

Partitioning the V segment into 5 classes of positions
Let B = (B 1 , . . ., B L ) and G = (G 1 , . . ., G L ) denote the amino acid sequence of the V segment of a given targeted antibody and its germline V gene.To fix ideas, suppose that the antibody is VRC01; its heavy chain is encoded by the VH1-2*02 allele (Fig. 8).Let V denote the set of amino acid positions of the bnAb V segment.This set may be partitioned into two subsets: the positions at which the germline and bnAb amino acids are identical and those at which they differ.We respectively denote these two subsets by V 1 and V 2 ; they satisfy denote the length of the V segment, and put L a = |V a |, a = 1, 2. For VRC01, L = 98, L 1 = 57, and L 2 = 41.Now, let S = (S 1 , . . ., S L ) denote the amino acid sequence of any BCR, to be compared against both its germline V gene and the bnAb.Let V 11 (S) = {p ∈ V 1 : S p = G p } denote the set of positions within V 1 at which S matches the germline and the bnAb.Define also ).The mutations accumulated in S that moved the BCR away from the bnAb are identified by V 12 (S), those shared with the bnAb belong to V 22 (S), whereas those that moved the BCR sequence S away from the germline V gene without changing its Hamming distance to the bnAb are in V 23 (S).See Fig. 8 for an illustrative example.Mutations may change the size of the subsets V ab (S); however, by construction, the subsets V 2b (S), b = 1, 2, 3 are invariant to mutations in V 1 ; likewise, mutations in V 2 do not change the subsets V 11 (S) and V 12 (S). Define (2, 1), (2, 2), (2, 3)}.The total number of mutations in S is N (S) = L 12 (S)+ L 22 (S)+ L 23 (S).A successful sequential immunization regimen might elicit BCR getting closer to targeted bnAbs as they accumulate mutations.A simple measure of the immune response toward this goal is the net mutational gain/loss relative to a bnAb, defined as (S) = L 22 (S) − L 12 (S), with a positive value indicating overall progression of the sequence toward the target antibody.The dynamic model formulated in the next three sections describes the evolution of L(S) when S is sampled from a tree generated by a multitype branching process.

The mutational model, absent of cell kinetics
Let S (0) = (S (0) 1 , . . ., S (0) L ) denote the amino acid sequence of the V segment of the BCR of any B cell at birth.We make the simplifying assumption that the two daughters of a dividing B cell inherit identical Ig gene sequences from their mother at division, represented by S (1) = (S (1)  1 , . . ., S (1) L ).Unless a synonymous mutation occurs, the sequences S (0) and S (1) differ at the mutated position.Let π denote the probability that mutations occur in the V segment.Allen et al. (1987) estimated that the BCR of one out of every other B cell inherits one mutation from their mother (Allen et al. 1987), such that we may set π ∈ (0; 0.5) when focusing on the V segment.In what follows, we formulate a mutational model that describes the accumulation of somatic hypermutations in the subsets V ab , assuming at most one mutation happens per cell cycle.
In Case 2.1, a non-synonymous mutation in V 21 (S (0) ) matches the bnAb amino acid at the mutating position, bringing the BCR one amino acid closer to the bnAb.In Case 2.2, a non-synonymous mutation also in V 21 (S (0) ) does not match the bnAb amino acid at the mutating position and grows V 23 (S (0) ) by one position.In Case 2.3, a non-synonymous mutation in V 22 (S (0) ) reverts the BCR to the germline amino acid.In Case 2.4, a non-synonymous mutation occurs in V 22 (S (0) ) which does not match the germline amino acid.In Case 2.5, a non-synonymous mutation in V 23 (S (0) ) reverts the BCR back to the germline amino acid.In Case 2.6, a non-synonymous mutation also in V 23 (S (0) ) matches the bnAb amino acid.In Case 2.7, a mutation, possibly synonymous, in V 2 matches neither the germline nor the bnAb amino acid, leaving the subsets V 21 (S (0) ), V 22 (S (0) ), and V 23 (S (0) ) unchanged.We note that mutations within V 21 (S (0) ) and V 22 (S ( 0) ) must necessarily be synonymous, whereas those within V 23 (S (0) ) can be either synonymous or non-synonymous.

A parametrization of the transition probabilities
We present a parametrization of the probabilities π S (0) i and π S (0)  i|a structured around the number of positions in V ab (S (0) ) to facilitate their interpretation and selecting their values in simulations.We note that, depending on S (0) , not all 10 transitions presented in Sect.6.4 are possible; for instance, Case 1.1 cannot happen when S (0) 1 = 0.These cases are handled by setting π S (0)  i|a = 0.The probability of a mutation occurring in V 1 compared to V 2 may depend on the relative sizes of these two sets, and we write π S (0) 2 for some constant ρ 1 ∈ (0, L L 1 ) which represents the selection bias in favor of (ρ 1 ∈ (1, L 1 /L)) or against (ρ 1 ∈ (0, 1)) bnAb-like mutations.When ρ 1 = 1, the odds of a mutation occurring in L 2 and depends solely on the size of V 1 relative to that of V 2 .Since mutation rates vary across positions along the Ig gene sequences, positions in V 1 may collectively have a higher (or lower) probability of mutation than those in V 2 , and we can adjust ρ 1 ∈ (1, L L 1 ) to increase the likelihood of the mutation occurring in V 1 relative to V 2 , and vice versa when ρ 1 ∈ (0, 1).

A branching process of convergent evolution
In a clone started from a naive B cell, the BCR sequence of the founder B cell is such that L(S) = (L 1 , 0, L 2 , 0, 0) because it is unmutated.As cells divide, mutations accumulate in the clonal family, and BCR sequences S accumulate mutations, and the collection of vectors L(S) represents the position of the clonal family with respect to the germline and bnAb sequences.Here we investigate whether the BCR repertoire converges to a bnAb, and whether convergence to a bnAb requires a highly favorable set of circumstances.
To describe the intra-clonal evolution of the distances L(S), we define a multitype age-dependent branching process in which every cell is assigned a type represented by a vector that takes values in the set where N = {0, 1 . ..} represents the set of non-negative integers.The entries of are interpreted just as those of L(S), as described in Sect.6.4.
For every (0) , ( 1) , (2) ∈ K × K × K, the conditional probability that any cell of type (0) divides into one cell of type (1) and one cell of type (2) , given it divides, is assumed to be q The first equation specifies the probability of no mutation occurring.The second one specifies the probability of a mutation occurring in V a according to Case a.i, as described in Sect.6.4, with both daughters inheriting the same set of mutations from their mother.The third equation reflects the assumption that both daughter cells cannot inherit different mutations from their mother.We assume that the probability of death and the lifespan distribution of every cell are independent of its type: p 0 = p 0 and G (•) ≡ G(•), ∈ K.This assumption is motivated by the aggregation of positions and amino acid substitutions in defining the five regions, as discussed in Sect.6.2.Together with the conditional probabilities q (0) (1) (1) , they model the combined outcome of mutation and antigen-mediated selection.Because it does not distinguish these two processes, the proposed model is partly phenomenological.The entries of M = (m (0) ( 1) ) ( 0) , ( 1) ∈K are given by m (1) ( 1) ) ( 0) , (1) ∈K .
See Sect.7.4 for a proof.The irreducibility assumption on Q holds, for example, if mutations in V 1 and V 2 are reversible.It follows from Proposition 5 that, conditional on survival (i.e., W > 0), the frequency of mutations within a clone converges a.s.over time; specifically, such that the entries of the vector v/K provide the asymptotic mean frequency of the number of positions of V falling into the subsets V ab , (a, b) ∈ {(1, 1), (1, 2), (2, 1), (2, 2), (2, 3)}, in expanded clones.

Drift of the BCR repertoire relative to HIV bnAbs
We randomly selected 1000 BCR sequences sampled from human subjects, none of which were living with HIV (Guo et al. 2019).All BCR were inferred to be encoded by the VH1-2*02 allele.We calculated the number of mutations relative to VH1-2*02 and the net mutational gain/loss against four HIV bnAbs: VRC01, 3BNC117, IOMA, DH270.1, all originating from VH1-2*02 (Fig. 8, middle row).As expected, BCR induced by antigens unrelated to HIV preferentially accumulate mutations that are distinct from those observed in the HIV bnAbs, as evidenced by the amount of mutational losses increasing with the number of mutations.Thus, sequential HIV vaccines may have to select for rare mutations to steer antibodies toward bnAbs.The data also suggest that inducing bnAb-specific mutations is more challenging for some bnAbs than other, as indicated by the faster drift of some BCRs away from IOMA and DH270.1 relative to 3BNC117 and VRC01.
To interpret this observation, we used our model to first examine the evolution of the BCR repertoire relative to VRC01.We set the probability of death ( p 0 ) to 4 10 to ensure the survival of clonal families while preventing excessive expansion for computational feasibility.The V segment of heavy chains accounts for slightly less than half of the variable domain, but is somewhat more prone to somatic hypermutations than other parts of the BCR.Therefore, taking Allen et al. (1987)'s estimate of the mutation rate into consideration, we set the probability of mutation in the V segment ( π) to 1 4 .The mutation rate may differ between V 1 and V 2 , partly because some bnAbspecific mutations may be much less frequent than non-bnAb specific mutations.These probabilities also depend on the antigens to which B cells are exposed.Thus, we considered several values of ρ 1 = 1 10 , 1 4 , 1 2 , 3 4 , 9 10 to perform a sensitivity analysis.The BCR sequencing data used in our analysis resulted from exposures to highly diverse antigens.Therefore, we do not expect the per-position mutation rates to differ between the sets V ab (S (0) ), and put ρ 1|1 = ρ 1|2 = ρ 2|2 = 1.As argued in Sect.6.5, we set σ 1|1 = 1 3 to specify the likelihood of a mutation in V 11 (S ( 0) ) being synonymous.To specify the probability of a mutation in V 12 (S (0) ) not reverting to a germline amino Fig. 9 Mutational paths toward a targeted antibody sequence.Each horizontal block of adjacent vertical bars represents one antibody sequence, partitioned into nucleotide positions.Arrows represent one-nucleotide mutations between antibody sequences.This diagram gives a simplified overview of the sets of mutational paths induced by somatic hypermutations started from an unmutated (germline) antibody sequence toward a targeted antibody sequence.Progression toward the targeted antibody may be hampered by the fact that as sequences accumulate targeted mutations, the number of candidate unfavorable mutations increases relative to those that are favorable.When the antibody is one mutation away from the target, the remaining mutation is up against all other potential mutations.This challenge may be collectively addressed by the multiple B cells from a clonal family that may express identical or similar antibody sequences, each offering a chance of successfully mutating toward the targeted sequence acid, we set σ 2|1 = 1 10 .This choice is based on several considerations.Firstly, in the best-case scenario where reverting to a germline amino acid requires a single mutation, most substitutions at the wobble base, which may account for approximately one-third of all mutations, would not qualify.Among the remaining two-thirds of potential mutations at the first and second bases, only those that induce the specific germline amino acid would produce a reverse mutation, acknowledging however that the (effective) average number of candidate amino acid mutations is likely lower than 19 (ignoring mutations into stop codons).Furthermore, as mutations accumulate, some positions may require two or even more simultaneous mutations to revert to the germline amino acid.With an effective average number of five amino acid mutations, we might expect σ 2|1 2 15 in the best-case scenario; accounting for cases requiring more than one mutations provides support for setting σ 2|1 = 1 10 or lower.We note that the choice of σ 2|1 has limited impact if L 12 (S (0) )/L 1 is small.Following a similar line of arguments, we set σ 12|1 = 1 10 to specify the probability of a mutation in V 12 (S (0) ) matching the bnAb amino acid, acknowledging that a lower value may be more appropriate to study drift of the B cell repertoire following exposure to HIVrelated antigens which may induce bnAb-like mutations with higher probability.Thus, we explored values of σ 12|1 = 1 5 and σ 12|1 = 1 2 in sensitivity analyses.Based on previous arguments regarding synonymous mutations, we set P S S (0)  21 |M S (0) 21 = 1 3 .Hence, P S S (0)  22 |M S (0)

21
= 2 3 which gives σ 13|2 = 2 − σ 12|2 after algebraic calculations, and we set σ 13|2 = 19 10 .Justified by a similar line of arguments, we set σ 21|2 = 1 10 and σ 23|2 = 19 10 .Finally, we assume that σ 31|2 = σ 32|2 = 1 10 using the line of arguments used to justify setting σ 12|2 = 1 10 .Results from simulations are displayed in Fig. 8, sorted by values of ρ 1 (columns) and σ 12|1 (rows).When σ 12|1 = 1 10 , overall, patterns observed in simulations are similar to those from BCR sequencing data, indicating a tendency of the BCR repertoire to drift away from the bnAb as mutations accumulate in the V segment, unless ρ 1 is close to 0; for example, ρ 1 ≤ 1 10 , reflecting that the per-position mutation rate is at least 10 times higher in V 2 than in V 1 .In the absence of selection bias (ρ 1 = 1), mutations would appear in V 2 with a probability of 41.8% (= 41 98 × 100%).Thus, even with almost 50% chance of inducing a mutation in V 2 , the BCR repertoire is expected to rapidly drift away from VRC01 (see panel with ρ = 1.0) because the bnAb amino acid is only one of the candidate mutations.Our simulations suggest a value of ρ 1 0.5 for VRC01, indicating that affinity maturation during exposures to HIV-unrelated antigens has a preferential bias for inducing mutations at positions at which VRC01 and VH1-2*02 do not match, and that the probability of this happening is about twice as large compared to if there was no selection bias (i.e., case where ρ 1 = 1).The results obtained by increasing σ 12|1 to 1 2 or 4 5 suggest that this parameter has some but relatively limited influence on the results, with its impact decreasing with ρ.
The value of ρ 1 appeared to differ between the four bnAbs.After accounting for the number of mutations in their respective V segment, simulations suggested the following values of ρ 1 : ρ 1 0.5 for 3BNC117 and ρ 1 ≥ 0.9 for both IOMA and DH270.1.Thus, induction of mutations in V 2 were less likely for IOMA and DH270.1 than for VRC01 and 3BNC117.One potential explanation is the fact that VRC01 and 3BNC117 are more mutated than IOMA and DH270.1, making the induction of VRC01-and 3BNC117-like mutations easier.Typically, as BCRs approach a targeted antibody, the challenge of accumulating additional advantageous mutations tends to increase because the probability of encountering unfavorable mutations eventually surpasses that of favorable ones, as illustrated in Fig. 9.This could result in a slowdown in progress as B cells approach the target.The point of equilibrium towards which members of a clonal family tend to gravitate is formally described in Proposition 5 and equation (23).

Proof of Proposition 2
We begin the proof with the Bray-Curtis index.Since Assumption (A4) holds and W > 0, we have for δ ≥ 0 and t large enough that ).
The first term in the right-hand side converges a.s. to P{W 1 > 0} by a Strong Law of Large Numbers for i.i.d.r.v..For the second one, notice that Since P{W k (t + δ − T k ) > 0 (t − f (t)) = } − P{W k > 0} is a monotone function in t for every k, , we can exchange limit and summation to obtain We that, for every δ ≥ 0, Proceeding similarly, we show that the right-hand side of (26) normalized by (t) also converges in probability to P{W 1 > 0}.Therefore, for every δ ≥ 0, a(t, δ)/ (t) P −→ P{W 1 > 0} as t → 0. We finally deduce from (25) that J Y (t, t + δ) P −→ 0 as t → ∞.Next, to prove the limit of BC Y (t, t + δ), we write = 1, To establish the limit of the proportion of similarity P Y (t, t + δ), we have that lim t→∞ P Y (t, t + δ) and the result follows, completing the proof.

Fig. 2
Fig.2Dynamics of germinal center B cells as captured by the proposed model.Founder B cells join the GC at random time points (T 1 , T 2 , T 3 , T 4 . ..).These arrival times are described by a point process, potentially time-inhomogeneous.Each founder B cell initiates a clonal family in which B cells either divide or die.Each cell is assigned a type at birth; when it divides it may produce offspring of different types (e.g., reflecting the impact of somatic hypermutations).In the example depicted here, the population includes two types of B cells, but more complex models with arbitrary number of types may be constructed

Fig. 4
Fig. 4 From left to right: median, mean, and ratio of medians (or means) of the time to first type-k founder joining the GC when the time-inhomogeneous Poisson process has local intensity r k (t) = g k r 0 t γ (r 0 = 1).These parameters are plotted as a function of γ ∈ (−1, 1)

Fig. 6
Fig. 6 Behavior of a two-type process, with and without immigration.Panels A1-A4 show results from 10 independent simulations of the model without immigration, with (m 11 , m 12 , m 21 , m 22 ) = (0.20, 0.50, 0.55, 0.20) and G(•) assumed exponential with parameter λ = 1/24.Each trajectory represents the evolution of one clone.Each of them begins with one type-2 (∼ high binding affinity) cell.Panel A1 (resp., A2): number of low (resp., high) binding affinity B cells per clone over time; Panel A3: count of high versus low binding affinity B cells which grow at the same rate; Panel A4: ratio of the numbers of low and high binding affinity B cells plotted over time; the black straight line indicates its a.s.limit 0.37.Panels B1-B4 are identical to panels A1-A4, except that simulations were conducted by setting (m 11 , m 12 , m 21 , m 22 ) = (0, 0.30, 0.55, 0.20).The a.s.limit is 0.31 in this setting.Panels C1-C4 are similar to those above, but the model now includes immigration.The parameters of the branching process were identical to those used for panels B1-B4; the immigration Poisson process was time-homogeneous and assumed that g 1 = g 2 = 0.5 and r 0 = 0.25

Fig. 7
Fig. 7 Temporal evolution of clonal dominance and clonal diversity.Top: Evolution of the absolute and relative size of individual clonal families (i.e., Z (k) 23|2 .Finally, we assume that P S S(0)  21 |M S (0) the bnAb sequence.Since P S S (0) 21 |M S