On the aging dynamics in an immune network model

arXiv:cond-mat/0305285v1 [cond-mat.stat-mech] 13 May 2003 On the Aging Dynamics in an Immune Network Model Mauro Copelli, Rita Maria Zorzenon dos Santos Laboratório de Fı́sica Teórica e Computacional Departamento de Fı́sica, Universidade Federal de Pernambuco Cidade Universitária, 50670-901, Recife, PE, Brazil Daniel Adrián Stariolo∗ Instituto de Fı́sica, Universidade Federal do Rio Grande do Sul CP 15051, 91501-970, Porto Alegre, Brazil November 21, 2018 Abstract Recently we have used a cellular automata model which describes the dynamics of a multi-connected network to reproduce the refractory behavior and aging effects obtained in immunization experiments performed with mice when subjected to multiple perturbations. In this paper we investigate the similarities between the aging dynamics observed in this multi-connected network and the one observed in glassy systems, by using the usual tools applied to analyze the latter. An interesting feature we show here, is that the model reproduces the biological aspects observed in the experiments during the long transient time it takes to reach the stationary state. Depending on the initial conditions, and without any perturbation, the system may reach one of a family of long-period attractors. The perturbations may drive the system from its natural attractor to other attractors of the same family. We discuss the different roles played by the small random perturbations (“noise”) and by the large periodic perturbations (“immunizations”). 1 Introduction In this paper we discuss a model for the evolution of the immune repertoire of B cells, which are responsible for the humoral immune responses. B cells belong to one of the main classes of white blood cells: the lymphocytes. These cells carry on their surface the order of 105 molecular receptors (proteins) and once activated they produce antibodies, which are copies of their molecular receptor. During the life time of an individual the immune system is able to produce the order of 1011 different antibodies or different populations of B cells. The antigen (virus, bacteria, poison, cellular residue, etc) is not recognized as a whole but by its epitopes, which are patches on its structure that may be recognized by specific sites of the antibody molecules. By pattern recognition different antibodies will mark the epitopes of a given antigen, therefore forming a complex that will be eliminated by macrophages (another class of white blood cells) [1, 2]. According to clonal selection theory [2, 1], elements that challenge the immune system will determine the populations (clones) of B cells that will proliferate: those populations will produce antibodies which will be able to recognize different epitopes of the specific antigen. The immune network theory [3, 2], however, is based on the fact that the antibodies (and molecular receptors) are able to recognize and to be recognized, and therefore during the immune response there are different types of interaction: antigen-antibodies and antibodies-B cells. In other words, when a given population of B cells is activated by the presence of a given antigen the produced antibodies will not only mark the specific antigens but also activate new B cell populations with complementary molecular receptors, which in turn will recognize them. The increase on the concentration of these complementary populations, on their turn, will maintain the proliferation of the antigen recognizing population, installing a feedback mechanism that will keep several populations activated. This kind of dynamics will generate a functional multi-connected network among different populations of B cells that will be dynamically regulated by mechanisms of activation and suppression. The network will then play an important role on the regulation of the immune responses. Although the immune network theory is part of the current immunological thinking, there are only few experiments supporting the interaction among clones with complementary receptors and the existence of such a network [4, 5]. According to these experimental findings, if the network exists only 10 − 20% of the populations will belong to it, the rest of the immunocompetent populations remaining at rest. Recently we have successfully used a mathematical model [16, 6, 7, 8] (inspired in a previous one proposed by de Boer, Segel and Perelson [17]) which takes into account the main features of Jerne’s immune network theory, ∗ Research Associate of the Abdus Salam International Center for Theoretical Physics, Strada Costiera 11, Trieste, Italy to simulate experiments on immunization and aging performed with mice [9] that could not be explained by the clonal selection theory . The simulations allowed to interpret the experimental results from the point of view of the immune network theory. The model allows to follow the evolution of the concentrations of the different populations of B cells in discrete shape space, a formalism which maps all possible molecular receptors of a given organism into points of a ddimensional space. To each point (receptor) we associate a clone that corresponds to the population of B cells and antibodies characterized by this molecular receptor. The concentration of each clone will be described by a three-state automaton representing low, medium and high concentrations, and the interaction among different clones is based on complementarity. As far as we know, the model [17, 16, 19, 6, 7] corresponds to the first successful attempt in describing the dynamics of the immune network as proposed by Jerne [3] and recently it could be used to reproduce experimental results performed with mice [9]. Besides the biological implications of the results obtained in Ref. [9], the dynamical behavior exhibited by this model in the biologically relevant parameter region is quite interesting by itself and should be better investigated. This region has been shown [6] to exist for dimensions d ≥ 2 and comprises a broad stripe near the transition between stable and chaotic behaviors, in which the model describes a multi-connected functional network [6, 7]. In this region we found the majority of the populations in the resting state (low concentration) while the activated ones may reach 10 − 20% of the total number of populations. The activated populations are aggregated in clusters of different sizes, which fuse and split as time passes following an aggregation and disaggregation dynamics. Therefore the largest cluster at each time step is found in a different region [7]. In the experiments described in ref. [9], 6 mice of the same linage are subjected to the immunization protocol as follows: the researchers inject ova by means of an intra-venal injection, wait for 14 days and measure the number of specific antibodies. Then they inject again the same amount of antigen, wait for 7 days and measure the amount of antibodies, and continue by repeating the same protocol every 7 days. In order to simulate the immunization protocol, the CA is subjected to specific perturbations by flipping chosen resting populations (low concentration) to their activated state (high concentration). Depending on the perturbation (damage) size, it may disappear after a few time steps or part of the damage (activated populations) may be incorporated to the network. In other words, the memory about the perturbations is due to the ability of the system to adapt (plasticity) and incorporate information about them. We have used two types of perturbations: random small ones, which correspond to the noise that mice are subjected to during the experiment, caused by the environment, and periodic large ones that will simulate the immunization protocol of multiple antigen perturbations under the same conditions [9]. After few presentations there is a saturation of the response of the system to the perturbation. During this process, the system incorporates new information and some of the previous information is lost, keeping the number of activated sites in the network almost unchanged. This kind of behavior was also observed experimentally in mice [10], where saturation is related to a refractory behavior of the immune system. We have also observed aging effects on the dynamics of this system [9]. An older system is more rigid: the network loses plasticity to incorporate new information [9, 11, 12]. A recent study [13] has shown that the distribution of cluster sizes during the time evolution of the system has a characteristic cluster size (exponential behavior), but the distribution of persistence times (the period during which a given population remains activated and belongs to the network) exhibits a power law behavior. While the existence of a characteristic cluster size may be related to the loss of plasticity, the power law behavior of the persistence time may be associated to the memory generated by the dynamics of the system. The question we address here refers to how the learning process takes place dynamically and what is the cause of the loss of plasticity. The slow dynamics observed in this system presents analogies with the physical aging effects observed and reported on glass studies [14, 15]: as the system gets older (ages) there is a loss of plasticity for structural or molecular relaxation and less changes are observed during the relaxation time. The mechanisms underlying the slow dynamics of glassy systems are the spatial (or geometric) disorder related with the difficulty to satisfy simultaneously all microscopic interactions, a characteristic called frustration. Glasses and spin glasses exhibit a rough energy landscape with many local minima which is responsible for slowing down the relaxation towards equilibrium. Their dynamics depends on the past history of the sample and under small perturbations they relax slowly. In our biologically motivated model the non-local interactions are based on complementarity (both perfect and slightly defective matches) and reflect the activation and suppression mechanisms in the complementary regions of the shape space. The analogue to frustration, in the network, is generated by the inability of the system to incorporate information and to satisfy the constraints of the activation and suppression mechanisms. We show that the CA dynamics gives rise to a large family of periodic attractors which are very robust, and could be regarded as the analogues of the minima in the glassy systems. The lost of plasticity and aging effects will therefore be related to the non-ergodicity in the phase space. In this paper the dynamics of the immunological responses in the network model are investigated using the common tools adopted in the study of glassy systems [15]. In particular we will focus on the relaxation of twotime autocorrelations of the B-cell populations, the structure of the attractors and effects of perturbations and immunizations. The paper is organized as follows: in section 2 we describe the model and the procedures adopted in order to measure the correlations among different configurations of the system; in section 3 we present and discuss the results obtained for small and large perturbations and in section 4 we present a summary, prospectives for future work and some conclusions. 2 The Model The model under study here is a modified version [6, 7, 8] of the model proposed by Stauffer and Weisbuch [16], which in turn was inspired in a previous model proposed by de Boer, Segel and Perelson [17] to describe the time evolution of the immune repertoire. It is a deterministic window cellular automaton model based on the shape space formalism [18], which describes the interactions of B-lymphocytes and antibodies, and the main mechanism underlying these interactions, which is pattern recognition (lock-key interaction). The dynamics of the system is regulated by specific interactions involving complementary molecular receptors of the different clones. The memory about the relevant antigens, presented to the system during its past history, emerges from the dynamics of the system, rather than being stored in a static registry. To each point of a d-dimensional discrete lattice we associate a given receptor, which in turn will be described by d-coordinates representing important physical-chemical features of the receptor that differentiate one from the other [18]. Since clones are classified according to their molecular receptor, to each point ~r of the discrete shape space we associate a three-state automaton B(~r, t) that will describe the concentration of its population over the time: low (B(~r, t) = 0), intermediate (B(~r, t) = 1) and high (B(~r, t) = 2). The time evolution of the cellular automaton is based in a non-local rule: population B(~r, t) at site ~r is influenced by the populations at site −~r (its mirror image or complementary shape) and its nearest-neighbors (−~r + δ~r) (representing defective lock-key interactions). The influence on the population at site ~r due to its complementary populations is described by the field h(~r, t): h(~r, t) = X B(~r ′ , t) (1) ~ r ′ ∈(−~ r+δ~ r) where for a given ~r the sum runs over the complementary shape ~r ′ = −~r and its nearest neighbors. Due to the finite number of states of the population B, the maximum value of the field h(~r) is hmax = 2(2d + 1). The updating rule is based on a window of activation, which describes the dose response function involved in B-cell activation [6, 16, 17]. There is a minimum field necessary to activate the proliferation of the receptor populations (θ1 ), but for high doses of activation (greater than θ2 ) the proliferation is suppressed. The updating rule may be summarized as: B(~r, t + 1) =  B(~r, t) + 1 B(~r, t) − 1 if θ1 ≤ h(~r, t) ≤ θ2 otherwise (2) but no change is made if it would lead to B = −1 or B = 3. We define the densities of sites in state i at time t as Bi (t) (i = 1, 2, 3). The initial configurations are randomly generated according to the following concentrations: B1 (0) = B2 (0) = x/2, while the remaining Ld (1 − x) sites are initiated with B(~r, 0) = 0. This model may exhibit stable or chaotic behaviors depending on the values of x, θ1 and θ2 − θ1 . However it is on the transition region between the two behaviors that the model behaves like a multi-connected network [6]. In order to simulate the immunization protocol performed in the mice experiments we have followed the procedure which is described in Ref. [9] and summarized below. We have adopted the scale of 5 time steps corresponding to 1 day [9]. While the system evolves according to the deterministic dynamics (eq. 2), small and large perturbations can be produced, by setting the state of the chosen sites at B(~r, t) = 2. Small Perturbations The small perturbations account for the immunological stimuli (noise) coming from the environment. The time interval between two consecutive small perturbations is a random number uniformly distributed between 1 and 100 time steps. Each perturbation corresponds to a random number of damages (from 1 to 3) introduced at regions of resting populations (B = 0) which are randomly drawn (at every perturbation). The size of each damage may also vary randomly from 1 to 3 (onion-like) concentric layers around a central site (containing 7, 25 or 63 populations, respectively, in 3 dimensions). Large Perturbations The large perturbations correspond to the immunization protocol which starts at a predetermined age of the mice. Like in the experiments, we stimulate the system periodically (every 35 steps ≃ 1 week), and always in the same region (which is initially chosen at random but kept unchanged along the simulation). The damage size in this case corresponds to six layers (377 populations) around an specific site. Previous results on this model have shown [9] that the response to the immunizations presents a strong dependence on the initial time (“age”) at which the periodic protocol starts (fitting experimental data extremely well for mice whose immunization protocol started at different ages). This has led us to the conjecture that the dynamical behavior of the system might be at least qualitatively similar to that of some glassy systems, despite the non-Hamiltonian nature of the CA dynamics. One of the quantities commonly used in the study of glassy systems is the two-time autocorrelation function between the system configurations at two given times t and tw . A common experiment in glassy systems consists in preparing the system at a high temperature and suddenly making a quench to a low temperature. Then the system is allowed to relax up to a waiting time tw , whose configuration is recorded. As the system continues to relax the autocorrelations between the instantaneous configurations at time t > tw and that at time tw are computed. The waiting time tw is called the age of the system in context of glassy systems. The monitoring of the two-time autocorrelations gives important insights on the relaxation process. The dynamics, whether stationary or not, can be readily recognized on the tw dependence of the autocorrelations, since in a stationary process two-time quantities depend only on time differences. Consequently by plotting correlations as a function of time difference t − tw for different waiting times it is possible to distinguish between an essential out of equilibrium process from a stationary one. The aging processes observed in glassy systems are then related to the lack of temporal invariance [15]. Inspired on this approach, we will analyze the multi-connected network dynamics by defining and analyzing quantities analogous to the two-time autocorrelations for the CA: Ctot (t, tw ) = N 1 X δ (B(~r, tw ), B(~r, t)) N (3) ~ r C22 (t, tw ) = 1 B2 (tw )N N X δ (B(~r, t), 2) , (4) ~ r|B(~ r ,tw )=2 where N = Ld is the total number of populations in the system. Ctot and C22 amount to normalized proximities (using Hamming distance as a measure) and from now on will be referred to simply as correlation functions. These quantities will be analyzed for different protocols: without any perturbation, with small perturbations (noise) and with large perturbations (immunization protocol) in order to differentiate the effects of the different kinds of perturbations. A complementary view of the long-term behavior of the system can be obtained by looking at the attractors to which the system evolves. In the present case it will consist mainly of cycles, which will be reflected by the periodicity of the correlations. For this purpose, we have obtained the return maps of the consecutive maxima of the correlation functions. Note that the period of the maxima thus obtained does not correspond to the period of the real system, which is at least twice as long [19]. 3 3.1 Results Without perturbations In previous works [19, 6, 7] it has been shown that depending on the initial concentration x of active populations the system may exhibit periodic or chaotic behavior. Systems with low initial concentrations of active sites (x < xc ) evolve to either fixed points or orbits with short periods, while for x > xc chaotic attractors appear. However, the biologically relevant region is in the transition region between these two behaviors, where the system reaches one of several periodic orbits (as will be seen below) with a very long period and after a long transient. From now on, all the results have been obtained using the same parameters adopted in Ref. [9]: d = 3, θ1 = hmax /3, θ2 = 2hmax /3 and x = 0.26 (on the transition region). Without any perturbation, the system evolves after a transient time towards a cycle with a large period, as shown in the return map of Fig. 1. We have also varied the waiting time from 10 to 100000 (not shown). The greater the waiting time, the closer to the attractor the system is, which is revealed by larger values of the autocorrelation functions. Once tw is larger than the transient, the time series for Ctot (t, tw ) (and also C22 (t, tw ) will include unity (see Fig. 1) and will not change for larger values of tw . A typical result for the time evolution of the correlation functions for tw = 100 is shown in Fig. 2. The upper panel corresponds to the time evolution of the densities B1 (intermediate concentrations) and B2 (high concentrations), while the lower panel corresponds to Ctot (t, tw ) (open circles) and C22 (t, tw ) (filled circles). Notice that while the concentrations relax to approximately constant values in a short time, Ctot and C22 take much longer to reach their attractors (note the logarithmic time scale). In order to study how the system relaxes towards the attractor it is more convenient to make use of C22 , since it measures the changes in the (more relevant) network of activated populations. In the case shown in Fig. 2 the transient time needed to attain the attractor is O(104 ) steps. It is necessary to point out, however, that this typical relaxation time is important only for the physical aspects of the dynamics. When mapped into the biological problem, it would correspond to ∼ 5.5 years, which is much longer than the average life time of the mice used in the kind of experiment we simulated. Therefore the relevant behavior, from the biological point of view, happens to be in the transient of the model and not in its stationary state, a result which is interesting on its own. 3.2 Random Small Perturbations How does the behavior of the system change when random small perturbations are produced on the parameter region used to simulate the real experiments performed with mice [9]? In order to investigate this issue, only the small perturbations described in Section 2 were produced, starting at time zero. In Fig. 3 we compare the time evolution of the densities and correlation functions obtained for single runs for the system without perturbations and with small random perturbations. Note that in the case with small perturbations the correlations initially follow those of the purely deterministic system, showing only small differences. After about 103 − 104 steps, however, they start decreasing faster, indicating some sort of cumulative effect that drives the system away from the region in phase space that it had approached until tw . These effects are more easily 1.0005 1 Ctot(t+1) 0.9995 0.999 0.9985 0.998 0.9975 0.9975 0.998 0.9985 0.999 Ctot(t) 0.9995 1 1.0005 Figure 1: Return map for maxima of the total correlation function Ctot (t) for L = 50 and tw = 10000 (single run). 0.03 B1 B2 B1, B2 0.025 0.02 0.015 0.01 0.005 0 1 10 100 1000 10000 100000 1e+06 10000 100000 1e+06 t 1 0.95 Ctot 0.9 0.85 0.8 0.75 Ctot c22 0.7 1 10 100 1000 t-tw Figure 2: Densities of intermediate and high concentrations vs. time (upper panel) and autocorrelations vs. t − tw for L = 50 and tw = 100, without any perturbation. noted for C22 . Moreover, the perturbations do not alter the behavior of the densities, as expected, since the number of activated populations is kept approximately constant by the self-regulatory mechanisms embedded in the dynamical rules. The changes are observed only in the populations that belong to the active network: in order to incorporate new information (new populations), older ones are deactivated. From the biological point of view, the realizations of the perturbations (small and large) differ from one individual to another, building the identity of each individual. At the end of life each individual will have a different history in terms of perturbations (antigen presentations) translated into the configuration of the populations belonging to the active network. This is nicely illustrated by Fig. 4, where two initially identical copies of a system undergo different realizations of the small perturbations according to the protocol described in section 2. The Hamming distance between them grows on a long time scale, revealing the mechanisms behind the behavior of the correlations in Fig. 3. Returning to Fig. 3, it is important to stress that the changes observed in the concentrations for very long times (upper panel) are due to finite size effects. They are caused by the fact that all small perturbations are produced on regions of resting populations. For a finite system, after a long time all the possibilities will have already been explored. Increasing the size of the system, the changes on the densities disappear. This is shown in Fig. 5, where we repeat the simulations, under the same conditions of Fig. 3, for a larger system (L = 100). 0.03 B1 B2 p33 B1 p33 B2 0.025 B1, B2 0.02 0.015 0.01 0.005 0 Ctot 1 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 10 100 1000 t 10000 100000 1e+06 1000 t-tw 10000 100000 1e+06 Ctot c22 p33 ctot p33 c22 1 10 100 Figure 3: Densities vs. time (upper panel) and correlations vs. t − tw for L = 50 and tw = 100; without perturbation (circles) and with small perturbations (triangles). 0.1 0.09 0.08 0.07 0.06 HD 0.05 0.04 0.03 0.02 0.01 0 1 10 100 1000 10000 100000 1e+06 t Figure 4: Normalized Hamming distance between two initially identical configurations subjected to different sequences of small perturbations (three samples). In order to test out the robustness of these results, we have varied the parameters controlling the protocol of the random small perturbations. For instance, by increasing the maximum number of different perturbations from 3 to 6 at each presentation, and/or by changing the maximum time interval between consecutive perturbations from 100 to 10, we observe the same qualitative behavior, with a faster decrease of the correlations — as expected, since in both cases the noise has been increased. What happens when the system’s autocorrelations decrease due to the small random perturbations after it has reached a periodic attractor? Is the system driven to another cycle? To answer this question we have studied the stability of the cycles using the following procedure: we let the system evolve without perturbation towards its attractor for tw = 10000 time steps, after which we perturb the system with random small perturbations during a time interval ∆t1 . Then we turn off the perturbations and allow the system to relax during another time interval ∆t2 after which we obtain the return map of the correlations in the following 200 time steps. If we produce only one perturbation at tw (∆t1 = 1) the system remains in its original cycle, yielding a return map similar to that of 0.25 B1 B2 p33 B1 p33 B2 B1, B2 0.2 0.15 0.1 0.05 0 1 10 100 1000 t 10000 100000 1e+06 10 100 1000 t-tw 10000 100000 1e+06 1 0.95 Ctot 0.9 0.85 0.8 Ctot c22 p33 ctot p33 c22 0.75 0.7 0.65 1 Figure 5: Densities vs. time (upper panel) and autocorrelations vs. t − tw for L = 100 and tw = 100 (the same conditions used in Fig. 3); without perturbation (circles) and with small perturbations (triangles). 1 0.9995 Ctot(t+1) 0.999 0.9985 0.998 0.9975 0.997 0.997 0.9975 0.998 0.9985 Ctot(t) 0.999 0.9995 1 Figure 6: Return map for maxima of total correlation function, for L = 50, tw = 10000, ∆t1 = 1000 and ∆t2 = 10000. Fig. 1. Then, under the same conditions we repeat the simulations adopting now ∆t1 = 1000 and ∆t2 = 10000. We observe in Fig. 6 that the return map has changed slightly when compared to Fig. 1. In particular, it no longer shows Ctot = 1 in the time series, but remains periodic, signaling that the system has shifted to a different cycle due to the perturbations. Apparently the periods observed in Figs 1 and 6 are the same. The distribution of periods and transient times are currently under investigation, results will be published elsewhere. Fig. 6 remains the same by increasing ∆t2 , which guarantees that the differences with respect to Fig. 1 are not a transient effect. According to these results, the role the noise plays, if allowed enough time to perturb the system significantly, is to drive the system from one attractor to a nearby one, which suggests that there is a family of periodic attractors which can be very close to one another in phase space. These effects, however, take place on a time scale which is longer [O(103 − 104 ) time steps] than the lifetime of the mice [O(102 − 103 ) time steps]. The small perturbations are therefore of little importance to the CA dynamics in the biologically relevant time scale, as suggested in previous work [9]. The behavior of the autocorrelation functions in Figs. 3 and 5 is reminiscent of what is observed in glassy systems. Here the CA approaches a periodic attractor, being thereafter deflected by the small perturbations. In 1 0.95 C(t, tw) 0.9 0.85 0.8 0.75 0.7 0.65 1 10 100 1000 t-tw 10000 100000 1e+06 Figure 7: C22 (t + tw ; tw ) for different waiting times before small perturbations. From bottom to top, tw = 100, 5000, 10000. glassy systems, temperature drives the system from one local minimum to another, preventing it from getting stuck in local minima of the potential energy surface. As the physical age of the system (characterized by the value of the waiting time tw ) grows, it finds itself exploring deeper and deeper regions of a rough potential energy surface, diffusing towards equilibrium [20]. Due to the roughness of the potential energy, equilibrium is only attained on very long time scales and relaxation is very slow. The older the system is the more it gets confined to a restricted region of phase space and the time scales for relaxation get longer and longer. Eventually, if we wait long enough, the system equilibrates and the dynamics becomes stationary, losing sensitivity to the waiting time. This picture of an aging physical system is reminiscent to the loss of plasticity for adaptation in a living organism as it gets older. In either case, one observes a strong dependence on the waiting time tw , evident when measuring two-time quantities like autocorrelation functions and responses. Results for the CA model are shown in Fig. 7, where we see the decay of autocorrelations for three systems with different ages or waiting times. Note that the horizontal axis is the time difference between the total time and the waiting time. The three curves should collapse in the case that the dynamical evolution is stationary. 3.3 Large immunizations What is the role of the large perturbations on the dynamics of the system? From previous studies we would expect the large perturbations to accelerate the aging process: while the random small perturbations would change the route to the natural attractor of the system, the large ones would reduce the transient time, a conjecture that would explain the loss of plasticity of the older mice [9]. The protocol adopted is the one described in Section 2, with six-layer perturbations every 35 time steps, always at the same sites. In Fig. 8 we compare the results obtained for the unperturbed system and those of the system subjected only to large immunizations starting at tw (note that the densities are now also plotted as functions of t − tw ). Somewhat surprisingly, the decrease of the correlations for large perturbations is small, when compared to the case of the small perturbations (compare with Fig. 3). In hindsight, however, this can be understood because the large perturbations are always produced in the same region. For very long times the correlations will attain a stationary regime whose active populations contain at least part of the region involved in the immunizations [9]. In other words, the driving produced by the large periodic perturbations is of a completely different nature than that of the small perturbations. The large perturbations seem to play a selective role: the cycles that the system can reach are restricted to those that contain at least part of the populations incorporated during the immunization. According to this picture, the aging effects observed in Ref. [9] could be simply related to the exploration of the phase space: the older the system the less possibilities of choosing new cycles it would have. In Section 4 we discuss this issue further. In Fig. 9 we compare the time evolution of the densities B1 and B2 and the autocorrelation functions for the system subjected only to large perturbations and for the system with both small perturbations and large immunizations. The results, as expected, confirm the dominance of the small perturbations, over the large ones, in driving the system faster to a different attractor. The increase of the densities around t ∼ 104 for L = 50 B1, B2 0.065 0.06 0.055 0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 B1 (p0 P16) B2 (p0 P16) B1 (p0 P00) B2 (p0 P00) Ctot, c22 1 10 100 1000 10000 100000 1000 10000 100000 t - tw 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 Ctot (p0 P16) c22 (p0 P16) ctot (p0 P00) c22 (p0 P00) 1 10 100 t-tw Figure 8: Densities vs. time (upper panel) and autocorrelations vs. t − tw for a system with large immunizations, L = 50 and tw = 100. Triangles: without perturbations of any kind. Circles: both large immunizations and small perturbations. 0.07 B1, B2 0.06 0.05 B1 (p33 P16) B2 (p33 P16) B1 (p0 P16) B2 (p0 P16) 0.04 0.03 0.02 0.01 ctot, c22 1 10 100 1000 10000 100000 1000 10000 100000 t - tw 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 ctot (p33 P16) c22 (p33 P16) ctot (p0 P16) c22 (p0 P16) 1 10 100 t-tw Figure 9: Densities of active populations vs. time (upper panel) and autocorrelations vs. t − tw for tw = 100 and L = 50. Triangles: large immunizations only. Circles: both large immunizations and small perturbations. corresponds to the same finite size effects occurring in Fig. 3, being associated to the active populations which have been incorporated into the network by means of the immunization protocol. 4 Concluding Remarks We have analyzed the dynamics of a cellular automata model for the B-cell repertoire which has reproduced experimental results of immunizations on mice. Our analyses provide a broader context in which memory and plasticity take place, as discussed in Ref. [9]. 55000 50000 18500 18400 18300 18200 18100 18000 17900 17800 17700 100 45000 HD 40000 35000 30000 25000 1000 10000 100000 20000 15000 1 10 100 1000 10000 100000 t Figure 10: Hamming distance between 100 pairs of random configurations as a function of time (average and standard deviation). Inset: evolution for t > 100 shows that the HD attains a stationary regime for sufficiently long times (standard deviation not shown). Defining quantities analogous to the autocorrelation functions used in glassy systems, we have shown that the biologically relevant phenomena takes place in the transient regime of the model. The dependence of these functions on both t and tw (as opposed to the difference t−tw only) is reminiscent of the “aging” behavior observed in glassy systems, despite the fact that the underlying dynamics of both systems are controlled by completely different mechanisms. Starting from different initial conditions the deterministic dynamics takes the system towards a family of long-period attractors. When subjected to random small perturbations (Fig 4), the system is driven towards a new attractor of the family, revealing that most of the noise is assimilated. However, only part of the large perturbations is incorporated, due to the mechanisms of activation and suppression, leading to a saturation of the learning process [9]. From the biological point of view, the history of the mouse (sample) will be written by the different antigen presentations (random perturbations), starting from its “birth” (initial condition). Since the system is large but finite there is a maximum amount of information it may incorporate. The closer it is to its “destiny” attractor, the less information it is able to learn, since the deeper it already is in a given basin of attraction. Therefore the biological aging corroborated by the experimental results may be simply a consequence of this dynamical feature. From the results obtained up to now there are evidences that the purely deterministic dynamics is therefore non-ergodic. Despite their large number, however, the family of periodic attractors occupy only a fraction of the phase space. Evidences of this compression in phase space has been obtained in the following computer experiment: selecting randomly two initial configurations (with the same initial concentration), we measured the Hamming distance between them as a function of time as both systems evolve without any perturbation. In Fig. 10 we show the evolution of the average HD for 100 pairs (error bars are standard deviations). During the first 10 time steps the HD decreases very rapidly, but for t > 100 we observe a quasi-stationary regime, reflecting the slow driving to the attractors (see inset). Note that this behavior is somewhat similar to that of the autocorrelation functions for small tw . Still focusing on the evidence of phase space compression, in Fig. 11 we present the distribution of HD sampled from 500(500 − 1)/2 pairs, for two consecutive time steps (t = 4000 and t = 4001). The distributions can be well described by a Gaussian, with a width that remains approximately constant (even for long times). Notice that the average value for t = 4001 is slightly larger than for t = 4000, reflecting the oscillatory behavior of the average HD as the pair of samples reaches their periodic attractors. It should be noted that, for large N , the Central Limit Theorem assures that randomly chosen initial configurations (with the same x) naturally give rise to a Gaussian distribution of the HD between any two of them, at time zero. Interestingly, the CA dynamics does not change the shape of the distribution, its only effect being to essentially shift the Gaussian towards lower mean values, on a long time scale. The spatio-temporal structure of the cycles, as well as their transients and basins of attraction, would be the object of further study and will be published elsewhere. This work was partially supported by CNPq, Faperj, Finep, Capes and Projeto Enxoval-UFPE. MC and RMZS acknowledge IF-UFF, where part of this work was done during their stay there. 90 80 No. of occurrences 70 60 50 40 30 20 10 0 15000 16000 17000 18000 19000 20000 21000 22000 23000 Hamming distance 80 t+1 No. of occurrences 70 60 50 40 30 20 10 15000 16000 17000 18000 19000 20000 21000 22000 23000 Hamming distance Figure 11: Distribution of Hamming distance sampled over 500(500-1)/2 configuration pairs at time t = 4000 and t = 4001. Vertical lines show the average value. References [1] Janeway, C. A, Travers, M. W. and Capra, J., Imunobiology: The Immune System in Health and Disease, 4th ed., Elsevier Science Ltda. Garland Publishing, New York, 2000. [2] Vaz, N. M. and de Faria, A. M. C., Guia Completo de Imunobiologia, Coopmed Editora, Belo Horizonte, Brazil (1993). [3] Jerne, N.K., Ann. Immuno. (Inst. Pasteur) 125 C, 373-389 (1974). [4] Holmberg, D., Anderson, Å, Carlsson, L. & Forsgren, S., Immunol. Rev. 110, 84-103 (1989). [5] Coutinho, A., Immunol. Rev. 110, 63-87 (1989). [6] Zorzenon dos Santos, R. M. and Bernardes, A. T., Physica A 219, 1-12 (1995). [7] Bernardes, A. T. and Zorzenon dos Santos, R. M., J. Theor. Biol. 186 173 (1997). [8] Zorzenon dos Santos, R. M., in Annual Reviews of Computational Physics VI, 159-202, ed. by D. Stauffer, World Scientific (1999). [9] Zorzenon dos Santos, R. M. and Bernardes, A. T., Phys. Rev. Lett. 81, 3034 (1998). [10] Verdolin, B., MSc. Thesis, Departamento de Bioquı́mica e Imunologia, Instituto de Ciências Biológicas, UFMG, Belo Horizonte (1997); Faria, A.M.C. et al., Braz. J. Med. Biol. Res. 31, 35-48 (1998); Vaz, N. et al., Scand. J. Immunol. 46, 225-229 (1997). [11] Faria, A. M. C. et al., Mech. Ageing Dev. 102, 67-80 (1998). [12] Lahmann, W. M., Menezes, J. S., Verdolin, B. A. and Vaz, N. M., Braz. J. Med. Biol. Res. 25 813-821 (1992). [13] Bernardes, A. T. and Zorzenon dos Santos, R. M., Int. J. Mod. Phys. C 12, 1 (2001). [14] Struick, L. C. E., Physical aging in amorphous polymers and other materials, (Elsevier, Houston, 1976). [15] Bouchaud, J.-P., Cugliandolo, L. F., Kurchan, J. and Mézard, M., in Spin Glasses and Random Fields, A.P. Young editor, (World Scientific, Singapore, 1998). [16] Stauffer, D. and Weisbuch, G., Physica A 180, 42-52 (1992). [17] De Boer, R. J., Segel, L. A. and Perelson, A. S., J. Theor. Biol. 155, 295-333 (1992). [18] Perelson, A. S. and Oster, G. F., J. Theor. Biol. 81, 645-670 (1979). [19] Zorzenon dos Santos, R. M., Physica A196, 12 (1993). [20] Debenedetti, P. G. and Stillinger, F. H., Nature 410, 259-267 (2001). On the Aging Dynami s in an Immune Network Model Mauro Copelli, Rita Maria Zorzenon dos Santos Laboratorio de Fsi a Teori a e Computa ional Departamento de Fsi a, Universidade Federal de Pernambu o Cidade Universitaria, 50670-901, Re ife, PE, Brazil Daniel Adrian Stariolo Instituto de Fsi a, Universidade Federal do Rio Grande do Sul CP 15051, 91501-970, Porto Alegre, Brazil May 13, 2003 Abstra t Re ently we have used a ellular automata model whi h des ribes the dynami s of a multi- onne ted network to reprodu e the refra tory behavior and aging e e ts obtained in immunization experiments performed with mi e when subje ted to multiple perturbations. In this paper we investigate the similarities between the aging dynami s observed in this multi- onne ted network and the one observed in glassy systems, by using the usual tools applied to analyze the latter. An interesting feature we show here, is that the model reprodu es the biologi al aspe ts observed in the experiments during the long transient time it takes to rea h the stationary state. Depending on the initial onditions, and without any perturbation, the system may rea h one of a family of long-period attra tors. The perturbations may drive the system from its natural attra tor to other attra tors of the same family. We dis uss the di erent roles played by the small random perturbations (\noise") and by the large periodi perturbations (\immunizations"). 1 Introdu tion In this paper we dis uss a model for the evolution of the immune repertoire of B ells, whi h are responsible for the humoral immune responses. B ells belong to one of the main lasses of white blood ells: the lympho ytes. These ells arry on their surfa e the order of 105 mole ular re eptors (proteins) and on e a tivated they produ e antibodies, whi h are opies of their mole ular re eptor. During the life time of an individual the immune system is able to produ e the order of 1011 di erent antibodies or di erent populations of B ells. The antigen (virus, ba teria, poison, ellular residue, et ) is not re ognized as a whole but by its epitopes, whi h are pat hes on its stru ture that may be re ognized by spe i sites of the antibody mole ules. By pattern re ognition di erent antibodies will mark the epitopes of a given antigen, therefore forming a omplex that will be eliminated by ma rophages (another lass of white blood ells) [1, 2℄. A ording to lonal sele tion theory [2, 1℄, elements that hallenge the immune system will determine the populations ( lones) of B ells that will proliferate: those populations will produ e antibodies whi h will be able to re ognize di erent epitopes of the spe i antigen. The immune network theory [3, 2℄, however, is based on the fa t that the antibodies (and mole ular re eptors) are able to re ognize and to be re ognized, and therefore during the immune response there are di erent types of intera tion: antigen-antibodies and antibodies-B ells. In other words, when a given population of B ells is a tivated by the presen e of a given antigen the produ ed antibodies will not only mark the spe i antigens but also a tivate new B ell populations with omplementary mole ular re eptors, whi h in turn will re ognize them. The in rease on the on entration of these omplementary populations, on their turn, will maintain the proliferation of the antigen re ognizing population, installing a feedba k me hanism that will keep several populations a tivated. This kind of dynami s will generate a fun tional multi- onne ted network among di erent populations of B ells that will be dynami ally regulated by me hanisms of a tivation and suppression. The network will then play an important role on the regulation of the immune responses. Although the immune network theory is part of the urrent immunologi al thinking, there are only few experiments supporting the intera tion among lones with omplementary re eptors and the existen e of su h a network [4, 5℄. A ording to these experimental ndings, if the network exists only 10 20% of the populations will belong to it, the rest of the immuno ompetent populations remaining at rest. Re ently we have su essfully used a mathemati al model [16, 6, 7, 8℄ (inspired in a previous one proposed by de Boer, Segel and Perelson [17℄) whi h takes into a ount the main features of Jerne's immune network theory,  Resear h Asso iate of the Abdus Salam International Center for Theoreti al Physi s, Strada Costiera 11, Trieste, Italy 1 to simulate experiments on immunization and aging performed with mi e [9℄ that ould not be explained by the lonal sele tion theory . The simulations allowed to interpret the experimental results from the point of view of the immune network theory. The model allows to follow the evolution of the on entrations of the di erent populations of B ells in dis rete shape spa e, a formalism whi h maps all possible mole ular re eptors of a given organism into points of a ddimensional spa e. To ea h point (re eptor) we asso iate a lone that orresponds to the population of B ells and antibodies hara terized by this mole ular re eptor. The on entration of ea h lone will be des ribed by a three-state automaton representing low, medium and high on entrations, and the intera tion among di erent lones is based on omplementarity. As far as we know, the model [17, 16, 19, 6, 7℄ orresponds to the rst su essful attempt in des ribing the dynami s of the immune network as proposed by Jerne [3℄ and re ently it ould be used to reprodu e experimental results performed with mi e [9℄. Besides the biologi al impli ations of the results obtained in Ref. [9℄, the dynami al behavior exhibited by this model in the biologi ally relevant parameter region is quite interesting by itself and should be better investigated. This region has been shown [6℄ to exist for dimensions d  2 and omprises a broad stripe near the transition between stable and haoti behaviors, in whi h the model des ribes a multi- onne ted fun tional network [6, 7℄. In this region we found the majority of the populations in the resting state (low on entration) while the a tivated ones may rea h 10 20% of the total number of populations. The a tivated populations are aggregated in lusters of di erent sizes, whi h fuse and split as time passes following an aggregation and disaggregation dynami s. Therefore the largest luster at ea h time step is found in a di erent region [7℄. In the experiments des ribed in ref. [9℄, 6 mi e of the same linage are subje ted to the immunization proto ol as follows: the resear hers inje t ova by means of an intra-venal inje tion, wait for 14 days and measure the number of spe i antibodies. Then they inje t again the same amount of antigen, wait for 7 days and measure the amount of antibodies, and ontinue by repeating the same proto ol every 7 days. In order to simulate the immunization proto ol, the CA is subje ted to spe i perturbations by ipping hosen resting populations (low on entration) to their a tivated state (high on entration). Depending on the perturbation (damage) size, it may disappear after a few time steps or part of the damage (a tivated populations) may be in orporated to the network. In other words, the memory about the perturbations is due to the ability of the system to adapt (plasti ity) and in orporate information about them. We have used two types of perturbations: random small ones, whi h orrespond to the noise that mi e are subje ted to during the experiment, aused by the environment, and periodi large ones that will simulate the immunization proto ol of multiple antigen perturbations under the same onditions [9℄. After few presentations there is a saturation of the response of the system to the perturbation. During this pro ess, the system in orporates new information and some of the previous information is lost, keeping the number of a tivated sites in the network almost un hanged. This kind of behavior was also observed experimentally in mi e [10℄, where saturation is related to a refra tory behavior of the immune system. We have also observed aging e e ts on the dynami s of this system [9℄. An older system is more rigid: the network loses plasti ity to in orporate new information [9, 11, 12℄. A re ent study [13℄ has shown that the distribution of luster sizes during the time evolution of the system has a hara teristi luster size (exponential behavior), but the distribution of persisten e times (the period during whi h a given population remains a tivated and belongs to the network) exhibits a power law behavior. While the existen e of a hara teristi luster size may be related to the loss of plasti ity, the power law behavior of the persisten e time may be asso iated to the memory generated by the dynami s of the system. The question we address here refers to how the learning pro ess takes pla e dynami ally and what is the ause of the loss of plasti ity. The slow dynami s observed in this system presents analogies with the physi al aging e e ts observed and reported on glass studies [14, 15℄: as the system gets older (ages) there is a loss of plasti ity for stru tural or mole ular relaxation and less hanges are observed during the relaxation time. The me hanisms underlying the slow dynami s of glassy systems are the spatial (or geometri ) disorder related with the diÆ ulty to satisfy simultaneously all mi ros opi intera tions, a hara teristi alled frustration. Glasses and spin glasses exhibit a rough energy lands ape with many lo al minima whi h is responsible for slowing down the relaxation towards equilibrium. Their dynami s depends on the past history of the sample and under small perturbations they relax slowly. In our biologi ally motivated model the non-lo al intera tions are based on omplementarity (both perfe t and slightly defe tive mat hes) and re e t the a tivation and suppression me hanisms in the omplementary regions of the shape spa e. The analogue to frustration, in the network, is generated by the inability of the system to in orporate information and to satisfy the onstraints of the a tivation and suppression me hanisms. We show that the CA dynami s gives rise to a large family of periodi attra tors whi h are very robust, and ould be regarded as the analogues of the minima in the glassy systems. The lost of plasti ity and aging e e ts will therefore be related to the non-ergodi ity in the phase spa e. In this paper the dynami s of the immunologi al responses in the network model are investigated using the ommon tools adopted in the study of glassy systems [15℄. In parti ular we will fo us on the relaxation of twotime auto orrelations of the B- ell populations, the stru ture of the attra tors and e e ts of perturbations and immunizations. The paper is organized as follows: in se tion 2 we des ribe the model and the pro edures adopted in order to measure the orrelations among di erent on gurations of the system; in se tion 3 we present and dis uss the results obtained for small and large perturbations and in se tion 4 we present a summary, prospe tives for future work and some on lusions. 2 2 The Model The model under study here is a modi ed version [6, 7, 8℄ of the model proposed by Stau er and Weisbu h [16℄, whi h in turn was inspired in a previous model proposed by de Boer, Segel and Perelson [17℄ to des ribe the time evolution of the immune repertoire. It is a deterministi window ellular automaton model based on the shape spa e formalism [18℄, whi h des ribes the intera tions of B-lympho ytes and antibodies, and the main me hanism underlying these intera tions, whi h is pattern re ognition (lo k-key intera tion). The dynami s of the system is regulated by spe i intera tions involving omplementary mole ular re eptors of the di erent lones. The memory about the relevant antigens, presented to the system during its past history, emerges from the dynami s of the system, rather than being stored in a stati registry. To ea h point of a d-dimensional dis rete latti e we asso iate a given re eptor, whi h in turn will be des ribed by d- oordinates representing important physi al- hemi al features of the re eptor that di erentiate one from the other [18℄. Sin e lones are lassi ed a ording to their mole ular re eptor, to ea h point ~r of the dis rete shape spa e we asso iate a three-state automaton B (~r; t) that will des ribe the on entration of its population over the time: low (B (~r; t) = 0), intermediate (B (~r; t) = 1) and high (B (~r; t) = 2). The time evolution of the ellular automaton is based in a non-lo al rule: population B (~r; t) at site ~r is in uen ed by the populations at site ~r (its mirror image or omplementary shape) and its nearest-neighbors ( ~r + Æ~r) (representing defe tive lo k-key intera tions). The in uen e on the population at site ~r due to its omplementary populations is des ribed by the eld h(~r; t): h(~ r; t) X = ~ r0 2( 0 B (~ r ; t) (1) ~ r+Æ~ r) where for a given ~r the sum runs over the omplementary shape ~r 0 = ~r and its nearest neighbors. Due to the nite number of states of the population B , the maximum value of the eld h(~r) is hmax = 2(2d + 1). The updating rule is based on a window of a tivation, whi h des ribes the dose response fun tion involved in B- ell a tivation [6, 16, 17℄. There is a minimum eld ne essary to a tivate the proliferation of the re eptor populations (1 ), but for high doses of a tivation (greater than 2 ) the proliferation is suppressed. The updating rule may be summarized as: B (~ r; t + 1) =  B (~ r ; t) + 1 B (~ r ; t) 1 if 1  h(~ r; t)  2 (2) otherwise = 3. We de ne the densities of sites in state i at time t but no hange is made if it would lead to B = 1 or B as Bi (t) (i = 1, 2, 3). The initial on gurations are randomly generated a ording to the following on entrations: B1 (0) = B2 (0) = x=2, while the remaining Ld (1 x) sites are initiated with B (~ r; 0) = 0. This model may exhibit stable or haoti behaviors depending on the values of x, 1 and 2 1 . However it is on the transition region between the two behaviors that the model behaves like a multi- onne ted network [6℄. In order to simulate the immunization proto ol performed in the mi e experiments we have followed the pro edure whi h is des ribed in Ref. [9℄ and summarized below. We have adopted the s ale of 5 time steps orresponding to 1 day [9℄. While the system evolves a ording to the deterministi dynami s (eq. 2), small and large perturbations an be produ ed, by setting the state of the hosen sites at B (~r; t) = 2. Small Perturbations The small perturbations a ount for the immunologi al stimuli (noise) oming from the environment. The time interval between two onse utive small perturbations is a random number uniformly distributed between 1 and 100 time steps. Ea h perturbation orresponds to a random number of damages (from 1 to 3) introdu ed at regions of resting populations (B = 0) whi h are randomly drawn (at every perturbation). The size of ea h damage may also vary randomly from 1 to 3 (onion-like) on entri layers around a entral site ( ontaining 7, 25 or 63 populations, respe tively, in 3 dimensions). Large Perturbations The large perturbations orrespond to the immunization proto ol whi h starts at a predetermined age of the mi e. Like in the experiments, we stimulate the system periodi ally (every 35 steps ' 1 week), and always in the same region (whi h is initially hosen at random but kept un hanged along the simulation). The damage size in this ase orresponds to six layers (377 populations) around an spe i site. Previous results on this model have shown [9℄ that the response to the immunizations presents a strong dependen e on the initial time (\age") at whi h the periodi proto ol starts ( tting experimental data extremely well for mi e whose immunization proto ol started at di erent ages). This has led us to the onje ture that the dynami al behavior of the system might be at least qualitatively similar to that of some glassy systems, despite the non-Hamiltonian nature of the CA dynami s. One of the quantities ommonly used in the study of glassy systems is the two-time auto orrelation fun tion between the system on gurations at two given times t and tw . A ommon experiment in glassy systems onsists in preparing the system at a high temperature and suddenly making a quen h to a low temperature. Then the system is allowed to relax up to a waiting time tw , whose on guration is re orded. As the system ontinues to relax the auto orrelations between the instantaneous on gurations at time t > tw and that at time tw are omputed. The waiting time tw is alled the age of the system in ontext of glassy systems. The monitoring of the two-time auto orrelations gives important insights on the relaxation pro ess. The dynami s, whether stationary 3 or not, an be readily re ognized on the tw dependen e of the auto orrelations, sin e in a stationary pro ess two-time quantities depend only on time di eren es. Consequently by plotting orrelations as a fun tion of time di eren e t tw for di erent waiting times it is possible to distinguish between an essential out of equilibrium pro ess from a stationary one. The aging pro esses observed in glassy systems are then related to the la k of temporal invarian e [15℄. Inspired on this approa h, we will analyze the multi- onne ted network dynami s by de ning and analyzing quantities analogous to the two-time auto orrelations for the CA: Ctot (t; tw ) C22 (t; tw ) = = 1 N X N Æ (B (~ r; tw ); B (~ r; t)) ~ r X (3) N 1 B2 (tw )N ~ rjB (~ r ;t w )=2 Æ (B (~ r; t); 2) ; (4) where N = Ld is the total number of populations in the system. Ctot and C22 amount to normalized proximities (using Hamming distan e as a measure) and from now on will be referred to simply as orrelation fun tions. These quantities will be analyzed for di erent proto ols: without any perturbation, with small perturbations (noise) and with large perturbations (immunization proto ol) in order to di erentiate the e e ts of the di erent kinds of perturbations. A omplementary view of the long-term behavior of the system an be obtained by looking at the attra tors to whi h the system evolves. In the present ase it will onsist mainly of y les, whi h will be re e ted by the periodi ity of the orrelations. For this purpose, we have obtained the return maps of the onse utive maxima of the orrelation fun tions. Note that the period of the maxima thus obtained does not orrespond to the period of the real system, whi h is at least twi e as long [19℄. 3 3.1 Results Without perturbations In previous works [19, 6, 7℄ it has been shown that depending on the initial on entration x of a tive populations the system may exhibit periodi or haoti behavior. Systems with low initial on entrations of a tive sites (x < x ) evolve to either xed points or orbits with short periods, while for x > x haoti attra tors appear. However, the biologi ally relevant region is in the transition region between these two behaviors, where the system rea hes one of several periodi orbits (as will be seen below) with a very long period and after a long transient. From now on, all the results have been obtained using the same parameters adopted in Ref. [9℄: d = 3, 1 = hmax =3, 2 = 2hmax =3 and x = 0:26 (on the transition region). Without any perturbation, the system evolves after a transient time towards a y le with a large period, as shown in the return map of Fig. 1. We have also varied the waiting time from 10 to 100000 (not shown). The greater the waiting time, the loser to the attra tor the system is, whi h is revealed by larger values of the auto orrelation fun tions. On e tw is larger than the transient, the time series for Ctot (t; tw ) (and also C22 (t; tw ) will in lude unity (see Fig. 1) and will not hange for larger values of tw . A typi al result for the time evolution of the orrelation fun tions for tw = 100 is shown in Fig. 2. The upper panel orresponds to the time evolution of the densities B1 (intermediate on entrations) and B2 (high on entrations), while the lower panel orresponds to Ctot (t; tw ) (open ir les) and C22 (t; tw ) ( lled ir les). Noti e that while the on entrations relax to approximately onstant values in a short time, Ctot and C22 take mu h longer to rea h their attra tors (note the logarithmi time s ale). In order to study how the system relaxes towards the attra tor it is more onvenient to make use of C22 , sin e it measures the hanges in the (more relevant) network of a tivated populations. In the ase shown in Fig. 2 the transient time needed to attain the attra tor is O(104 ) steps. It is ne essary to point out, however, that this typi al relaxation time is important only for the physi al aspe ts of the dynami s. When mapped into the biologi al problem, it would orrespond to  5:5 years, whi h is mu h longer than the average life time of the mi e used in the kind of experiment we simulated. Therefore the relevant behavior, from the biologi al point of view, happens to be in the transient of the model and not in its stationary state, a result whi h is interesting on its own. 3.2 Random Small Perturbations How does the behavior of the system hange when random small perturbations are produ ed on the parameter region used to simulate the real experiments performed with mi e [9℄? In order to investigate this issue, only the small perturbations des ribed in Se tion 2 were produ ed, starting at time zero. In Fig. 3 we ompare the time evolution of the densities and orrelation fun tions obtained for single runs for the system without perturbations and with small random perturbations. Note that in the ase with small perturbations the orrelations initially follow those of the purely deterministi system, showing only small di eren es. After about 103 104 steps, however, they start de reasing faster, indi ating some sort of umulative e e t that drives the system away from the region in phase spa e that it had approa hed until tw . These e e ts are more easily 4 1.0005 1 Ctot(t+1) 0.9995 0.999 0.9985 0.998 0.9975 0.9975 0.998 Figure 1: Return map for maxima of the total 0.9985 0.999 Ctot(t) orrelation fun tion 0.9995 1 1.0005 Ctot (t) for L = 50 and tw = 10000 (single run). 0.03 B1 B2 B1, B2 0.025 0.02 0.015 0.01 0.005 0 1 10 100 1000 10000 100000 1e+06 10000 100000 1e+06 t 1 0.95 Ctot 0.9 0.85 0.8 0.75 Ctot c22 0.7 1 10 Figure 2: Densities of intermediate and high for 100 1000 t-tw on entrations vs. time (upper panel) and auto orrelations vs. L = 50 and tw = 100, without any perturbation. t tw noted for C22 . Moreover, the perturbations do not alter the behavior of the densities, as expe ted, sin e the number of a tivated populations is kept approximately onstant by the self-regulatory me hanisms embedded in the dynami al rules. The hanges are observed only in the populations that belong to the a tive network: in order to in orporate new information (new populations), older ones are dea tivated. From the biologi al point of view, the realizations of the perturbations (small and large) di er from one individual to another, building the identity of ea h individual. At the end of life ea h individual will have a di erent history in terms of perturbations (antigen presentations) translated into the on guration of the populations belonging to the a tive network. This is ni ely illustrated by Fig. 4, where two initially identi al opies of a system undergo di erent realizations of the small perturbations a ording to the proto ol des ribed in se tion 2. The Hamming distan e between them grows on a long time s ale, revealing the me hanisms behind the behavior of the orrelations in Fig. 3. Returning to Fig. 3, it is important to stress that the hanges observed in the on entrations for very long times (upper panel) are due to nite size e e ts. They are aused by the fa t that all small perturbations are produ ed on regions of resting populations. For a nite system, after a long time all the possibilities will have already been explored. In reasing the size of the system, the hanges on the densities disappear. This is shown in Fig. 5, where we repeat the simulations, under the same onditions of Fig. 3, for a larger system (L = 100). 5 0.03 B1 B2 p33 B1 p33 B2 0.025 B1, B2 0.02 0.015 0.01 0.005 0 Ctot 1 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 0.55 0.5 10 100 1000 t 10000 100000 1e+06 1000 t-tw 10000 100000 1e+06 Ctot c22 p33 ctot p33 c22 1 10 100 Figure 3: Densities vs. time (upper panel) and orrelations vs. t tw for L = 50 and tw = 100; without perturbation ( ir les) and with small perturbations (triangles). 0.1 0.09 0.08 0.07 0.06 HD 0.05 0.04 0.03 0.02 0.01 0 1 10 100 1000 10000 100000 1e+06 t Figure 4: Normalized Hamming distan e between two initially identi al on gurations subje ted to di erent sequen es of small perturbations (three samples). In order to test out the robustness of these results, we have varied the parameters ontrolling the proto ol of the random small perturbations. For instan e, by in reasing the maximum number of di erent perturbations from 3 to 6 at ea h presentation, and/or by hanging the maximum time interval between onse utive perturbations from 100 to 10, we observe the same qualitative behavior, with a faster de rease of the orrelations | as expe ted, sin e in both ases the noise has been in reased. What happens when the system's auto orrelations de rease due to the small random perturbations after it has rea hed a periodi attra tor? Is the system driven to another y le? To answer this question we have studied the stability of the y les using the following pro edure: we let the system evolve without perturbation towards its attra tor for tw = 10000 time steps, after whi h we perturb the system with random small perturbations during a time interval
t1 . Then we turn o the perturbations and allow the system to relax during another time interval
t2 after whi h we obtain the return map of the orrelations in the following 200 time steps. If we produ e only one perturbation at tw (
t1 = 1) the system remains in its original y le, yielding a return map similar to that of 6 0.25 B1 B2 p33 B1 p33 B2 B1, B2 0.2 0.15 0.1 0.05 0 1 10 100 1000 t 10000 100000 1e+06 10 100 1000 t-tw 10000 100000 1e+06 1 0.95 Ctot 0.9 0.85 0.8 Ctot c22 p33 ctot p33 c22 0.75 0.7 0.65 1 Figure 5: Densities vs. time (upper panel) and auto orrelations vs. t tw for L = 100 and tw = 100 (the same onditions used in Fig. 3); without perturbation ( ir les) and with small perturbations (triangles). 1 0.9995 Ctot(t+1) 0.999 0.9985 0.998 0.9975 0.997 0.997 0.9975 0.998 0.9985 Ctot(t) 0.999 0.9995 1 Figure 6: Return map for maxima of total orrelation fun tion, for L = 50, tw = 10000,
t1 = 1000 and
t2 = 10000. Fig. 1. Then, under the same onditions we repeat the simulations adopting now
t1 = 1000 and
t2 = 10000. We observe in Fig. 6 that the return map has hanged slightly when ompared to Fig. 1. In parti ular, it no longer shows Ctot = 1 in the time series, but remains periodi , signaling that the system has shifted to a di erent y le due to the perturbations. Apparently the periods observed in Figs 1 and 6 are the same. The distribution of periods and transient times are urrently under investigation, results will be published elsewhere. Fig. 6 remains the same by in reasing
t2 , whi h guarantees that the di eren es with respe t to Fig. 1 are not a transient e e t. A ording to these results, the role the noise plays, if allowed enough time to perturb the system signi antly, is to drive the system from one attra tor to a nearby one, whi h suggests that there is a family of periodi attra tors whi h an be very lose to one another in phase spa e. These e e ts, however, take pla e on a time s ale whi h is longer [O(103 104 ) time steps℄ than the lifetime of the mi e [O(102 103 ) time steps℄. The small perturbations are therefore of little importan e to the CA dynami s in the biologi ally relevant time s ale, as suggested in previous work [9℄. The behavior of the auto orrelation fun tions in Figs. 3 and 5 is reminis ent of what is observed in glassy systems. Here the CA approa hes a periodi attra tor, being thereafter de e ted by the small perturbations. In 7 1 0.95 C(t, tw) 0.9 0.85 0.8 0.75 0.7 0.65 1 Figure 7: 10000. 10 100 1000 t-tw 10000 100000 1e+06 C22 (t + tw ; tw ) for di erent waiting times before small perturbations. From bottom to top, tw = 100, 5000, glassy systems, temperature drives the system from one lo al minimum to another, preventing it from getting stu k in lo al minima of the potential energy surfa e. As the physi al age of the system ( hara terized by the value of the waiting time tw ) grows, it nds itself exploring deeper and deeper regions of a rough potential energy surfa e, di using towards equilibrium [20℄. Due to the roughness of the potential energy, equilibrium is only attained on very long time s ales and relaxation is very slow. The older the system is the more it gets on ned to a restri ted region of phase spa e and the time s ales for relaxation get longer and longer. Eventually, if we wait long enough, the system equilibrates and the dynami s be omes stationary, losing sensitivity to the waiting time. This pi ture of an aging physi al system is reminis ent to the loss of plasti ity for adaptation in a living organism as it gets older. In either ase, one observes a strong dependen e on the waiting time tw , evident when measuring two-time quantities like auto orrelation fun tions and responses. Results for the CA model are shown in Fig. 7, where we see the de ay of auto orrelations for three systems with di erent ages or waiting times. Note that the horizontal axis is the time di eren e between the total time and the waiting time. The three urves should ollapse in the ase that the dynami al evolution is stationary. 3.3 Large immunizations What is the role of the large perturbations on the dynami s of the system? From previous studies we would expe t the large perturbations to a elerate the aging pro ess: while the random small perturbations would hange the route to the natural attra tor of the system, the large ones would redu e the transient time, a onje ture that would explain the loss of plasti ity of the older mi e [9℄. The proto ol adopted is the one des ribed in Se tion 2, with six-layer perturbations every 35 time steps, always at the same sites. In Fig. 8 we ompare the results obtained for the unperturbed system and those of the system subje ted only to large immunizations starting at tw (note that the densities are now also plotted as fun tions of t tw ). Somewhat surprisingly, the de rease of the orrelations for large perturbations is small, when ompared to the ase of the small perturbations ( ompare with Fig. 3). In hindsight, however, this an be understood be ause the large perturbations are always produ ed in the same region . For very long times the orrelations will attain a stationary regime whose a tive populations ontain at least part of the region involved in the immunizations [9℄. In other words, the driving produ ed by the large periodi perturbations is of a ompletely di erent nature than that of the small perturbations. The large perturbations seem to play a sele tive role: the y les that the system an rea h are restri ted to those that ontain at least part of the populations in orporated during the immunization. A ording to this pi ture, the aging e e ts observed in Ref. [9℄ ould be simply related to the exploration of the phase spa e: the older the system the less possibilities of hoosing new y les it would have. In Se tion 4 we dis uss this issue further. In Fig. 9 we ompare the time evolution of the densities B1 and B2 and the auto orrelation fun tions for the system subje ted only to large perturbations and for the system with both small perturbations and large immunizations. The results, as expe ted, on rm the dominan e of the small perturbations, over the large ones, in driving the system faster to a di erent attra tor. The in rease of the densities around t  104 for L = 50 8 B1, B2 0.065 0.06 0.055 0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 B1 (p0 P16) B2 (p0 P16) B1 (p0 P00) B2 (p0 P00) Ctot, c22 1 10 100 1000 10000 100000 1000 10000 100000 t - tw 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 Ctot (p0 P16) c22 (p0 P16) ctot (p0 P00) c22 (p0 P00) 1 10 100 t-tw Figure 8: Densities vs. time (upper panel) and auto orrelations vs. L = 50 and tw = 100. t tw for a system with large immunizations, Triangles: without perturbations of any kind. Cir les: both large immunizations and small perturbations. 0.07 B1, B2 0.06 0.05 B1 (p33 P16) B2 (p33 P16) B1 (p0 P16) B2 (p0 P16) 0.04 0.03 0.02 0.01 ctot, c22 1 10 100 1000 10000 100000 1000 10000 100000 t - tw 1 0.95 0.9 0.85 0.8 0.75 0.7 0.65 0.6 ctot (p33 P16) c22 (p33 P16) ctot (p0 P16) c22 (p0 P16) 1 10 100 t-tw Figure 9: Densities of a tive populations vs. time (upper panel) and auto orrelations vs. L = 50. t tw for tw = 100 and Triangles: large immunizations only. Cir les: both large immunizations and small perturbations. orresponds to the same nite size e e ts o urring in Fig. 3, being asso iated to the a tive populations whi h have been in orporated into the network by means of the immunization proto ol. 4 Con luding Remarks We have analyzed the dynami s of a ellular automata model for the B- ell repertoire whi h has reprodu ed experimental results of immunizations on mi e. Our analyses provide a broader ontext in whi h memory and plasti ity take pla e, as dis ussed in Ref. [9℄. 9 55000 50000 18500 18400 18300 18200 18100 18000 17900 17800 17700 100 45000 HD 40000 35000 30000 25000 1000 10000 100000 20000 15000 1 10 100 1000 10000 100000 t Figure 10: Hamming distan e between 100 pairs of random on gurations as a fun tion of time (average and standard deviation). Inset: evolution for t > 100 shows that the HD attains a stationary regime for suÆ iently long times (standard deviation not shown). De ning quantities analogous to the auto orrelation fun tions used in glassy systems, we have shown that the biologi ally relevant phenomena takes pla e in the transient regime of the model. The dependen e of these fun tions on both t and tw (as opposed to the di eren e t tw only) is reminis ent of the \aging" behavior observed in glassy systems, despite the fa t that the underlying dynami s of both systems are ontrolled by ompletely di erent me hanisms. Starting from di erent initial onditions the deterministi dynami s takes the system towards a family of long-period attra tors. When subje ted to random small perturbations (Fig 4), the system is driven towards a new attra tor of the family, revealing that most of the noise is assimilated. However, only part of the large perturbations is in orporated, due to the me hanisms of a tivation and suppression, leading to a saturation of the learning pro ess [9℄. From the biologi al point of view, the history of the mouse (sample) will be written by the di erent antigen presentations (random perturbations), starting from its \birth" (initial ondition). Sin e the system is large but nite there is a maximum amount of information it may in orporate. The loser it is to its \destiny" attra tor, the less information it is able to learn, sin e the deeper it already is in a given basin of attra tion. Therefore the biologi al aging orroborated by the experimental results may be simply a onsequen e of this dynami al feature. From the results obtained up to now there are eviden es that the purely deterministi dynami s is therefore non-ergodi . Despite their large number, however, the family of periodi attra tors o upy only a fra tion of the phase spa e. Eviden es of this ompression in phase spa e has been obtained in the following omputer experiment: sele ting randomly two initial on gurations (with the same initial on entration), we measured the Hamming distan e between them as a fun tion of time as both systems evolve without any perturbation. In Fig. 10 we show the evolution of the average HD for 100 pairs (error bars are standard deviations). During the rst 10 time steps the HD de reases very rapidly, but for t > 100 we observe a quasi-stationary regime, re e ting the slow driving to the attra tors (see inset). Note that this behavior is somewhat similar to that of the auto orrelation fun tions for small tw . Still fo using on the eviden e of phase spa e ompression, in Fig. 11 we present the distribution of HD sampled from 500(500 1)=2 pairs, for two onse utive time steps (t = 4000 and t = 4001). The distributions an be well des ribed by a Gaussian, with a width that remains approximately onstant (even for long times). Noti e that the average value for t = 4001 is slightly larger than for t = 4000, re e ting the os illatory behavior of the average HD as the pair of samples rea hes their periodi attra tors. It should be noted that, for large N , the Central Limit Theorem assures that randomly hosen initial on gurations (with the same x) naturally give rise to a Gaussian distribution of the HD between any two of them, at time zero. Interestingly, the CA dynami s does not hange the shape of the distribution, its only e e t being to essentially shift the Gaussian towards lower mean values, on a long time s ale. The spatio-temporal stru ture of the y les, as well as their transients and basins of attra tion, would be the obje t of further study and will be published elsewhere. This work was partially supported by CNPq, Faperj, Finep, Capes and Projeto Enxoval-UFPE. MC and RMZS a knowledge IF-UFF, where part of this work was done during their stay there. 10 90 80 No. of occurrences 70 60 50 40 30 20 10 0 15000 16000 17000 18000 19000 20000 21000 22000 23000 Hamming distance 80 t+1 No. of occurrences 70 60 50 40 30 20 10 15000 16000 17000 18000 19000 20000 21000 22000 23000 Hamming distance Figure 11: Distribution of Hamming distan e sampled over 500(500-1)/2 on guration pairs at time t = 4000 and = 4001. Verti al lines show the average value. t Referen es [1℄ Janeway, C. A, Travers, M. W. and Capra, J., Imunobiology: The Immune System in Health and Disease, 4th ed., Elsevier S ien e Ltda. Garland Publishing, New York, 2000. [2℄ Vaz, N. M. and de Faria, A. M. C., Guia Completo de Imunobiologia, Coopmed Editora, Belo Horizonte, Brazil (1993). [3℄ Jerne, N.K., Ann. Immuno. (Inst. Pasteur) 125 C, 373-389 (1974). [4℄ Holmberg, D., Anderson,  A, Carlsson, L. & Forsgren, S., Immunol. Rev. 110, 84-103 (1989). [5℄ Coutinho, A., Immunol. Rev. 110, 63-87 (1989). [6℄ Zorzenon dos Santos, R. M. and Bernardes, A. T., Physi a A 219, 1-12 (1995). [7℄ Bernardes, A. T. and Zorzenon dos Santos, R. M., J. Theor. Biol. 186 173 (1997). [8℄ Zorzenon dos Santos, R. M., in Annual Reviews of Computational Physi s VI, 159-202, ed. by D. Stau er, World S ienti (1999). [9℄ Zorzenon dos Santos, R. M. and Bernardes, A. T., Phys. Rev. Lett. 81, 3034 (1998). [10℄ Verdolin, B., MS . Thesis, Departamento de Bioqumi a e Imunologia, Instituto de Ci^en ias Biologi as, UFMG, Belo Horizonte (1997); Faria, A.M.C. et al., Braz. J. Med. Biol. Res. 31, 35-48 (1998); Vaz, N. et al., S and. J. Immunol. 46, 225-229 (1997). [11℄ Faria, A. M. C. et al., Me h. Ageing Dev. 102, 67-80 (1998). [12℄ Lahmann, W. M., Menezes, J. S., Verdolin, B. A. and Vaz, N. M., Braz. J. Med. Biol. Res. 25 813-821 (1992). [13℄ Bernardes, A. T. and Zorzenon dos Santos, R. M., Int. J. Mod. Phys. C 12, 1 (2001). [14℄ Strui k, L. C. E., Physi al aging in amorphous polymers and other materials, (Elsevier, Houston, 1976). [15℄ Bou haud, J.-P., Cugliandolo, L. F., Kur han, J. and Mezard, M., in Spin Glasses and Random Fields, A.P. Young editor, (World S ienti , Singapore, 1998). 11 [16℄ [17℄ [18℄ [19℄ [20℄ Stau er, D. and Weisbu h, G., Physi a A 180, 42-52 (1992). De Boer, R. J., Segel, L. A. and Perelson, A. S., J. Theor. Biol. 155, 295-333 (1992). Perelson, A. S. and Oster, G. F., J. Theor. Biol. 81, 645-670 (1979). Zorzenon dos Santos, R. M., Physi a A196, 12 (1993). Debenedetti, P. G. and Stillinger, F. H., Nature 410, 259-267 (2001). 12