Modeling and analysis of the meta-population dynamics of lymphocyte repertoires

Journal of Computational and Applied Mathematics 184 (2005) 223 – 241 www.elsevier.com/locate/cam

Modeling and analysis of the meta-population dynamics of lymphocyte repertoires Ramit Mehr∗ Faculty of Life Sciences, Bar-Ilan University, Ramat-Gan 52900, Israel Received 22 February 2004; received in revised form 10 July 2004

Abstract The vertebrate adaptive immune system, with its learning and memory capabilities, is a particularly rich example of the complexity of biological systems. The most difficult challenge in immunological research is the understanding of lymphocyte repertoires—the generation of their diversity and the forces that shape the ever-changing dynamics of lymphocyte clones. The explosively growing quantities of immunological information, which are rapidly becoming available from experimental and clinical studies, necessitate the use of theoretical tools, namely, mathematical and computational modeling and bioinformatical methods of analysis. In this paper, I review several examples of problems in lymphocyte repertoire modeling and analysis, demonstrate the types of solutions employed, and highlight the contribution of these theoretical studies to immunological research. © 2005 Elsevier B.V. All rights reserved. Keywords: Lymphocyte repertoires; Lymphocyte development; Mathematical model; Computer simulation; Bioinformatics; Sequence analysis; Lineage tree analysis

1. Introduction One of the most intriguing challenges to theoretical biologists is presented by the adaptive immune system—one of the only two biological systems capable of continuously learning and memorizing its experiences. Within this rich source of theoretical problems, one particularly difficult area is the study of lymphocyte repertoires—the generation of antigen receptor repertoire diversity, the dynamics of lymphocyte development under normal or immune-deficient conditions, the forces of selection and interaction ∗ Tel.: +972 3 531 7990; fax: +972 3 535 1824.

E-mail address: [email protected]. 0377-0427/$ - see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2004.07.033

224

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

with antigen that shape lymphocyte repertoires. The difficulty in modeling lymphocyte repertoires does not result from the non-linear behavior of each component in the system, nor from the astronomic numbers of lymphocyte clones in the human patient or the experimental animal. Even though actual modeling of comparable number of clones would require a much higher computing power than is available to the average theoretician, studies on a smaller number of clones can often give sufficient insight into the behavior of the whole system. The main difficulty stems from the need to formulate a description of the system on several levels: genetic, molecular, cellular and systemic. The various aspects of immune system cell behavior are not tractable by the classical mathematical models of cell populations which have traditionally been used. For example, when studying the selection of developing T lymphocytes in the thymus, there is only so much one can do with a model that is limited to the level of cell populations. Before long, the rearrangement of T cell receptor (TCR) genes, and the interaction of their protein product with major histocompatibility (MHC)-peptide complexes, must be taken into account. Hence it is not only the population dynamics within each clone (i.e., cell division, differentiation and death processes) that must be modeled, but also the meta-dynamics: creation, selection and elimination of whole clones in the population, and the molecular interactions or genetic changes that form the basis for these meta-dynamics. The solution to the demand for increasing complexity depends on the problem at hand. No model can integrate all aspects of the real system—nor should it attempt to do so, because models must be simple enough to be tractable. Hence it is up to the researcher to decide which aspects of the system to focus on, and which aspects to neglect, in each part of the study. A modeler may choose to construct a simulation of a number of individual cells, and follow them in time, as we have done in our studies of isotype switch [48], but this is feasible only if the processes are very simple and the number of cells that must be simulated in order to receive statistically significant results is not large. Alternative solutions may be to follow one cell at a time, or follow whole clones of cells rather than individual cells in time, as demonstrated below. In any case, computational challenges abound, not only in formulating models and running simulations thereof, but also in creating novel methods of analysis of the exponentially growing amounts of data generated daily by both models and experiments.

2. Combining genetic and cellular processes: antigen receptor gene rearrangement How does one deal with a problem requiring a simulation of large numbers of lymphocyte clones, with many parameters linked to each clone? We shall demonstrate this using the problem of randomness versus sequentiality in antigen receptor gene rearrangement. This series of studies addressed the apparent contradiction between continuous antigen receptor gene rearrangement in lymphocytes and allelic exclusion. In contrast to the paradigm of “one gene–one protein” which rules most of biology, the genes for the antigen receptors of B and T lymphocytes are rearranged from segments taken from separate libraries in the gene locus. Gene rearrangement is an irreversible process occurring during the development of a B or T lymphocyte. The B cell receptor (BCR) light chain variable region gene, and the T cell receptor (TCR) alpha chain variable region gene, are rearranged from two segments, V (variable) and J (joining), taken from two consecutive libraries in the BCR light chain or the TCR alpha chain gene locus. Similarly, the B cell receptor (BCR) heavy chain variable region gene, and the T cell receptor (TCR) beta chain variable region gene, are each rearranged from three segments, V, D (diversity) and J, taken from three consecutive libraries in the corresponding gene locus. Together with the fact that cutting and joining of

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

225

gene segments are imprecise, and that specialized enzymes also perform additional deletion or addition of random nucleotides at the joining points, this combinatorial generation of receptor genes results in an enormous diversity, giving a repertoire of receptors able to recognize every conceivable biological molecule. Rearrangement may occur on each of the two alleles of every receptor chain locus, until a productive rearrangement (which can be transcribed and translated into a protein) is achieved and the cell can express a complete and functional receptor. However, the long-accepted concept of allelic exclusion states that antigen receptor chain genes in B and T lymphocytes are expressed from only one of the two alleles in each cell, ensuring that each lymphocyte expresses a single antigen receptor specificity on the cell surface. Recent evidence suggests, however, that in the BCR light chain [28], and in the TCR
chain [36], rearrangement may not stop after a productive receptor gene has been formed and expressed. Secondary light chain rearrangements occur not only after non-productive rearrangements, but also after productive rearrangements that render a B cell autoreactive (“receptor editing”) [37], resulting in one of the mechanisms of central tolerance. However, if the choice of allele for secondary rearrangements is random, it is (at least theoretically) possible that a cell will rearrange and then simultaneously express two different light chains. How then, if at all, is allelic exclusion maintained in the face of continued rearrangement? For
 T cells, the situation is even more complex [36]. Rearrangement and expression of TCR
chain genes does not stop after the expression of the first rearranged
chain, but appears to continue until the cell is either positively selected, or dies. Due to the lack of allelic exclusion in TCR
, a T cell may not only contain two productively rearranged TCR
alleles, but also simultaneously express the two resulting TCRs. Indeed, the frequency of T cells simultaneously expressing two different V
genes has been found [34] to vary between 10−3 –10−4 . Only the cell surface expression of three V
genes was monitored, which means that the frequency of T cells expressing any pair of V
genes may be orders of magnitude higher. Independently, 26% of various T cell clones were found to contain two productive V
–J
rearrangements [24,25]. This was an alarming discovery, because an
 T cell that matures in the thymus expressing two different TCRs may be positively selected on one of them, while the other TCR may be autoreactive. Why, then, is allelic exclusion not maintained in TCRs? The observations of allelic inclusion in TCR
raised several more questions. Can allelic inclusion be fully accounted for by multiple rearrangements alone? Do these rearrangements occur completely at random, or is there some underlying order? What is the role of positive and negative selection in driving, or limiting, the process of TCR gene rearrangement? Several previous models, which attempted to explain experimental observations such as the fractions of cells containing two productive TCR
rearrangements, did not sufficiently account for TCR gene organization, which limits secondary rearrangement, and for the effects of subsequent thymic selection. The work described in the following was designed to address these questions. 2.1. The simulations We first created a stochastic simulation of BCR gene rearrangement. The general structure of the simulation is given in Fig. 1. It simulates the following steps in a developing B cell’s life: • Cell “birth”—the simulation starts each new cell with all its light chain genes in the un-rearranged (germline) state.

226

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

Fig. 1. A schematic representation of our simulation of immunoglobulin light-chain gene rearrangement. Actions are depicted in rectangles, and decisions are depicted in diamonds. Pd is the death probability of cells that have not been subjected to negative selection, while Pdas is the death probability of “anti-self” cells that have been negatively selected (Pdas > Pd ).

• Light chain choice: the cell decides whether to rearrange  or  next. • BCR  rearrangement: choice of which allele to rearrange, and on this allele, which J  segment to use. J segments 5 to the chosen segment are deleted. • BCR  rearrangement—similar to  rearrangement, but without deletion. • Selection: probabilistic decision whether the rearrangement is in frame, whether the resulting light chain can pair with the existing heavy chain, and whether the resulting BCR is autoreactive (anti-self). Rearrangement is repeated until a productive, H/L matched, non-autoreactive BCR is produced, or until the cell dies. Cells containing an anti-self rearrangement are assigned a high death probability. Cells containing only out of frame V–J joins, H/L mismatched heavy-light chain pairs, or germline , are assigned a moderate death probability. If any of the latter types of cells do not die, they may attempt secondary rearrangements. Cells containing one in frame, H/L matched, non-autoreactive rearrangement, are allowed to mature. The parameters for the simulation, that is, probabilities of death in each stage, probabilities for isotype or allele choices, etc, were varied between simulation runs within the biologically plausible ranges.

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

227

For example, the probability for a V–J rearrangement to be productive is one-third, because this is the probability that V will be joined to J in the correct reading frame. On the other hand, the death rate of autoreactive cells is not known, and hence a whole range of values was studied in the simulations. We have then extended the simulation to include all variable region segments, V, (D), and J, for the heavy as well and the light chain [17]. The extended simulation allows for cell division, and follows all the cells in each lineage. The simulation of TCR gene rearrangement [36] was similar, with two extensions: (a) it records all variable region segments, V, D, and J, for the four possible chains—
, ,  and . (b) It has three types of selection: one at the pre-TCR stage checking that the  chain is functional, and then positive and negative selection of cells expressing an
 TCR. The description of the rearrangement status of each allele of each receptor chain, including the probability for further rearrangement associated with each individual segment, necessitate the use of a large number of parameters associated with each cell. Fortunately, we could reduce the amount of computation involved by reducing the description on the level of cell population dynamics to that of a single clone of cells, its divisions, death and so on. This is because the experimental data to which simulation results were compared consists of repertoire statistics, and such statistics could be collected by successive, rather than simultaneous, simulations of many clones. This was fortunate, because we found that we had to simulate very large numbers of cells (up to 105 in the BCR simulations) in order for the properties of the repertoire to be stable between simulations. In the TCR rearrangement simulations, with high cell division probabilities, an order of 104 individual cells may all belong to a small number of clones. Thus it was necessary to include a large number of independent clones (each possibly containing many cells) in each simulation, for parameters such as the ratio to stabilize. The number of clones was considered to be sufficient when both inter- and intra-simulation variabilities were small (less than 10% of the initial variability). This has been achieved for all quantities measured by generating 4000 independent clones in each simulation. The variability criterion thus defined helped to identify the proper division probabilities (probability of division following -selection and following positive selection of
 T cells) to be used in the simulations. A further simplification was that each event of selection was described by a single probability parameter, neglecting the details of the molecular interactions. 2.2. Results: BCR gene rearrangement The observed ratio between  light chain- and  light chain-bearing B cells in the murine serum is about 20:1; even in immature murine bone marrow cells it is larger than 10:1. It was controversial whether the :  ratio can be explained solely on the basis of the higher potential for multiple rearrangements in the  locus, combined with immature B cell death due to negative selection, without assuming preferential expansion of  B cells over  B cells. In addition to explaining the high :  ratio, the postulated ability of a B cell to make secondary rearrangements on a single  allele was suggested to explain why about 70% of mouse splenic  B cells have only one rearranged  locus, the other one remaining unrearranged [4]; and why in mice that have only one functional  locus,  B cell production is about 70% (rather than one-half) of that in wild-type mice [1]. Furthermore, it has been suggested [1,12,47] that rearrangement may proceed sequentially rather than stochastically. That is, that 5 J segments are used before 3 J segments. Hence, the following question arose: can ordered rearrangement account for the observations on allele bias, J usage, and the :  ratio? We tested the hypothesis of allele bias (preference to rearrange a rearranged allele first) as follows. Simulations were performed with either no bias (choice of the allele for the next rearrangement is

228

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

completely random) or strong bias (the cell continues rearranging the already-rearranged allele until there are no more J segments available for rearrangement on that allele, and only then it may switch to the other allele). The results were compared to experimental observations. Similarly, we tested the hypothesis of sequentiality in J usage, using simulations in which the choice of the next J segment for rearrangement was either random, strictly sequential (first J1, then J2, and so on) or quasi-sequential (an intermediate case with partial bias, that is, preference of J1 and J2 over J4 and J5). We have also checked the influence of parameters such as the probability of choosing  over , and the stringency of negative selection as expressed by death probabilities, on the resulting repertoires. The final conclusion of the series of studies was that BCR gene rearrangement is ordered at some levels but not others. First, our simulations support the notion that  light chain rearrangement precedes  light chain rearrangement in most, if not all, cases. On the other hand, our results did not support the hypothesis of allele preference. Third, our results show that experimental observations are consistent with the hypothesis of quasi-sequential rearrangement within each  allele. That is, when J1 and/or J2 are available for rearrangement, the cell is more likely to choose one of these segments over J4 or J5; but not with completely random segment choice. The effect of negative selection was determined by the death probability of cells that have not succeeded in producing a productively rearranged, non-autoreactive BCR. We found that the higher the death probability, the larger the :  ratio. Cells are allowed few rearrangement attempts and hence are likely to die before exhausting  and proceeding to  rearrangements. Our simulations showed that in order for the results to be consistent with experimental observations, we must assume that negative selection limits the number of rearrangement attempts to 2–3 per cell, in agreement with previous results [32]. Therefore, we refer to BCR rearrangement as a negative selection-limited process. Taken together, the above results reveal our proposed answer to the question of allelic exclusion in B cells. We propose that B cell allelic exclusion results from a certain degree of order in gene rearrangement and a stringent process of negative selection. The quasi-ordered use of isotypes and segments allows a B cell to maximize the number of rearrangements attempts. Because negative selection limits the number of rearrangement attempts to 2–3 per cell, it is extremely unlikely for two productive rearrangements to exist simultaneously on two light chain alleles. This results in the almost complete absence of cells expressing two different BCRs from the repertoire, i.e. in effective allelic exclusion. 2.3. Results on TCR gene rearrangement TCR gene rearrangement uses the same mechanisms as BCR gene rearrangement, and differs from it only in the number of segments of each type and in the selection processes operating on the developing cells. One of the measurable quantities which has received much attention in the literature is the
:  ratio in thymocytes and mature T cells, which can be 20:1 or even higher, depending on the tissue studied. Theoretical predictions based on models that do not include multiple rearrangements fell around 2:1, which is far from the experimentally observed range of values. The difference was attributed to cell division. Our simulations showed that both multiple rearrangements and cell divisions are required in order to explain this ratio, in contrast to the situation in B cells where there was no need to invoke cell divisions. The ratio increases with the number of cell divisions after  selection and after thymic selection. Conversely, the ratio decreases when we increase the death probabilities of cells that fail  selection or
 selection. Additionally, the
:  ratio increases with the probability that  rearrangement will precede  rearrangement.

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

229

The measure for the extent of allelic exclusion, or rather allelic inclusion, in T cells, is the fraction of TCR
“double-productives”—T cells that carry productive TCR
rearrangements on both alleles. Theoretical models predict that secondary rearrangements are necessary to explain the experimentally observed fraction of up to 26% TCR
“double-productives”. Our simulations examined the dependence of this quantity on the degree of order in TCR
rearrangement and on selection probabilities. Fractions of TCR
double-productives higher than 20%, as observed experimentally, are obtained in our simulations only when we allow multiple TCR
rearrangements but assume they are non-ordered and unbiased. At the least, if there is order in TCR
rearrangement, it is masked by the large number of rearrangements per cell (10–12 rearrangements in our simulations). A novel quantity defined in this study, for which no observations existed, is the fraction of cells with an auto-reactive receptor among the TCR
“double-productives”. This fraction is independent of multiple rearrangements, because it depends only on the last rearrangements on the two alleles. However, we found that the fraction of autoreactive “double-productives” is highly sensitive to the death rate of autoreactive thymocytes. If this rate is low, as it must be to get 26% “double-productives”, then the fraction of autoreactive double-productives can be as high as 70%. This value is only an upper bound, since it was obtained for the case in which the cell is selected only according to its last rearrangement. Otherwise, this number will be lower, and will also depend on the relative expression levels of the two receptors, which are not addressed by the current model. More experimental data would be beneficial for settling this issue, which may help elucidate instances of escape from central tolerance in T cells. In spite of the strong similarities revealed in our studies between the way rearrangement seems to operate in B and T cells, it is worthwhile to note a crucial difference between the development of T cells to that of B cells. While BCR rearrangement seems to be limited by negative selection only, T cell development, on the other hand, seems to be limited by positive rather than by negative selection: developing T cells in the thymus are allowed a much more generous time window for continued TCR
rearrangement, so that multiple rearrangements on both alleles become the rule rather than the exception. The probability that a cell will contain more than one productive rearrangement of the TCR
chain, and even express two TCR
chains simultaneously, is far from negligible. Moreover, positive selection may rescue the cell from death, and allow it to mature, based on the virtues of only one of its expressed receptors, as long as the other receptor is not so extremely auto-reactive as to cause immediate deletion of the cell. This is a potentially dangerous situation, because the second receptor may still be weakly autoreactive, or, worse, may be specific, with high affinity, to a self-peptide that is not presented in the thymus. In spite of the existence of peripheral mechanisms of self-tolerance which safeguard against improper activation of T cell, such improper activation does sometimes happen. Thus, understanding TCR gene rearrangement, selection and editing is key to understanding autoimmunity.

2.4. Molecular interactions and their effect on clonal dynamics: T cell selection In the above-described examples, we focused on the genetic processes (gene rearrangement or mutations) which drive population meta-dynamics mostly through the creation of new clones. On the other hand, selection processes, which decide the fate of the clones—survival or deletion—were condensed into a single parameter representation. In a different set of studies [5,6], we took the complementary approach, of focusing in detail on the molecular interactions underlying clonal selection (of thymocytes in this case) while neglecting the genetic details for the sake of simplicity.

230

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

The initial T cell repertoire, after its creation through gene rearrangement, is further shaped in the thymus by two selection processes. Positive selection promotes the differentiation to a further developmental stage of cells bearing receptors with a sufficiently large affinity for peptides presented on molecules of the major histocompatibility complex (MHC) expressed in the thymus. This confers to T cells the property of self MHC restriction: they recognize peptides presented in the groove of host MHC molecules, but ignore them when presented on foreign MHC. Then, negative selection deletes cells whose TCRs bind thymic MHC–peptide complexes (pMHC) with very high affinity, thus preventing the emergence of self reactive T cells. Overall, only about 3% of thymocytes have the intermediate affinity needed to fully mature. In addition to self MHC restriction, the mature repertoire is also characterized by a high alloreactivity. Typically, 1–24% of T cells react against the product of a given foreign MHC allele. This high response frequency is hard to reconcile with the fact that only one T cell in 104 –106 of the naive repertoire recognizes a given pathogen. What should the quantitative properties of the processes driving TCR generation and selection be in order to produce the experimentally observed levels of self restriction, alloreactivity, and antigen response? We have addressed this issue using a computer simulation based on a bit-string model of TCR-pMHC interaction. The model showed that the observed levels of selfrestriction, alloreactivity, and antigen response do not contradict each other, as they can all be achieved within the same—though narrow—subset of the parameter space. We found that affinity-driven selection is compatible with experimental estimates of these latter quantities only if (i) TCRs see more peptide residues than MHC polymorphic residues, (ii) the majority of positively selected clones are deleted by negative selection, (iii) between 1 and 3.6 clonal divisions occur on average in the thymus after completion of TCR rearrangement, and (iv) selection is driven by 103 –105 self peptides. T cell development has been investigated from a variety of perspectives encompassing, on the one hand, molecular events underlying T cell selection and activation, and on the other hand, dynamics of thymic cell populations. However, the two levels of observation are not independent: molecular processes control population dynamics, and both determine T cell repertoire properties. A thorough understanding of T cell selection cannot be achieved without revealing the quantitative relationships between these two levels of description of thymic selection. For example, alloreactivity frequency is, by definition, a global property of the T cell repertoire, not the property of any particular TCR/MHC–peptide ternary complex. Our model was a first step in bridging the gap between microscopic molecular events controlling T cell selection and activation, and their effect on macroscopic properties of the repertoire. This study revealed how quantitative parameters inferred from data on cellular and molecular level experiments contribute to the shaping of the mature T cell repertoire. In the work described above, however, T cell clones were again dealt with one at a time, rather than simultaneously. Hence we could not study the effects of competition, or of a changing MHC-peptide “landscape”, on the development of the TCR repertoire. These remain open issues for future work. We have, however, addressed simultaneous repertoire development in our work on the humoral immune response, described below.

3. Multi-clonal competition driven by mutations: the repertoire shift problem When studying the humoral immune response, with its features of clonal selection of B lymphocytes and affinity maturation of their antigen receptors, it is crucial to take into account the dynamics of a number of simultaneously developing, competing clones. Hypermutation is a genetic-level process, which stands in the basis of the population meta-dynamics through the creation of new clones, so that

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

231

we need at least a two-level model. The intermediate level, that of molecular interactions between the receptor and its nominal antigen, has the important role of determining the fate of B cell clones. For the purpose of many models, however, receptor–antigen interactions can and has been modeled using a single parameter—the affinity of the interaction—because affinity is by far the main factor affecting the fate of each B lymphocyte during the response. This is the strategy that was employed in our study of the problem of repertoire shift [42]. Repertoire shift is a phenomenon in which antibodies dominating a primary response to antigen are severely depleted, or even absent, in the secondary response. The antibodies that dominate the secondary response use variable (V) region genes that are different from those of the primary response. This shift is the rule rather than the exception in responses to T-dependent antigens, where hypermutation and affinity maturation occur. Repertoire shift stands in apparent contradiction to the classical affinity-maturation paradigm, which predicts that dominant clones in the primary response would be most readily recruited into memory responses. Why would a B cell clone, which is dominant in the primary response, be missing from the secondary response? When one considers the role of somatic hypermutation, a way to reconcile repertoire shift with the affinity maturation paradigm becomes apparent: the B cell clone expressing the primary-dominant antibody may fail to dominate the memory response because it is more likely than lower-affinity clones to suffer deleterious mutations. The process of affinity maturation is generally thought to improve the affinity of antibodies to antigen through mutation and selection. If, however, a B cell clone begins with a relatively high affinity to the antigen, a random mutation would be more likely to decrease rather than increase its affinity. At the same time, other, previously minor clones, which can improve through mutation, may become dominant clones. This hypothesis can best be visualized by assuming that, in sequence space, the “landscape” of antibody affinity to the antigen is quite rugged [35], with many local “peaks” and “valleys”. The process of affinity maturation at most takes each antibody from where it started on this landscape only to the nearest local peak, which may not be the point with absolute highest affinity to the antigen in question. The antibody that dominated the primary response may thus end up in a local affinity peak that is lower than the final peaks reached by antibodies that started in different, initially lower, points on the landscape. Indication that the affinity landscape is indeed “rugged” comes from findings of single mutations that confer a large (up to 10-fold) increase in affinity. As only the end results of mutation and clonal selection can be observed experimentally, we chose to use mathematical modeling and computer simulations of the dynamics of B cell clones to test the above hypotheses. Our model consists of an array of B cell sub-clones, which can proliferate, differentiate and mutate in response to antigenic stimulation (Fig. 2). Each sub-clone, i, is assigned a number, i , that represents the affinity of interaction between the cell’s antigen receptor and its nominal antigen. Mutations occur only in the germinal center compartment and are expressed as changes in the value of i . Our “landscape” hypothesis states that the better the receptor is to begin with, the less likely it is to improve through mutation. To express this hypothesis in the model, we originally introduced an artificial distribution of mutation-derived affinity changes [42]. Currently, however, we use a bit string model of the mutation process, with one bit string representing the antigen and another representing the clone’s BCR (Fig. 3). However, when the affinity between the antigen and the BCR was defined as the total number of matches between the two respective strings, as in Fig. 3, the distribution of affinity values resembled a normal distribution with an average matching of L/2 bits. This does not reflect the actual distribution in vivo, in which high-affinity cells are much more rare than low-affinity ones. Hence it also does not give the correct distribution of the results of mutations, because in this scheme each bit flip can

232

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

β

µn

Naive η*ma η*na

µm

Centroblast

memory

η*ap

ρ*g

* ηbc * ηcb

* ηcm

µa

Activated η*ab ρa

mutation rate: v

Centrocyte

* ηcp

µc

plasma

µp

Fig. 2. A schematic representation of the humoral response simulation. The simulation is structured as an array of B lymphocyte clones, each associated with a number representing its relative affinity to the antigen. Each clone contains the following B cell subsets: naive, activated, plasma, centroblast (CB) and centrocyte (CC) germinal center cells, and memory cells. The transitions between the subsets are depicted here as arrows, and marked with an asterix (∗ ) if they are antigen-dependent. Activated and GC cells may also proliferate (denoted by small arrows which return to the same compartment they exit from). Differentiation pathways (solid arrows) and death (dotted arrows) are also shown.

String length = L = 12 Antigen string :

1

0

0

0

1

1

0

1

0

Antibody/BCR string :

1

0

1

0

0

0

1

1

0

0

1

0

Matching :

0

0

1

0

1

1

1

0

0

1

0

1

0

1

1

0

0

1

0

1

1

1

Total affinity = number of matches/L = 6/12 = 0.5 1 0 1 Mutation = change of a bit in the antibody string :

0

1

Fig. 3. A simple version of the Ag/Ab binding bit string model, where affinity is modeled by the number of bit-to-bit matches between the strings, divided by the string’s length. Mutations are modeled by single bit flips.

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

233

only change the affinity by 1/L, while in the real system mutations may change the affinity by varying degrees, from common mutations that make little difference to rare mutations that change the affinity by an order of magnitude or more. In the current model, we emulate this phenomenon by grouping the bits into groups of several bits (the number of bits in a group is a parameter, and can vary between 2 and 20). In each group, we require that the number of matching bits exceed a certain threshold (e.g. 3 bits with a group size of 5), and the affinity is defined as the number of groups in which the number of matching bits exceeds the threshold. This model does give the correct mutation distribution. Members of each clone may be naive (with the number of naive cells represented by Ni ), activated (Ai ), centroblast (CBi ), centrocyte (CCi ), memory (Mi ) or plasma (Pi ) cells. There is constant input of naive cells from the bone marrow, while further differentiation (emigration of cells to the next compartment) depends on interaction with antigen, proportionally to the affinity of the cell’s receptor to the antigen and to the amount of antigen present. Proliferation is also antigen-dependent, and assumed to occur only in the activated and germinal center compartment. Cell death may occur in all compartments. The transition from germinal center to memory cells is assumed to be antigen-independent, and memory cells can also be re-activated by antigen. Antigen depletion results from consumption by activated and plasma cells. The temporal evolution of cells in each clone can be summarized by the following equations, where j denotes the daily proliferation rate of subset j, j representing n (naive), a (activated), b (germinal centroblast), c (germinal centrocyte), m (memory) or p (plasma). Similarly, j k is the daily rate of transition from subset j to subset k, and j is the daily cell death rate of subset j. Thus, naive B cells are being “born” out of the bone marrow with a rate of  cells per clone per day, die with a rate j , or become activated with a rate na . Once activated, B cells proliferate with a rate a , die with a rate a , become plasma cells with a rate ap , or go to the germinal centers and become centroblast with a similar rate ab . Similar equations hold for activated, germinal center centroblast and centrocyte, memory and plasma cells. dNi =  + (n − n − R i na )Ni , dt dAi = R i na Ni + R i ma Mi + [R a − a − R i ( acb + ap )]Ai , dt dCBi = R i acb Ai + RG i cb CCi + (b RG i − bc )CBi , dt dCC i = bc CBi + [(c − cb − cp − cm )RG i − c ]CCi , dt dMi = RG i cm CCi + (m − m − ma R i )Mi , dt dPi = RG i cp CCi + R i ap Ai − m Pi . dt Note that our model is based on the assumption that memory and plasma cells belong to the same lineage. Additionally, the only inter-clone interaction in this model is competition for antigen. No other effects (e.g. idiotypic regulation) are modeled here, for simplicity. R denotes the relative antigen concentration outside of the GCs, (measured relative to its initial concentration, so t hat the maximum concentration

234

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

equals 1), and it decreases with time according to
k (Ai + f P i )i dR = −R i=1 , dt Bmax where f is the factor denoting the relative efficiency of plasma cells, compared to activated cells, in eliminating the antigen (f ?1), and Bmax is a scaling factor. RG , the antigen concentration inside the GCs, which is assumed to decrease much more slowly than that of free antigen due to sequestration on follicular dendritic cells, varies with time according to dRG =  R −  RG dt that is, antigen from outside the GCs flows inside with a rate  day−1 , and is lost with a rate  day−1 . A version of the model in which antigen consumption in the GC is proportional to the number of GC cells is currently being tested. Most transitions (naive to activated, activated to CB or to plasma, memory to activated), and the proliferation of activated or germinal center cells, depend on antigen concentration in a linear manner. Antigen-dependent transitions are multiplied by i , in order to distinguish between clones with different strengths of interaction with the antigen. The transition from germinal center to memory cells is antigen-independent. Memory cells can also be re-activated by antigen. The rates of all these processes were taken from the literature or from previous models. Mutations were modeled by taking cells from the GC compartment and allowing them to “mutate”, that is, randomly obtain new values of their affinity number. The rate of mutation in the simulations was adjusted to fit what is known from experimental systems, taking into account the probabilities that a mutation will be silent, lethal, neutral or advantageous. For a single clone, a mathematical analysis of the differential equations describing the above model may have been sufficient. However, we were interested in the much more complex dynamics of affinity maturation, which emerge when a large number of clones compete for a given antigen. Additionally, we included in our model the genetic-level process of hypermutation, which is discrete and stochastic. Hence we turned to computer simulations of this model. The simulations integrate, in each “time step” and for each sub-clone, the equations representing the number of cells in each subset within this sub-clone. The stochastic component—that of mutation—is included in the model as follows. Within each clone, the number of centroblasts which have divided in each time step is recorded. For each of the daughter cells (that is, for twice the number of divisions), a “mutation attempt” is performed on the bit-string representing the receptor gene of the current sub-clone. The resulting “mutated” string is subjected to “selection”, that is, it is compared to the antigen string to determine the change in affinity to antigen which has resulted from the mutation. The mutated cell can only form a new sub-clone if the change in affinity is positive and passes a selection criterion defined for that run—for example, if it represents a change of over 10% in affinity. Cells that have mutated their receptors but have failed the selection may die as a result—that is, the number of failures is deducted from the total number of centroblasts. Simulations which did not take into account the above-described “landscape” hypothesis, assuming instead that production of memory cells is proportional to the clone’s initial affinity, generated the picture predicted by the affinity maturation paradigm: clones that dominated the primary response continued to dominate the secondary response, that is, there was no repertoire shift. However, the total numbers of cells at the peak of each response, the timing of the peak and the decline of the response, all agreed with experimental observations [18,23], confirming that we have chosen a biologically reasonable set

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

235

of parameter values. Hence we proceeded to incorporate alternative hypotheses concerning the cause of repertoire shift into our model. (Since, on one hand, the experimental observations were only rough estimates of cell numbers, and, on the other hand, the order of magnitude of most parameter values was known, a manual data fitting has sufficed here; the details and parameter values are given in [42]). To incorporate the “landscape” hypothesis in the model, we have set mutation probabilities such that most mutation-induced affinity changes are small, while large “leaps” in affinity are relatively rare, as described above. Selection was applied in our model in the following way. If the post-mutation affinity is smaller than the pre-mutation affinity, the cell fails selection and dies. A cell is positively selected, forming a new sub-clone, only if mutation has improved its affinity by a certain factor, which is a parameter in the simulation. If the improvement is smaller, the mutation is ignored. Simulations with mutation and selection thus implemented give a rich and realistic picture of repertoire shift, evident in all relevant cell subsets—the activated, plasma, GC and memory B cells. If we extend our simulations to include a third and even a fourth antigenic challenge, the clones that dominate the secondary response continue to dominate subsequent responses, in agreement with experimental observations. This means that the optimal receptor for most antigens can be found by the immune system following the first encounter with the antigen, even though it does not exist in the germline repertoire. Repertoire shift is more pronounced if the distribution of mutations is made narrower, that is, if advantageous mutations are rare. This supports our key point, that repertoire shift results from the destructive aspect of hypermutation on antibodies that had a high affinity to the antigen prior to mutation. Second, as expected, a similar effect is observed when increasing the relative advantage conferred by mutations, or the stringency of the selection operating on newly generated mutants. A related question is whether the average affinity of antibodies to the antigen increases between the primary and the secondary response. Experimental studies [16] concluded that the main effect of somatic hypermutation on the primary antibodies must be destructive, making way for the appearance of new clones, which supports our present conclusion. Studies of the humoral immune response have traditionally focused on the positive aspect of somatic hypermutation. The paradigm of affinity maturation, stating that B cells taking part in the secondary response are “improved versions” of those that dominated the primary response, is often confused with somatic hypermutation. One should keep in mind that somatic hypermutation is only the mechanism generating candidates for affinity improvement, on which selection then operates. One must remember that affinity-improving mutation is more likely to be a rare event in a multitude of deleterious mutations. A thorough understanding of the humoral immune response must acknowledge the destructive aspect of hypermutation. From the point of view of modeling, it is noteworthy that the rich dynamics of hypermutation, affinity maturation and repertoire shift were obtained in the model using only ten initial B lymphocyte clones. This is comparable to what was found in in vivo studies [18,45]. The key points were the realistic description of B cell subsets within each clone, combined with a model of mutation and selection, which captured the essential features of the “landscape” hypothesis.

4. Lineage tree analysis of mutated immunoglobulin genes The generation of “lineage trees” or “dendrograms” to visualize the lineage relationships of B cell mutants in the GCs has been used in the past to confirm the role of the GC as the location of somatic hypermutation [14,22,26], to identify lineage relationships between cells from independent GCs [46] or

236

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

0 Split Nodes 1 2

7

4

5

8

9 10

Root This tree as a list: 0 1 Trunk 2 3 Pass-Through Nodes 1 2 4 5 3 3 6 6 4 7 8 9 5 10 11 6 11 11 12 12 12 13 13

Leaves

Fig. 4. An example of a lineage tree.

different tissues [10,11] and from additional processes of diversification such as gene conversion in the rabbit [38–41]. The experimentally generated lineage trees reflect multiple rounds of mutations for each germline V gene that participate in the primary response. We believe that much information about the dynamics of antigen-driven clonal selection during the immune response is contained in the shape of lineage trees deduced from the final responding clones [42]. For example, trees generated from clones during the peak of the primary response are much more “bushy” [13], but trees become less “bushy” as the response progresses [15]. The “pruned” shape of these trees has been referred to as evidence for the destructive character of somatic hypermutation. Other examples of lineage trees drawn to illustrate various aspects of the germinal center reaction, or differences in this reaction under varying circumstances, abound in the literature. So far, however, lineage tree classification has been based only on a qualitative, intuitive assessment of the most obvious shape characteristics. Hence we set out to explore whether the information embedded in the mathematical shape characteristics of lineage trees can in any way be quantifiable, and whether it can be shown to correlate with the dynamics of the underlying immune response. The objective of our first study was to develop a rigorous computer-aided algorithm for extracting the information contained in lineage trees, using the tools of mathematical graph theory. The algorithm we developed characterizes trees according to their various graph-theoretical measures, as described below, and we have found correlations between these measures and the dynamical parameters of the GC response that generated the trees. A lineage tree is defined, graphically, as a rooted tree where the nodes correspond to B cell receptor gene sequences (Fig. 4). For two nodes X and Y, we say that Y is a son of X if the sequence corresponding to Y is a mutant of the sequence corresponding to X, which differs from X by only one mutation, and is one mutation further than X away from the original (germline) gene, that is, the root. Two B cells with identical receptor genes will thus correspond to the same node. A lineage tree depicts the maturation process of a B cell at a certain moment of observation—it consists only of the Ig gene sequences of cells that were sampled at that moment and their ancestors back to the root. Most of these cells’sequences are represented by leaves in the tree; others are internal nodes (see below). The ancestors are not necessarily sampled at the time of observation. Nodes in the tree can be either the root node (always the node numbered zero in our definition), leaves (sequences of cells that had no sons at the time of observation—see nodes 7, 8, 9, 10 and 13 in Fig. 4), or internal nodes. Internal nodes can be either split nodes—those with more than one son (nodes 1, 2 and 4); or pass-through nodes—those with exactly one son (nodes 3, 5, 6, 11 and 12). The trunk is the distance between the root to the first split node, which equals 1 in this case.

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

237

Table 1 Parameters measured on each tree Parameter

Abbreviation

Total number of leaves Total number of nodes, including the root Number of internal nodes Number of pass-through nodes Trunk length Minimum, maximum and average (per tree) path length Minimum, maximum and average distance from a leaf to the last (closest to leaf) split node Minimum, maximum and average distance from a leaf to the first (closest to root) split node Minimum, maximum and average distance from the root to any split node Minimum, maximum and average distance between adjacent split nodes Distance between the root to the maximal (in terms of number of sons) split node Minimum, maximum outgoing degree Outgoing degree averaged over all nodes Outgoing degree averaged over all SPLIT nodes The root’s outgoing degree

L N IN PTN T MinPL, MaxPL, AvgPL MinDLSN, MaxDLSN, AvgDLSN MinDLFSN, MaxDLFSN, AvgDLFSN MinDRSN, MaxDRSN, AvgDRSN MinDASN, MaxDASN, AvgDASN DRMSN MinOD, MaxOD AvgOD AvgOD2 RootD

The shape of lineage trees is quantified by measuring a number of properties of the graph describing the tree. For quite a few of these properties, the maximum, minimum and average values per tree may be measured. The complete list of parameters measured is given in Table 1. The meaning and usefulness of each of these parameters is discussed below. Our tree-measurement computer program reads a tree in the format described above and measures the graphical parameters, creating a text report. We intend to create, in the near future, a user-friendly interface for this program and make it available on the web for any researcher interested in analysis of their lineage trees; for the time being we are willing to analyze data sent to us upon request (see http://repertoire.ls.biu.ac.il/TREES for details). Note that, for the purpose of our analysis, we are not interested in the properties of the individual cells or clones represented by the lineage tree, but rather in the overall characteristics of the lineage tree as a graphical entity. We demonstrated in the first papers [2,9] that the information extracted using our algorithm is indeed valuable in revealing the dynamics of hypermutation and antigen-driven selection in germinal centers. We then proceeded to examine the influence of experimental parameters on the graphical properties of lineage trees. Lineage tree parameters depend on the following four factors: (1) the number of different sequences within each clone obtained in a given experiment—the “pick size”. This in turn is affected by the method of cell collection for DNA analysis, which also influences data resolution, as discussed below; (2) the intrinsic potential for improvement via hypermutation possessed by each antigen receptor gene—the closer the original receptor is to the optimal binder of the particular antigen, the fewer the mutations that could lead to improvement of binding, and the higher the probability that a mutation will be deleterious or reduce the receptor’s affinity to the antigen [42]; (3) the rate and duration of the hypermutation process itself; and (4) the stringency of selection operating on the mutating cells. The latter two factors depend on the tissue in which the response occurs, the antigen, and all other factors regulating response dynamics. It is important to understand the dependencies of the measured parameters on these four factors, as our

238

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

goal is to minimize the influence of arbitrary factors such as pick size and individual gene potentials, and learn as much as possible about the intrinsic processes of hypermutation and selection. Thus, the goal of our next few studies [7,8] was to deepen our insight about the meaning of each parameter measured on Ig gene phylogenetic trees. We have done so by analyzing all trees available in the literature from a defined system—the murine primed immune response—and comparing them to trees derived from human tissues. Many interesting relationships between the measured parameters were revealed, leading to some biological insights. In particular, we addressed the artifacts introduced by experimental methods of cell sampling, the meaning of each parameter measured on the tree graphs, and the differences between the dynamics of the humoral immune response in different lymphoid tissues. Our findings may be summarized as follows. First, the cell sampling method used for DNA analysis is crucial in determining the resolution of the data. Whole tissue samples give data whose resolution is much lower than single cell or even whole GC analysis, and should probably be avoided. The best way is probably single-cell or small sample analysis, provided a sufficient number of cells are picked and analyzed. Even when GCs are not available (e.g. when only patient blood samples are available) single cell or small sample analysis is probably the only way to perform a meaningful analysis. Second, for subtle differences in mutation rate, duration, or selection to be discovered, analysis of larger data sets with high resolution is required. Third, the choice of parameters for analysis depends on the question being asked. For example, if we are looking at rate or duration of the diversification process (as in comparing patients in different disease stages), the pick size is irrelevant, as we will be looking at path lengths along the lineage tree; however, one must assume that the rate of mutation is the same in all cases in order to compare duration, or vice versa, as both will affect path lengths in the same way. Lineage tree study may, at this point, be developed in various directions. In one study [27], we have examined the generation of immunoglobulin diversity in species, such as rabbits and chickens, which use gene conversion much more extensively than rearrangement during the primary diversification of their immunoglobulin repertoire, and combine it with hypermutation during primary secondary diversification. This method may have a great potential for studying the dynamics of immune responses in various disease situations—chronic infectious diseases, autoimmune diseases and lymphomas of GC B cell origin, to name a few—by comparing the trees generated in these situations to those generated in normal immune response dynamics. However, one should be cautious in making such comparisons. The cell sampling method used for DNA analysis, the tissue from which B cells were taken, the species of experimental animal, the intensity of the diversification process as reflected by the rate of hypermutation, and the stringency of selection exerted (if at all) on the mutating cells, may all affect the measured tree parameters in various ways, as shown above. In addition, there are theoretical questions that may be addressed via simulations. One question is the effect of selection parameters on tree properties. Various theoretical models were suggested for how selection operates on GC B cells [3,19–21,29–31,33,42–44]. We may address these models, and the effect of model parameters on lineage tree characteristics, via simulation studies. A different question is whether there are additional parameters that may be defined on trees and yield additional insight into response dynamics. We intend to carry on the research along these lines.

5. Conclusions In this paper, I have reviewed several examples of theoretical models and analysis tools of lymphocyte repertoire generation and selection. My aim was to point at the difficulties posed by this area of theoretical

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

239

work, and demonstrate various types of solutions. The main message I would like to convey is that computational tools, in particular computer simulations and bioinformatical analysis, are necessary and must complement classical mathematical models for grasping, visualizing and analyzing lymphocyte repertoire meta-population dynamics. This is true provided the simulation and analysis tools are geared to deal with the immense complexity of the immune system. One must, of course, bear in mind that no model can—or should—attempt to be a complete description of reality. However, over-simplification must be avoided. Appropriately formulated models and methods of analysis are very helpful whenever the need arises for complex, quantitative “thought experiments” to accompany experimental investigations.

Acknowledgements The author gratefully acknowledges collaborations with V. Detours, D.K. Dunn-Walters, S. Litwin, R. Mage, A. Perelson, H. Piper and M. Shannon, co-authors of some of the studies reviewed above. The studies reviewed here were supported in parts by an Israel Science Foundation Grant number 759/01-1, and an Yigal Alon Fellowship (to R.M.).

References [1] H. Arakawa, T. Shimizu, S. Takeda, Re-evaluation of the probabilities for productive rearrangements on the  and  loci, Internat. Immunol. 8 (1996) 91–99. [2] M. Banerjee, R. Mehr, A. Belelovsky, J. Spencer, D. Dunn-Walters, Age and tissue-specific differences in human germinal centre B cell selection, Eur. J. Immunol. 32 (2002) 1947–1957. [3] F. Celada, P.E. Seiden, Affinity maturation and hypermutation in a simulation of the humoral immune response, Eur. J. Immunol. 26 (1996) 1350–1358. [4] C. Coleclough, R.P. Perry, K. Karjalainen, M. Weigert, Aberrant rearrangements contribute significantly to the allelic exclusion of immunoglobulin gene expression, Nature (London) 290 (1981) 372–378. [5] V. Detours, R. Mehr, A. Perelson, A quantitative theory of affinity-driven T cell repertoire selection, J. Theoret. Biol. 200 (1999) 389–403. [6] V. Detours, R. Mehr, A. Perelson, Deriving quantitative constraints under which T-cell selection operates from data on the mature T-cell repertoire, J. Immunol. 164 (2000) 121–128. [7] D. Dunn-Walters, A. Belelovsky, H. Edelman, M. Banerjee, R. Mehr, The dynamics of germinal centre selection as measured by graph-theoretical analysis of mutational lineage trees, Clin. Dev. Immunol. 9 (2003) 233–245. [8] D. Dunn-Walters, H. Edelman, R. Mehr, Immune system learning and memory quantified by graphical analysis of Blymphocyte phylogenetic trees, BioSystems 76 (2004) 141–155. [9] D.K. Dunn-Walters, M. Banerjee, R. Mehr, Age effects on antibody affinity maturation, Biochem. Soc. Trans. 31 (2003) 447–448. [10] D.K. Dunn-Walters, L. Boursier, P.J. Ciclitira, J. Spencer, Immunoglobulin genes from human duodenal and colonic plasma cells are mutated, Biochem. Soc. Trans. 25 (1997) 324S. [11] D.K. Dunn Walters, P.G. Isaacson, J. Spencer, Sequence analysis of human IgVH genes indicates that ileal lamina propria plasma cells are derived from Peyer’s patches, Eur. J. Immunol. 27 (1997) 463–467. [12] S.J. Foster, H.P. Brezinschek, R.I. Brezinschek, P.E. Lipsky, Molecular mechanisms and selective influences that shape the  gene repertoire of IgM+ B cells, J. Clin. Invest. 99 (1997) 1614–1627. [13] J. Jacob, G. Kelsoe, In situ studies of the primary immune response to (4-hydroxy-3-nitrophenyl) acetyl. II. A common clonal origin for periarteriolar lymphoid sheath-associated foci and germinal centers, J. Exp. Med. 176 (1992) 679–687. [14] J. Jacob, G. Kelsoe, K. Rajewsky, U. Weiss, Intraclonal generation of antibody mutants in germinal centres, Nature 354 (1991) 389–392.

240

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

[15] J. Jacob, J. Przylepa, C. Miller, G. Kelsoe, In situ studies of the primary immune response to (4-hydroxy-3nitrophenyl)acetyl. III. The kinetics of V region mutation and selection in germinal center B cells, J. Exp. Med. 178 (1993) 1293–1307. [16] U. Kalinke, E.M. Bucher, B. Ernst, A. Oxenius, H.-P. Roost, S. Geley, R. Kofler, R. Zinkernagel, H. Hengartner, The role of somatic mutation in the generation of the protective humoral immune response against vesicular stomatitis virus, Immunity 5 (1996) 639–652. [17] G. Kalmanovich, R. Mehr, Models for anitgen receptor gene rearrangement. III. Heavy and light chain allelic exclusion, J. Immunol 170 (2003) 182–193. [18] G. Kelsoe, Life and death in germinal centers (redux), Immunity 4 (1996) 107–111. [19] T.B. Kepler, A.S. Perelson, Cyclic re-entry of germinal center B cells and the efficiency of affinity maturation, Immunol. Today 14 (1993) 412–415. [20] C. Kesmir, R.J. De Boer, A mathematical model on germinal center kinetics and termination, J. Immunol. 163 (1999) 2463–2469. [21] S.H. Kleinstein, J.P. Singh, Toward quantitative simulation of germinal center dynamics: biological and modeling insights from experimental validation, J. Theoret. Biol. 211 (2001) 253–275. [22] C. Kocks, K. Rajewsky, Stepwise intraclonal maturation of antibody affinity through somatic hypermutation, Proc. Natl. Acad. Sci. 85 (1988) 8206–8210. [23] M.H. Kosco-Vilbois, H. Zentgraf, J. Gerdes, J. Bonnefoy, To ‘B’ or not to ‘B’ a germinal center?, Immunol. Today 18 (1997) 225–230. [24] M. Malissen, J. Trucy, E. Jouvin-Marche, P.-A. Cazenave, R. Scollay, B. Malissen, Regulation of TCR
and  gene allelic exclusion during T-cell development, Immunol. Today 13 (1992) 315–322. [25] C.A. Mallick, E.C. Dudley, J.L. Viney, M.J. Owen, A.C. Hayday, Rearrangement and diversity of T cell receptor  chain genes in thymocytes: a critical role for the  chain in development, Cell 73 (1993) 513–519. [26] T. Manser, Evolution of antibody structure during the immune response, J. Exp. Med. 170 (1989) 1211–1230. [27] R. Mehr, H. Edelman, D. Sehgal, R. Mage, Analysis of mutational lineage trees from sites of primary and secondary immunoglobulin gene diversification in rabbits and chickens, J. Immunol. 172 (2004) 4790–4796. [28] R. Mehr, M. Shannon, S. Litwin, Models for antigen receptor gene rearrangement. Biased receptor editing in B cells—implications for allelic exclusion, J. Immunol. 163 (1999) 1793–1798. [29] M. Meyer-Hermann, A mathematical model for the germinal center morphology and affinity maturation, J. Theoret. Biol. 216 (2002) 273–300. [30] M. Meyer-Hermann, A possible role of chemotaxis in germinal center formation. A possible role of chemotaxis in germinal center formation, Internat. Immunol. 14 (2002) 1369–1381. [31] M. Meyer-Hermann, A. Deutsch, M. Or-Guil, Recycling probability and dynamical properties of germinal center reactions, J. Theoret. Biol. 210 (2001) 265–285. [32] D. Nemazee, Theoretical limits to massive receptor editing in immature B cells, in: G. Kelsoe, M.F. Flajnik (Eds.), Current Topics in Microbiology and Immunology, vol. 229, Somatic Diversification of Immune Responses, Springer, Berlin, Germany, 1998, p. 163. [33] M. Oprea, A.S. Perelson, Somatic mutation leads to efficient affinity maturation when centrocytes recycle back to centroblasts, J. Immunol. 158 (1997) 5155–5162. [34] E. Padovan, G. Casorati, P. Dellabona, S. Meyer, M. Brockhaus, A. Lanzavecchia, Expression of two T cell receptor
chains: dual receptor T cells, Science 262 (1993) 422–424. [35] A.S. Perelson, S.A. Kauffman (Eds.), Molecular evolution on rugged landscapes: proteins, RNA and the immune system. Santa Fe Institute Studies in the Sciences of Complexity, vol. 9, Addison–Wesley, California, 1991, pp. 73–176. [36] H. Piper, S. Litwin, R. Mehr, Models for antigen receptor gene rearrangement T Cell receptor editing: allelic exclusion or inclusion?, J. Immunol. 163 (1999) 1799–1808. [37] M.Z. Radic, M. Zouali, Receptor editing, immune diversification and self-tolerance, Immunity 5 (1996) 505–511. [38] E. Schiaffella, D. Seghal, A.O. Anderson, R.G. Mage, Gene conversion and hypermutation during diversification of VH sequences in developing splenic germinal centers of immunized rabbits, J. Immunol. 162 (1999) 3984–3995. [39] D. Seghal, H. Obiakor, R.G. Mage, Distinct clonal diversification patterns in young appendix compared to antigen-specific splenic clones, J. Immunol. 168 (2002) 5424–5433. [40] D. Seghal, E. Schiaffella,A.O.Anderson, R.G. Mage,Analysis of single B cells by PCR reveals rearrangedVH with germline sequences in spleen of immunized rabbits: implications for B cell repertoire maintenance and renewal, J. Immunol. 161 (1998) 5347–5356.

R. Mehr / Journal of Computational and Applied Mathematics 184 (2005) 223 – 241

241

[41] D. Seghal, E. Schiaffella, A.O. Anderson, R.G. Mage, Generation of heterogeneous rabbit anti-DNP antibodies by gene conversion and hypermutation of rearranged VL and VH genes during clonal expansion of B cells in splenic germinal centers, Eur. J. Immunol. 30 (2000) 3634–3644. [42] M. Shannon, R. Mehr, Reconciling repertoire shift with affinity maturation: the role of deleterious mutations, J. Immunol. 162 (1999) 3950–3956. [43] M. Shlomchik, P. Watts, M. Weigert, S. Litwin, Clone: a Monte-Carlo computer simulation of B cell clonal expansion, somatic mutation, and antigen-driven selection, Curr. Top. Microbiol. Immunol. 229 (1998) 173–197. [44] B. Sulzer, L. van Hemmen, A.U. Neumann, U. Behn, Memory in idiotypic networks due to competition between proliferation and differentiation, Bull. Math. Biol. 55 (1993) 1133–1182. [45] Y. Takahashi, R.P. Dutta, D.M. Cerasoli, G. Kelsoe, In situ studies of the primary immune response to (4-hydroxy-3nitrophenyl)acetyl. V. Affinity maturation develops in two stages of clonal selection, J. Exp. Med. 187 (1998) 885–895. [46] K.A. Vora, K. Tumas-Brundage, T. Manser, Contrasting the in situ behaviour of a memory B cell clone during primary and secondary immune responses, J. Immunol. 163 (1999) 4315–4327. [47] D.L. Wood, C. Coleclough, Different joining region J elements of the murine  immunoglobulin light chain locus are used at markedly different frequencies, Proc. Natl. Acad. Sci. 81 (1984) 4756–4760. [48] B. Yaish, R. Mehr, Models for the dynamics and order of immunoglobulin isotype switch, Bull. Math. Biol. 67 (2005) 15–32.