Stochastic approximation to the specific response of the immune system
Abstract
We develop a stochastic model to study the specific response of the immune system. The model is based on the dynamical interaction between Regulatory and Effector CD4+ T cells in the presence of Antigen Presenting Cells inside a lymphatic node. At a mean field level the model predicts the existence of different regimes where active, tolerant, or cyclic immune responses are possible. To study the model beyond mean field and to understand the specific responses of the immune system we use the Linear Noise Approximation and show that fluctuations due to finite size effects may strongly alter the mean field scenario. Moreover, it was found the existence of a certain characteristic frequency for the fluctuations. All the analytical predictions were compared with simulations using the Gillespie’s algorithm.
keywords:
Immunology, Linear Noise Approximation, Stochastic modelling, T cellsPacs:
02.70.Lq, 75.10.Hk, 75.40.Mg, 75.70.Kw1 Introduction
The purpose of the immune system is to detect and neutralize the molecules, or cells, dangerous for the body without damaging the healthy cells. A fundamental process of the immune system is the maintenance of selftolerance, i.e., the prevention of harmful immune responses against body components [1]. The biological significance of this process becomes very patent upon its failure during pathological conditions known as autoimmune diseases.
The risk of autoimmunity cannot be dissociated from the capacity of the immune system to cope with diverse and fast evolving pathogens. The latter is achieved by setting up a vast and diverse repertoire of antigen receptors expressed by lymphocytes, which as a whole is capable of recognizing any possible antigen. Most lymphocytes have a unique antigen receptor (immunoglobulin in Bcells and TCR in T cells) that is encoded by a gene that results from somatic mutation and random assortment of gene segments in lymphocyte precursors. The randomness in the generation of antigen receptors makes it unavoidable that lymphocytes with receptors recognizing body antigens are also made. These autoreactive lymphocytes can potentially cause autoimmune diseases if their activation and clonal expansion are not prevented. The question is how is this avoided in healthy individuals?
Regulatory CD4+ T cells, which express forkhead box protein 3 (FoxP3) are enriched in the CD25 pool of healthy individuals have been gaining increasing relevance in immunology [2]. Many lines of evidence indicate that these cells play a key role in the development of natural tolerance and in the prevention of autoimmune pathologies, by controlling the activation and proliferation of other autoreactive lymphocytes. The functional significance of these cells has broadened, as they were shown to modulate the immune response against pathogens, preventing the associated immunopathology , and the rejection of transplants.
Some years ago, León et al [3] proposed a simple model, consistent with these results able to explain the mechanism of immunological selftolerance. The main hypothesis behind this model is that regulatory cells inhibit the proliferation of effector T cells, but depend on the latter for their own proliferation. The interaction between these cells is mediated by the presence of Antigen Presenting Cells (APC).
The model, was developed to study the response of macroscopic populations of lymphocytes, i.e.,the lymphocytes that may be present in a single lymph node and was technically approached studying the stability or numerically solving the corresponding differential equations [3]. This approach however, might be inappropriate to study the specific clonal responses in the immune system.
The lymphocytes that exist in a lymph node, are really divided in approximately or different clones[4]. I.e. sets of lymphocytes which recognize the same antigen and participate in their own specific, rather independent, immunological response. Therefore, the number of cells involved in any particular clonal response is far smaller: or T cells, which interact specifically with a small subset of all the antigen presenting cells available in the lymph node. One must naturally expect to found clones of different sizes, perhaps reflecting the ubiquitous nature of the antigen recognized by the T cells (i.e. the amount of APCs which are recognized). However, the small size expected for most T cells clones would made differential equations an inappropriate tool to describe their dynamics.
In this work, we first reformulated the model of Kalet and co. [3] within a more general stochastic frameworks of interactions between regulatory (R) and effector (E) Tcells. From this formulation we derive, on one hand a mean field description similar but not identical to the ones already known in the literature for the response of the entire immune system [5]. On the other, a selfconsistent approach to unveil the role of the fluctuations in the specific response of the system.
The rest of the paper is organized as follows. In the next section, that may be bypassed in a first reading of the manuscript, we present the main mathematical techniques involved in our work. Section 3 introduces and motivates the model. In section 4 we first compare, at the meanfield level, the predictions of our model with those of similar approaches in the literature. Then, we discuss how fluctuations may affect these predictions highlighting their relevance for the clonal response of the immune system. Finally, the conclusions and some remarks suggesting possible extensions of our work appear in Section 5.
2 Mathematical Techniques
When the effects of fluctuations are relevant, as it is often the case in biological systems, a modeling approach based on ordinary differential equations is not appropriate.
This is well understood when the source of noise are the thermal fluctuations. But, while less mentioned, it is also important when the noise is induced by the finite number of elements in the populations, i.e., molecules, cells, or individuals [6].
A more convenient framework to describe these problems is a stochastic approach taking into account the probabilistic nature of every process. The natural procedure is to write down a differentialdifference equation that keeps track in time of the probability that the system would have a certain configuration (). Under the Markov assumption it is called the Master equation [7]
(1) 
where is a vector whose components are the numbers of elements of every different species in the system, is the probability of having the state at time , and is the probability of the transition per unit time.
Unfortunately, the Master equation can only be solved in special cases, usually, when the transition probabilities are linear or constant in the state variables [8]. In more general situations we must, instead, resort to approximation methods. These approximation methods may be divided into two broad classes: numerical simulation algorithms, and perturbative calculations. Under the first class the Gillespie [9] algorithm has become the standard tool in the community (see appendix for a short description of the algorithm). It is a Monte Carlo’s algorithm that simulates a single trajectory consistent with the unknown probability distribution characterized by the Master equation. It guarantees to find the correct solution of the Master equation, but it is often time consuming and it is not suitable for parameter exploration. On the contrary, perturbative methods lack the confidence given by a full and formal solution of the Master equation, but may provide important insight on the behavior of the system and on the relevance of the parameters.
In this work we present two different perturbative solutions. The Linear Noise Approximation ( LNA) [10] which is an expansion around the large size solution of the corresponding Master equation. This gives us a tool to understand the relative size of the fluctuations and their relations with the parameters of the model. On the other hand, the Effective Stability Approximation ( ESA) [11], which is an extension of the LNA, provides a quantitative characterization of the effect of the fluctuations in the stability of the phases predicted by the corresponding mean field model.
2.1 The Linear Noise Approximation
An important difficulty in solving the Master equation arises from the discreteness of the variables involved (). As long as the number of elements increases the system evolution turns more regular and the meanfield equations give a more accurate description. The Linear Noise Approximation is a systematic approximation method that rests upon the assumption that the deterministic evolution of the concentrations in the system can be meaningfully separated from the fluctuations, and that these fluctuations scale roughly as the square root of the system size .
Moreover, in systems where the size of the populations differs in orders of magnitudes, it is expected that the respective concentrations fluctuate within different scales. For this reason, we consider as many parameters as species are involved in the reactions, and write the numbers of individuals per species proportional to these values. Each will be a characteristic size scale for every particular specie.
Therefore, under the LNA, the population numbers per specie are written as:
(2) 
where is the deterministic prediction for the concentration of the species with respect to the parameter , and measures the fluctuations around . Note that formalizing the LNA in this way, there isn’t any a priory assumption about the system size and concentrations higher than one are allowed. If for every specie the standard formulation is recovered.
Now we shall use the continuous variables instead of the integers to write the probability distribution . Let us group the magnitudes into the macroscopic concentrations vector , and the into the fluctuations vector . The matrix is a diagonal matrix whose principal elements are the extensive parameters with the units of the volume to which concentrations are refered.
Let us define the following variables:
(3) 
thus, we have
(4) 
Taking (4) into account to expand the Master equation (1) in powers of , at the leading order the following system of differential equations appears:
(5) 
This corresponds to the deterministic rate equations that are often used to describe the dynamics of such systems at a meanfield level.
Here the time scale has been modified following and is the stoichiometry coefficient, the number in what the species varies when the process occurs. These coefficients will correspond with the element of row and column of the stoichiometry matrix S.
In the next order approximation one finds:
(6) 
where
(7) 
and
(8) 
and are respectively the elements () of the matrices A and D, where, A can be identified with the Jacobian matrix of the system of differential equations (5).
Note that equation (6) is a FokkerPlank equation that characterizes the probability distribution for the fluctuations , centered on the macroscopic trajectory .
The matrices A and D are independent of , which appears only linearly in the drift term. As a consequence, the distribution will be Gaussian for all time [12]. In particular, at equilibrium () the fluctuations are distributed with density
and variance determined by
(9) 
where A and D are evaluated on the studied fixed point [12].
The steadystate time correlation function is
(10) 
Moreover, if fluctuations are important it may be useful to study their properties. With this aim, a very powerful tool is the Fourier analysis of the equations that govern it. In our case it is the FokkerPlank equation (6), who gives us all the information about the temporal behavior of fluctuations. Given the inconvenient form of this equation for our purpose, it is more reasonable to write it in a completely equivalent formulation more benevolent to investigation using Fourier transforms [13]. The problem can then be formulated as the set of stochastic differential equations of the Langevin type [12]:
(11) 
where is a Gaussian noise with zero mean and correlation function given by
(12) 
The power spectrum of the fluctuations can be found for every specie of the system by averaging the square modulus of the Fourier transformation of . In this way,
(13)  
Being and [13]. With this expression we have an analytic measure of the contribution of every frequency to the oscillatory behavior of the concentrations around each fixed point. Through it, the properties of the fluctuations can be studied.
2.2 Effective Stability Approximation
It is required just a bit of intuition to realize that the noise can affect the predictions of the deterministic rate equations models [11]. This is particularly clear, once we question about the stability of the fixed point predicted by a meanfield model.
In the deterministic analysis, the stability of the model to small perturbations is found by linearizing about the equilibrium point: and the eigenvalues of the Jacobian provide the decay rate of the exponential eigenmodes [14].
If we now consider the fluctuations around the small perturbation , the concentrations remain . The Jacobian matrix , will be expressed in terms of the fluctuations around the steadystate. In the limit , we can linearize J with respect to and the new Jacobian becomes
(14) 
Therefore, new effective eigenvalues , including the influence of the intrinsic noise, characterize the stability of the phases [11].
(15) 
where is proportional to , i.e., inversely proportional to the system size.
When the system size is small, fluctuations are relevant, and in many cases the corrections can not be ignored. This small correction may, indeed, change the sign of the eigenvalues (now ) and as a consequence, to change the stability of the phase.
3 The crossregulatory model
In this section we motivate and present a stochastic model to describe the dynamics of two (possible small) populations of T cells E and R of the immune system.
Some clues on how the regulatory CD4+ T cells suppress the response of other cells have been derived from wellcorrelated in vitro and in vivo experiments. These studies suggest that Tcellmediated suppression is not mediated by soluble factors and require celltocell contact [15][16]. Moreover, regulatory CD4+ T cells can only suppress the response of other cells if the ligands of both cells are expressed by the same Antigen Presenting Cell (APC). It also has been seen that regulatory T cells do not proliferate in vitro when cultivated alone in the presence of APCs.
Despite these strong requirements, it is not yet clear what the nature of this mechanism is. The model requires a set of postulates that summarize the life cycle of these cells, and their interactions. Inspired by the success of León et al.[3] we are going to assume that the proliferation of T cells occurs via its conjugation with an APC. An E (effector) cell could duplicate only if conjugated with an APC which has no R (regulatory) cell. On the other hand, an R cell can duplicate only if conjugated with an APC where at least one E cell is also conjugated. In that way, regulatory T cells exist only in the presence of effector T cells, but the growth of the E cells is regulated by the presence of the R cells. A picture illustrative of these postulates is shown in Figure 1.
In vivo the total number of APCs and their capacities to stimulate T cells may change, and it is perhaps a function of the T cells themselves, but we consider them as a constant population in which every APC has a finite and fixed number of conjugation sites. Then, each site can be empty or occupied by a regulatory or effector T cell and all the cells that belong to the same phenotype (E or R) will have the same probability of being conjugated with any free site of any APC.
Within these postulates we have to deal with the populations of free and conjugated cells for each specie, and the transition between the different states. This makes the model very cumbersome at the stochastic level, and almost intractable. Fortunately the processes of formation and dissociation of the conjugates are relatively fast compared to the overall dynamics of the T cell populations. While a T cell remains conjugated with its APC for a few hours, the T cell mitotic cycle lasts on average 12 hours and the T cell lifespan is at least several days. Then, it is reasonable to assume that only free lymphocytes die, and that the numbers of conjugated T cells of both phenotypes are in quasisteady state.
Therefore, if is the total number of APCs, and , the number of conjugation sites per APC we may define:
Probability per unit time of combination of an E cell and an APC.
Probability per unit time of dissociation of E to its conjugate.
Here and are the numbers of free and conjugated E cells respectively, and and are positive parameters that characterize the processes of occurrence and disappearance of the conjugates of this cellular species.
In the same way, for R cells these probabilities per unit time are given by and respectively.
If now we set and , the stationary state for the formation and dissociation of conjugates for every species reads:
(16) 
After this approximation, our analysis involves only four significant processes, the birth and death of effector and regulatory cells, and is far more tractable. The four processes are represented in Figure 2.
Each event is characterized by an occurrence probability per unit time from the state to the posterior one . It is,
and the stoichiometry matrix corresponding to this system, whose element (i,j) contains the changes produced on the specie when it occurs the process is
(17) 
Now, following the postulates of our model, a new E cell appears in the system if there already exists a combined E cell with an APC that has no R cells simultaneously conjugated. Then, the probability of the increase of in is proportional to the number of APCs (), times the probability to find in it a conjugated E cell (), times the probability that this APC has no R cells simultaneously combined. The later has the form of an Hypergeometric distribution [3], but to gain in simplicity of the expressions we can manage it approximately as , valid when [17]. This takes us to the expression
(18) 
The analysis for the growth of the R population in unit is similar. A new regulatory cell appears in the system with certain probability if there exists already both R and E cells conjugated in the same APC during the same interval of time. The probability of the existence of an E cell conjugated in a given APC where in addition exists space for an R cell can be written as . Multiplying this by the probability of randomly find that specific APC with an R cell combined in it gives the probability for this birth:
(19) 
In equations (18) and (19), and are parameters for the effector and regulatory species that characterize the processes of formation and dissociation of the different conjugates.
Finally, we assume that the death probabilities for each specie are proportional to the numbers of individuals, with proportionality factors for effector and for regulatory.
Let us set the concentration of E cells with respect to the quantity , and will be the concentration of R cells relative to , where and are the characteristic mean values. The probability vector per unit time , that contains the information of the four processes of birth and death is formed by the elements:
4 Results and Discussion
4.1 Meanfield approximation
Once we know all the elements of the transition probabilities per unit time, and the stoichiometry of the events, we are in conditions to write the deterministic rate equations that correspond to our model (see eq. (5) ).
(24)  
(25) 
where the time has been measured according to .
In the steadystate, , the solutions for the deterministic system are given by a fixed pair () of numbers. Depending on the parameters, the stationary solutions with biological sense may be of three types: (), () and ().
Attending to the stability of these solutions we define five phases that are sketched in Figure 3.
Phase 0: In this phase only the trivial fixed point () exists. It is characterized by the relation and within this phase, this fixed point is always stable.
Phase I: In this phase, only the point () is stable.
Phase II: Only the point () is stable.
Phase III: It is a bistable phase. Both points () and () are stable
Phase IV: In this phase, the fixed points presented above are unstable. A limit cycle develops.
This phase diagram is very similar to the one reported in [3]. However, for completeness and because the Phase IV was not predicted before we discuss in detail below the biological interpretation of the five regions.
Since E cells themselves do not mediate the effective immune response, but they trigger it, we interpret the state dominated by regulatory T cells () as tolerance or unresponsiveness, and the state dominated by effector () as an effective immune response or autoimmunity. Such interpretation arises from the fact that when R cells exist in the system, this population competes with the effectors for the APCs, limiting the growth of the E population.
Particularly relevant from the immunological point of view is Phase III. Operating in this region, the model can switch between the states dominated by regulatory or effector T cells in a reversible way [3].
For the fixed points of the kind (), , the eigenvalues of the Jacobian matrix are:
(26) 
and
(27)  
the first one of these eigenvalues is always negative and its corresponding eigenvector is . It guarantees that, if not perturbed along the direction, the fixed point is stable. This means that in those phases where both species can coexist (, and ), once the R population is extinguished, to take the system away from the () state it is necessary that at least one R cell comes from outside the system. The absence of such a perturbation can be interpreted as the lost of tolerance, and the appearance of an effective immune response. It explains the fact that thymectomized animals are more susceptible to procedures that induce autoimmunity that euthymic animals [18].
It is also well known that if in a tolerant organism the number of APCs is increased enough, the immune system may turn on an immunological response. To study under which conditions our model reproduces this effect it is useful to understand that to change the parameters space from a state given by () to another characterized by () is equivalent, in terms of the deterministic equations, to change the state () to () along the straight line that matches them. So, the increase in the number of APCs () may be interpreted as a motion along a straight line with slope that passes through the origin in the vs map.
In our model, a tolerant organism is located in phases or . In the case of a system that starts in region , raising the number of APCs makes the system first, to evolve into region (see the line with points in Figure 3). Biologically, it may be interpreted as the entrance into a region where the tolerance disappears and the animals develops immunological response. On the other hand, when the organism starts in the phase , increasing the number of APCs drifts the system into a cyclic immune response (see the dashed line in Figure 3).
It is worth to highlight that our meanfield predictions compare very well with those obtained by K. León et al.[3] for a similar model of RE interactions regulated by APCs. However, our phase diagram is a bit richer. Exploring a parameter region where the bistable phase III may be now bounded by a phase where the biologically meaningful fixed points are unstable. Using parameters in this zone we have numerically solved the system of differential equations (24) and (25). As Figure 4 shows, in this zone appear very well defined closed orbits, representative of limit cycles.
The existence of this limit cycle allows the model to reproduce the emergence of a cyclic immune response in pathologies such as Multiple Sclerosis [19], [20]. In this way, changing the parameters one can observe two different kind of immune response, one stationary and another cyclic. Both patologies are observed in the Multiple Sclerosis, perhaps because genetic variations between individuals turn on the immune response within a different set of parameters. To our knowledge, no other mathematical model based on the dynamics of T cells regulated by R cells had predicted it together with effective response and tolerance.
Moreover, it is useful to compare the previous meanfield results with our original stochastic model. The reasons are twofold, first it should confirm that the meanfield approach is correct when studying large systems and second, it should shed some light on the behavior of smaller systems, like those characteristics of specific responses. To do this, we run Gillespie’s algorithm. In Figure 5 we show the results of typical runs for parameter values near the different fixed points, while in Figure 6 the numerical solutions of the deterministic equations for and are compared with the Gillespie’s simulations for the cyclic regime. These graphs demonstrate that, within a stochastic approach, strong fluctuations appear around the values predicted by the meanfield equations. Below, we show that if the size of the population is small enough these fluctuations become relevant and may qualitatively change the results of the mean field approximation.
a  b  c 
a  b  c 
4.2 The role of the fluctuations
In smaller systems the noise may induce qualitative changes in the mean field dynamics discussed above. This fact might be particularly relevant, since it could imply that the dynamics of small T cells clones would be significantly different from the one described here and in previous works [3]. We have carried out a series of Gillespie’s simulations with our stochastic model using parameters sets that belong to region III when classified according to the meanfield results. Our results illustrated that steady state stabilities, specially for the steady state dominated by regulatory T cells (tolerant steadystate), could be significantly affected by stochastic fluctuations. Figure 7 shows how in simulations the T cell concentrations remain oscillating around the meanfield value of the tolerant steadystate where the simulations began. However after a while fluctuations drive the system away into a state dominated by effector T cells (immune steady state). In other cases simulations shows that the tolerant steady state became unstable under stochastic fluctuations but it dynamically evolve into a trivial steady state with zero effector and regulatory T cells.
A new phase diagram (predicted with the help of (15) is represented in Figure 8, and has some differences when compared with Figure 3. The more striking feature of this map is the apparition of a new zone embedded in region with the characteristics of phase . In this region of the parameters, instead of having a cyclic behavior, the system evolves into a tolerant state (region ). This conclusion can be corroborated by simulations, and is in contradiction with what the normal analysis of eigenvalues predicts.
This result sustains the idea that, increasing the number of APCs, a cyclic immune response may be turned into a tolerant response (see the dashed line joining regions and in Figure 8). Along this line, when the number of APCs is increased, simulations show a stretching of the amplitude of the cycles until they reduce to a mere fluctuating behavior around the deterministic fixed point as in Figure 5 c.
On the other hand, fluctuations, also reduced the bistable region favoring the immune response of the system. Note however, how small spots of bistability persist within region .
Moreover, we study the effect of these fluctuations when the system is subject to external perturbations. For example, in the phase III, the system, depending on the initial conditions may stay in the tolerant state or in the state . In the mean field model, to switch the system from one state to the other, it must be strongly perturbed and moved away from the respective basin of attraction of the fixed point. In the case of small systems, small perturbations may switch the state of the system, provided they are frequent enough.
In figure 9 we show the probability () of switching from the tolerant state to the state as a function of the number of perturbations. In this case, with a period the number of cells at a given time was reduced an 80 percent. is the size of the simulations (different symbols mean three different ). As can be easily see if the number of perturbations is small enough , the system do not switch from one state to the other. It is always tolerant. On the other hand, when , i.e. when the period of the perturbations is small enough, the system goes to the responsive state with probability 1. The existence of a clear threshold for the minimum number of perturbations needed to change the state also signatures the importance of the fluctuations for the occurrence of this transition.
All these results confirm our previous intuition. The role of the fluctuations may be relevant and may alter significantly the predictions of the meanfield models. In particular, regions that were considered bistable may now develop immune response, and regions in which a cyclic pattern was predicted may become fully tolerant. Moreover, also the response of the system to external perturbations change. Perturbations, that may be produce only transient behaviours in the deterministic model, may become relevant when fluctuations are taken into account. Of course, the magnitude of these results are parameter dependent and probably richer diagrams may be found exploring further the model. But our results already prove that the interpretation of the specific immune response in base of mean field models should be done with care.
4.3 Power spectrum
Any single run of the stochastic system (as panels b and c in Figure 5) shows that coherent oscillations are sustained in parameter regimes in which the deterministic equations approach a fixed point. It evidences the fact that fluctuations can not be safely ignored in systems composed of a few thousands of elements.
Therefore, it is important to understand, not only the consequences of these fluctuations as done above, but also their properties. To do this, we analyze the characteristics of these fluctuations around the steady state, using equation (13):
(28) 
keeping in mind that for our model
and
Then, evaluating A and D in the corresponding fixed point, and using them to calculate eq. (28) we can determine the contribution of every frequency to the composition of the oscillations. Then, the discrete Fourier transform of a data of concentration and time that comes from the simulations can be compared with the power spectra curve calculated with (28). It would give us a way to know how good the made approximations were.
In Figure 10 we show the power spectra calculated by (28) and the one obtained by averaging the Fourier transforms of 6000 different Gillespie’s simulations. Both have been normalized in a way that and it is evident that the coincidence is remarkably good.
The peak in the power spectrum indicates the existence of a resonance frequency. In it, the fluctuations of the concentrations that define the state of the system amplify their values as a function of the parameters of the model. This effect may be interpreted as the induction by the internal noise of a a kind of stochastic resonance that is responsible for the emerging of amplified oscillations in systems which, when described through the deterministic rate equations for a given set of parameters lack a cyclic behavior.
From the technical point of view, this characteristic frequency arises because of the existence of complex eigenvalues with nonzero imaginary parts for the Jacobian at the fixed point. In the standard stability analysis, the imaginary part of the eigenvalues are responsible for oscillatory transient behaviors near the fixed points. However, here the stochasticity due to the smallness of the system results in persistent perturbations away from this fixed point, so that both features together result in an overall oscillatory effect.
The existence of this characteristic frequency could be useful to induce transition between two states of the system just by making a parameter (such like the number of APCs) to fluctuate with a small amplitude but with the appropriate frequency. Simulations corroborate this when it is induced a change from a tolerant phase to an immune one by making to fluctuate near a  region border. Work is in progress to check the reliability of this approach to develop new therapeutic techniques.
5 Conclusions
We present a stochastic model of interactions between Effector and Regulatory CD4+ T cells in the presence of Antigen Presenting Cells that is able to show active, tolerance, or cyclic immune responses. We characterize the phase diagram in the mean field approximation, and then analyze the corrections, due to finite size effects, to this phase diagram in the presence of fluctuations. This corrected phase diagram may be understood as a first attempt to characterize a stochastic model for the specific response of the immune system. In fact, we prove that the fluctuations may strongly alter the mean field predictions turning for example tolerant zones in responsive or cyclic responses into tolerant, and therefore that any analysis of the clonal response of the immune system based on mean field predictions must be taken with caution. Moreover, we show that our model present sustained oscillations at a characteristic frequency that may be relevant to understand, or treat these clonal responses.
Appendix A: Gillespie’s Algorithm
The occurrence of every process that can take place into our system has a very strong stochastic component. Every reaction depends on a lot of variables that can only be described in a probabilistic way, so in a given state, both the next process that would take place, and the moment in which it will occur are random variables. Even when it is not possible to find an analytic solution for the Master equation that governs the system, microscopic simulations of the processes defined by the relations (20) to (23) can be carried out using the algorithm originally proposed by Gillespie [9].
The algorithm keeps the random fashion in which processes occur. The backbone of such a simulation consists in determining which one will be the next process, and when will it occur. Every run of the simulation will give a different temporal sequence of concentration that is, according with the Master equation of the system.
The procedure can be synthesized in three main steps:

Calculate the occurrence probabilities () per unit time for every process that can take place in the system.

Using the and values obtained, increase the time by , and adjust the cellular population levels to reflect the occurrence of the process.
References
 [1] R. E. Langman and M. Cohn. The ET (elephanttadpole) paradox necessitates the concept of a unit of Bcell function: the protection. Mol. Immunol., 24(7):675–697, 1987.
 [2] S. Sakaguchi et al. Immunologic tolerance maintained by CD25+ CD4+ regulatory T cells: their common role in controlling autoimmunity, tumor immunity, and transplantation tolerance. Immunological Reviews, 182:18–32, 2001.
 [3] K. León et al. Modelling Tcellmediated suppression dependent on interactions in multicellular conjugates. Journal of Theoretical Biology, 207:231–254, 2000.
 [4] A. S. Perelson and G. Weisbuch. Immunology for physicists. Reviews of Modern Physics, 69(4):1219–1267, 1997.
 [5] J. Carneiro et al. When three is not a crowd: a crossregulation model of the dynamics and repertoire selection of regulatory CD4+ T cells. Immunological Reviews, 216:48–68, 2007.
 [6] J. Elf and M. Ehrenberg. Fast evaluation of fluctuations in biochemical networks with the linear noise approximation. Genome Research, 13:2475–2484, 2003.
 [7] L. E. Reichl. A Modern Course in Statistical Physics. 2nd ed. John Wiley and Sons, INC, 1998.
 [8] C. Gardiner. Handbook of stochastic methods. Springer Verlag, Berlin, 1985.
 [9] D. T. Gillespie. Exact stochastic simulation of coupled chemical reactions. The Journal of Physical Chemistry, 81:2340–2361, 1977.
 [10] N. G. van Kampen. Stochastic Processes in Physics and Chemistry. Elsevier, Amsterdam, 1981.
 [11] M. Scott, T. Hwa, and B. Ingalls. Deterministic characterization of stochastic genetic circuits. PNAS, 104:7402–7407, 2007.
 [12] H. Risken. The FokkerPlanck Equation: Methods of Solution and Applications. Springer Verlag, Berlin, 1989.
 [13] A. McKane, T. Newman, J. Nagy, and M. Stefanini. Amplified biochemical oscillations in cellular systems. Journal of Statistical Physics, 128(12):165–191, 2007.
 [14] S. H. Strogatz. Nonlinear Dynamics and Chaos: With applications to physics, biology, chemistry, and engineering. Addison Wesley, Cambridge, MA, 1994.
 [15] T. Takahashi, Y. Kuniyasu, M. Toda, N. Sakaguchi, M. Itoh, M. Iwata, J. Shimizu, and S. Sakaguchi. Immunologic selftolerance maintained by CD25+CD4+ naturally anergic and suppressive T cells: induction of autoimmune disease by breaking their anergic/suppressive state. Int. Immunology, 10:1969–1980, 1998.
 [16] A. M. Thornton and E. M. Shevach. CD4+ CD25+ immunoregulatory T cells suppress polyclonal T cell activation in vitro by inhibiting interleukin 2 production. The Journal of Experimental Medicine, 188:287–296, 1998.
 [17] K. León, K. García, J. Carneiro, and A. Lage. How regulatory CD25+CD4+T cells impinge on tumor immunobiology? On the existence of two alternative dynamical classes of tumors. Journal of Theoretical Biology, 247:122–137, 2007.
 [18] E. M. Shevach. Regulatory T cells in autoimmmunity. Annual Review of Immunology, 18:423–449, 2000.
 [19] A. Compston and A. Coles. The natural history of relapses in multiple sclerosis. Lancet, 359:1221–1231, 2002.
 [20] T. Vollmer. The natural history of relapses in multiple sclerosis. Jounal of the Neurological Sciences, 256 Suppl 1:S:5–13, 2007.