Is the immune network a complex network?
Some years ago a cellular automata model was proposed to describe the evolution of the immune repertoire of B cells and antibodies based on Jerne’s immune network theory and shape-space formalism. Interactions among different B-cell clones, which may be found in low, intermediate and high concentrations, occur depending on the complementarities of their characteristic proteins and are regulated by activation and suppression mechanisms. Depending on the region of the parameter space, the model exhibits either stable (ordered) or chaotic behavior, but it is in the transition region of the parameter space between both regimes where we obtain the complex behavior of a self-regulated network, much like Jerneâs idea of the immune network. In this region, the network maintains the memory of the large perturbations, which simulate antigen presentations and reproduce immunization and ageing experiments performed with mice. Here we investigate if the networks generated by this model in the different regimes can be classified as complex networks. We have found that in the chaotic regime the network has random characteristics with large, constant values of clustering coefficients, while in the ordered phase, the degree distribution of the network is exponential and the clustering coefficient exhibits power law behavior. In the transition region we observed a mixed behavior (random-like and exponential) of the degree distribution as opposed to the scale-free behavior reported for other biological networks. Randomness and low connectivity in the active sites allow for rapid changes in the connectivity distribution of the immune network in order to include and/or discard information and generate a dynamic memory. However it is the availability of the low concentration nodes to change rapidly without driving the system to pathological states that allow the generation of dynamic memory and consequently a reproduction of immune system behavior in mice. Although the overall behavior of degree correlation is positive, there is an interplay between assortative and disassortative mixing in the stable and transition regions regulated by a threshold value of the node degree, which achieves a maximum value on the transition region and becomes totally assortative in the chaotic regime.
The main task of the immune system is to protect the integrity and identity of the body against any harm. The main cells of the immune system are macrophages and lymphocytes, the former being mostly responsible for innate immune responses and lymphocytes, which are responsible for adaptive or cell-mediated immune responses.
All immune cells carry a large number of molecular receptors (proteins) on the surface and the immune system works based on pattern recognition. There are two main classes of lymphocytes: T and B cells. While T cells are involved in signaling and functional activities for the majority of the pathogens, B cells are mainly responsible for the production of antibodies, which in general, function as markers for the pathogens to be phagocyted and improve the efficiency of the immune response. The antibodies produced by any B-cell population are copies of its molecular receptors and according to the clonal selection theory Janeway () proposed by Burnet in 1959, the antigen (or its binding sites), by pattern recognition chooses the B-cell clones (population of B cell and antibodies) that will proliferate.
According to estimates lymphocytes carry the order of molecular receptors on its surface and the human immune system is able to express the order of different receptors during its lifetime. Such large numbers allow for the recognition of any antigen presented to the immune system and for the completeness of the immune repertoire. If the repertoire is complete we should expect the elements of the immune system to recognize and be recognized by other elements; the same mechanism of recognition should also work for both antibody-antigen and antibody-antibody reactions. In 1974, Jerne Jerne (), taking these ideas into account, suggested that when the antigen is presented to the organism it will activate a set of B-cell clones and the production of specific antibodies (complementary to the antigen binding sites) that in turn would activate the receptors of other clones, and so on. Due to the interplay of activation and suppression mechanisms, the reaction chain would be finite, preventing the percolation of information through the entire system. This kind of dynamics generates a multi-connected network of cell populations that regulates the immune response, and at any time reflects the dynamic memory of the system regarding its previous history in terms of antigen presentation. In other words, the immune response to different pathogens (virus, bacteria, etc) is regulated by dynamics involving the complementary molecular receptors of the different B-cell clones that create a dynamic memory, which includes new and already existing information concerning the previous immunization process. This memory allows for a fast response to new presentations of previously seen antigens.
Since its proposal, little evidence has been put forward to support the existence of the immune network theory Coutinho89 (); Holmberg-Forsgren (); Lundkvist-Holmberg (), but research has suggested that if the network exists, then only of the lymphocytes will be activated, while the rest of the clones will form a pool of immunocompetent lymphocytes, which are able to recognize any antigen.
In 1992, Stauffer and Weisbush physA18042 () proposed a cellular automata model to describe the immune network based on Jerne’s ideas, shape-space formalism and a previous model introduced by De Boer and Perelson DeBoer-Segel-Perelson (). This discrete model was later modified by Zorzenon dos Santos and Bernardes physa2191 () and hereafter will be referred to as the BSP model. Using shape-space formalism, the model allows for the simulation of the large immune repertoire and the complementary interactions between B-cell clones. Each B-cell population is associated to a three-state automatum representing low, intermediate and high concentrations, and the interactions amongst the different populations are described by activation and suppression mechanisms. The dynamics of the model leads to stable (ordered) or chaotic behaviors that correspond to pathological rather than regular behaviors of the immune system. However, in the transition region between the two regimes physA-93 () during very long periods, the model describes aggregation-disaggregation dynamics with clusters splitting and fusing throughout time, as a multi-connected network of populations jtb186173 (). This network exhibits self-regulation and despite changes within the active populations, only 10-20% of the populations remain active. This multi-connected network subjected to multiple perturbations was used with mice to simulate and reproduce the behavior of their immune systems under multiple antigen presentations and the effects of ageing on generating these responses prl813034 ().
The dynamic behavior of the BSP model on the transition region was further investigated regarding different aspects. It has especially been observed that the memory of the system is dynamically allocated in order to incorporate new information regarding the history of the system without losing previously acquired important information Zorzenon-Copelli03 (). As the system becomes older, so less information is incorporated, which is a behavior associated to the loss of plasticity, as observed in other physical systems epjb34119 ().
Although there have been many studies on the parameter space and dynamics of the BSP model, it has not yet been characterized from the viewpoint of a complex network. Since it reproduces the behavior of real immune systems it is interesting to investigate its topological properties using graph theory formalism (see rmp7447 () and references therein) adopted in the characterization of complex networks found in nature nature41141 (); rmp7447 (). From previous studies on the dynamic behavior of the immune network it becomes clear that the regulation of the dynamics leading to stable, complex or chaotic behaviors emerges from the interplay between subnetworks of populations with low, intermediate and high activations. Our aim in this work is to investigate whether the networks generated by the BSP model in the different regimes (stable, transitional and chaotic) may be classified as complex networks. Therefore, we have focused on studying and characterizing the properties of these subnetworks for each phase with regard to degree distribution, the behavior of clustering coefficients, the existence of hierarchical structures and the correlations within the neighborhood (assortative or disassortative mixing behaviors) of the sites.
In what follows, in section II, we introduce the SW cellular automata model proposed to describe the immune network, in section III we undertake a short review of the main concepts involved in the study of the above-mentioned properties, in section IV we present and discuss the results obtained and in section V we make our final remarks.
Ii The cellular automata model
In the model introduced by Stauffer and Weisbuch physA18042 () shape-space formalism is adopted to describe the repertoire of B-cell clones and the interaction between different clones jtb81645 (). Each clone is represented by its molecular receptor, which corresponds to a point in a d-dimensional space. Thus, each receptor is characterized by properties corresponding to different characteristics, as for instance, the number of nucleotides, charge, hydrophobicity, etc. The nearest clone neighbors differ only by one property and the lock-key interactions occurring among populations with complementary molecular receptors is described by a non-local rule, where an on-site clone is influenced by those located at and its nearest neighbors (representing slightly defective interactions). According to estimates, if the shape-space notion is relevant in a continuous mathematical approach, we obtain jtb81645 () but if it is relevant in a discrete approach, then it is jtb186173 (). The population associated to each molecular receptor is represented by a three-state automaton describing its concentration at any given time: low (), intermediate () and high (). The influence on the population of site caused by its complementary populations is described by the field :
where for each the sum runs over the complementary shape and its nearest neighbors. Due to the finite number of B-cell population states, the maximum value of the field is . The rules describing the changes of B-cell populations due to their interactions with other populations are based on an activation window, which is inspired by a log-bell-shaped proliferation function associated to the receptor cross-linking involved in B-cell activation physA18042 (); physa2191 (); jtb186173 (). In this window there is a minimum field necessary to activate the proliferation of the receptor populations (), but there is also an upper limit for activation since for high doses of activation (greater than ) the proliferation is suppressed. The updating rules are summarized as follow:
but no change is made if it would lead to or . At each time step we define the densities of sites in state as , where = 0, 1 and 2 correspond to low, intermediate and high concentrations, respectively.
The initial configurations are randomly generated depending on the parameter that controls the initial concentrations: = =, while the remaining sites are initiated with low concentrations.
This model exhibits stable and chaotic regimes for physa2191 (), depending on the activation threshold () and the width of the activation window (). These different dynamic behaviors are separated by a transition region where the model exhibits complex behavior during cycles of very long periods physA-93 () in which clusters of different B-cell populations split and fuse jtb186173 () and the activated populations behave like a multi-connected network. In 1998 prl813034 () this multi-connected network was used to reproduce the results of experiments performed with mice and showed refractory behavior under multiple antigen presentations. Antigen presentation is simulated by the introduction of perturbations flipping the state of low concentration populations to those of high concentration, and computing the changes on the multi-connected network with respect to the initial state before the introduction of the perturbation. The same study has shown that the aging effects observed in the immune responses of mice can be explained by the loss of plasticity and the ability to bring about new changes in order to include new information in the system, i.e., the older the system, the more rapidly it saturates and the less intense is its response. It has also been found Bernardes-Zorzenon01 () that there is a characteristic cluster size associated to the loss of plasticity and a power law distribution for the permanence times (the time interval that each population remains activated or belongs to the multi-connected network). As expected, the absence of scale indicates that there is no typical permanence time, a feature that may be understood since the memory of the system is allocated dynamically at each time step depending on the interactions among the populations. Since the aging effects observed in the multi-connected network have similarities with physical glassy systems, auto-correlation functions were obtained epjb34119 (). As previously mentioned, the system with no perturbation is driven to a long-limit cycle-attractor after a long transient time. When subjected to small, random perturbations however, the very notion of a transient becomes fuzzy and the results in Ref. epjb34119 () show that the system is deflected from its attractor by small perturbations after time steps. Therefore, small perturbations will cause the system to change attractors from time to time due to their cumulative effects, and is reflected in the decreasing auto-correlation functions. Large perturbations would not accelerate the de-correlation process. In fact, large perturbations lead to a much weaker (slower) de-correlation, since they are always produced on the same sites in shape-space to simulate antigen presentations. Small perturbations can be more easily absorbed by the system than larger ones since they involve only local changes. The various studies performed until now jtb186173 (); epjb34119 (); Zorzenon-Copelli03 () to understand the dynamics of the multi-connected network have suggested that different behaviors (stable, complex and chaotic) emerge from the interplay between the subnetworks of low, intermediate and high concentrations due to the activation of suppression mechanisms. In this work, we have investigated the properties of these inter-dependent sub networks from the viewpoint of complex networks using graph theory formalism.
Iii The complex network properties investigated
In what follows we will consider the multi-connected network and the subnetworks as graphs composed of nodes and links. The key quantity on the study of the topological properties of graphs and complex networks is the degree distribution of the nodes. The connectivity or degree associated to the node is defined by the number of links that connects node to other nodes on the graph. In other words,
where the sum runs over all nodes not equal to and corresponds to connectivity or adjacent matrix elements. If there is a link between and , , otherwise . The degree distribution obtained for the network indicate the topological class of universality (regular, random or scale-free) to which the complex network belongs and whether there is an existence hidden variables that may deviate its behavior from one of the main universality classes.
Due to the complementary interaction that regulates the interactions between the B-cell populations, another quantity of interest in the study of the multi-connected network is the clustering coefficient of a given node, which is a measure of the degree to which nodes in a graph tend to cluster together. The clustering coefficient of a node () is defined as:
where is the number of connections between the neighbors of site . The average of all nodes gives the network clustering coefficient . While is a local property, is global. The global clustering coefficient gives an overall indication of the clustering in the network, whereas the local coefficient describes the embedding of single nodes. The behavior of C(K) as a function of K may indicate the presence of hierarchical structures in the network when it exhibits power-law behavior pre67026112 (). Hierarchical structures go beyond simple clustering by including the simultaneous organization of all scales in a network. Such structures are represented by trees or dendrograms, in which closely related pairs of vertices (nearest neighbors) have the lowest common ancestors than more distantly related pairs. In the case of the immune network we should expect to find some hierarchical structures, since it is built through complementary interactions using shape-space formalism, where similar shapes should have common ancestors.
To complete the characterization of the subnetworks we investigate the existence of assortative or disassortative mixing on the different subnetworks, a measurement that will define if the highly-connected nodes in each subnetwork tend to have more connections with other highly-connected nodes or with nodes with a lower number of connections, respectively. In order to investigate the assortativity of the network we calculate the degree-correlation of the complex network () knn_book (). When there is a tendency of the nodes to connect to other nodes with a similar degree, the correlation is assortative, while when there is a prevalence of links between nodes with dissimilar degrees, the correlations are disassortative. The average degree between nearest neighbors knn_book () is defined as:
The sum runs over all nodes with degree , is the number of nodes of degree and is the average nearest neighborâs degree of vertex :
where the sum runs over all nearest neighbors of node . When is the constant, the degrees of neighboring nodes are uncorrelated.
Another quantity of interest is the mean shortest path () that corresponds to the average number of connections between any two nodes and calculated over all pairs of nodes. For regular hypercubic lattices in d-dimension and in the case of random graphs grows logarithmically with the number N of nodes (. The small world effect corresponds to the case in which any pair of nodes is connected by the shortest path distance knn_book ().
By investigating these properties it was possible to better understand the dynamics of the immune network from a topological viewpoint and the interdependence of the three subnetworks.
All the results reported in this section were obtained using the following parameters: , , sites, and , where and they are representative of the behaviors obtained on the different regions analyzed. Thus, when varying the parameters we should expect variations in the localization and size of these regions in the parameter space, but the overall behavior of the quantities analyzed here in the different regions would be the same. For the sake of clarity we will refer to the subpopulations or subnetworks of different states of activation as: for low, for intermediate and for high. By subnetwork we mean the network formed by all nodes in the same state. In the present study, for the maximum number of neighbors of a given node is 24 and by neighbors of node we considered its nearest and next-nearest neighbors. Despite the fact that the dynamics of the model is based on complementary interactions and involve mirror images, the neighbors of a given site in a subnetwork correspond to the populations in the same state belonging to its neighborhood, as defined above. The average values used for the majority of graphs as well as the distributions, over 1000 samples were obtained, discarding the 1000 initial time steps, in order to guarantee that the samples correspond to thermalized states on the different regimes.
The transition region between stable to chaotic behaviors for this set of parameters is located in the vicinities of : the stable phase corresponded to the region of , and the chaotic associated to the region of . Figure 1 shows that the stable phase is dominated by populations in low concentration (), while in the chaotic phase the density of the three types of populations are of the same order of magnitude and therefore, of the nodes are activated ( and ) corresponding to a pathological state as previously suggested physA-93 (). All values shown in Figure 1 correspond to the values of the densities calculated after thermalization.
The degree distributions of the three subnetworks (, and ) in the stable phase are finite and exponential distributions, as shown in Figure 2. The positive exponent for distribution indicates a high probability that any given node belonging to this subnetwork will be linked to almost all its neighbors, thus generating large domains of low-activity nodes on the 3-D lattice. The negative slope obtained for and distributions indicate that in both subnetworks the active nodes have a high probability of being linked to very few neighbors. The stable regime would correspond to atypical states of the immune system with very few active populations, a feature that makes it difficult to maintain a memory regarding the perturbation history of the system, and thus explains why the mice experiments prl813034 () could not be reproduced in this region of the parameter space.
As shown in Figure 2, in the chaotic regime the degree distributions are similar, finite and Poissonian (or bell-shaped) with an average number of connections equal to , and for , and subnetworks, respectively. Therefore, theses subnetworks are random networks ErdosRenyi () as will be the entire network; this regime, as mentioned above physA-93 () would also correspond to a pathological state of the immune system, in which the order of of the populations are randomly activated corresponding to an overactive state equivalent to that of septicemia.
In the transition region () we find a completely different behavior of the distribution with respect to and distributions, as shown in Figure 3. The degree distribution indicates mixed behavior, exhibiting characteristics of both regimes, stable and chaotic. For the distribution is Poissonian and for it is exponential; thus, in the transition region the nodes have a large probability of being linked to more than 50% of its neighbors, forming large clusters. However, embedded in this structure there is a random network of nodes with low connectivity () exhibiting a characteristic connectivity or an average value of . As far as we know, this is the first time that such a type of mixed behavior of the degree distribution of a network has been reported on the literature. The degree distributions for and shown in Figure 3 although showing log-normal behavior can be identified as Poisson-like distributions with long tails, corresponding to random subnetworks with few highly connected nodes (long tails). These random subnetworks have active nodes with low connectivity, and therefore, under perturbation would be able to change more easily than highly connected nodes to incorporate information into the subnetwork. The average degree for these distributions ( and for and distributions respectively) is smaller than the average K obtained in the distributions of the chaotic regime, which in the case of is almost double.
As reported in the literature, many biological networks are scale-free and this structure would guarantee the identity and robustness of the network rmp7447 (); nature41141 (); LNP2004 (). The same behavior was obtained by Burns and Ruskin PA365549 () for the degree distribution of a different immune network model also based on shape-space formalism that was proposed to describe the development of the repertoire of immune cells. The BSP model is a different model describing only B-cell population interactions based on non-local rules inspired by the functioning of the immune system. The networks it generates under certain conditions are able to reproduce the dynamical memory and behavior observed in real immune systems. Poissonian behavior predominates in the degree distribution for the subnetworks in transition and the chaotic phases. These results are very interesting since they suggest that random structures are always present in the active subnetworks of these regimes and would be responsible for the dynamic allocation of memory observed in the transition region. The randomly attached active sites of low connectivity allow for rapid changes to add or discard information. The plasticity, which permits the memory of previous antigen (perturbation) history to be maintained, only observed on the transition region, will also depend on the behavior of the degree distribution of low concentration subnetworks exhibiting a mixed behavior (exponential and pure Poissonian) 3. In order to achieve the plasticity necessary to maintain the history of the perturbation (antigens) memory, it is necessary for the low connectivity subnetwork to have a random structure embedded on one that is very connected and almost regular: low connectivity sites would allow for rapid changes following the active populations while the high connected ones would keep the robustness and the integrity of the system. In the chaotic regime the fact that all distributions are Poissonian with a relatively low average connectivity allows for rapid changes in all subnetworks, and since the densities of different populations are of the same order of magnitude rapid changes would very easily lead to the over-activation of the populations. Therefore, the randomness on the degree distribution is an important feature for the immune network Jerne () to maintain a memory of the previous history of antigen presentations. The exponential behavior observed on the degree distribution for all subnetworks in the stable regime do not allow the memory to be maintained for two reasons: overall, the high density of low concentration sites is highly connected, which is an unfavorable aspect regarding the necessary changes to incorporate information, and the very low density of active sites does not allow for the accumulation of information.
We have also studied the behavior of the average degree of the neighbors of the nodes with degree . Overall, the results indicate that there is a positive correlation between the neighbors or assortative mixing properties, i.e., nodes with the same type of connectivity (low and high) tend to cluster together. Moreover, results indicate that there is an overall linear growth of as a function of K, but behavior is different for all three phases. In the stable regime grows almost linearly (positive correlation) for , and for , and distributions, respectively, but decreasing very rapidly (negative correlation) afterwards in the last two distributions (Fig. 4). In other words the positive correlation grows up to a critical value () of connectivity and becomes negative for . At the transition we observed that the average degree for distribution behaves linearly for all K, although with different slopes in at least two regions; the behavior of and persists as in the ordered phase, but in this case they grow linearly for and respectively, with a subsequent decrease, as in the previous regime (Fig. 4). The existence of cut-off degrees is characteristic of finite random networks, however the linear growth has not previously been observed. In the chaotic region () grows linearly with ((see Figure 4) for all subnetworks, as expected, since in this region, the degree distribution is Poissonian and the densities of nodes for all subnetworks are of the same magnitude. The positive correlation would allow for the over-activation of nodes as observed.
For the same set of parameters , in the transition and chaotic regions ( and respectively) we have estimated the mean shortest path for the three subnetworks. Due to the computational cost it is difficult to perform this calculation with the necessary statistics, which is why here we present the results obtained for a simple sampling. At the transition region we have obtained for , for and for and for the chaotic regime () the results are: for , for and for . In all cases we observe that (number of nodes in the subnetwork) but overall , i.e. the behavior between regular lattice and random networks. Therefore, is not necessarily a signature of small-world properties as suggested by other biological networks nature393440 (), but in our case it is a signature of the existence of very connected nodes. We should expect small-word properties as in other random networks, however due to the computational cost we have not looked for evidence of such properties.
When we compute the clustering coefficient as a function of for the three subnetworks (shown in Figure 5), we observe that for the subnetwork is approximately constant, with a slight increase from to in the transition region. This result indicates a strong aggregation between the populations at a low concentration, i.e., on average around percent of the neighbors are neighbors among themselves. However, the behavior of the clustering coefficients for subnetworks and differ completely from that of : for small in both cases the clustering coefficient is close to zero and start to increase around , reaching its maximum value () for . The increase of the clustering coefficient implies an increase in the aggregation of these subnetworks, reflecting an increase in the concentrations of active sites with intermediate and high concentrations as increases (Figure 1). We have also observed that the cohesiveness between different populations is always less significant than the cohesiveness among the elements of the same population (results not shown).
While the average clustering coefficient is almost constant as a function of in the chaotic and transition regions, it exhibits a power-law behavior in the stable region (), as shown in Figure 6. This behavior indicates that the subnetworks in the stable regime have hierarchical structures pre67026112 () with exponents that are smaller than unity ( for and for ) as observed in the case of the internet VPSV (). Hierarchical structures go beyond simple clustering by including organization at all scales in the network simultaneously. Trees, in which the nodes are connected, generally possess common ancestors that are topologically closer to them than nodes, which are not closely related. Since the numbers of populations with intermediate or high concentration are small in the stable regime, this would explain why under perturbation the system changes, incorporating information regarding perturbation and losing information on the previous perturbation, do not maintain any memory regarding its history epjb34119 (); Zorzenon-Copelli03 ().
Figure 6 shows the constant behavior of the average clustering coefficient as a function of K in the chaotic region. Constant behavior is characteristic of random networks where the clustering coefficient is where is the probability of connection between two nodes. For random networks the small aggregation is due to small ErdosRenyi (). Previous results indicate that in the chaotic region the subnetworks are random networks.
In this work, from the viewpoint of a complex network we have analyzed an immune network model, which was able to reproduce the behavior of the immune systems of mice. Indeed, the immune network is complex, and contrary to what has been reported in the literature regarding other biological networks, this biological network is characterized by random behavior rather than scale-free rmp7447 ().
In fact, the networks obtained in the three dynamical regimes observed in the model are composed of three different subnetworks of low, intermediate and high concentrations. In this work we have investigated the behavior of the three subnetworks in the different regimes. We observed the existence of random characteristics in the three regions: in the ordered phase the degree distributions are purely exponential; in the transition region while mix exponential () and Poisson (), the and are Possonian; and in the chaotic region the distributions are purely Poissonian. In reality, the immune network seems regulated by random subnetworks of active populations ( and ). The randomness of these structures formed by nodes with low average connectivity (with respect to the maximum number of neighbors) allow for rapid changes that add or discard network information. However the necessary plasticity that would allow for adaptation and maintainance of the networkâs identity and memory regarding its history in terms of perturbations (or antigen presentations) depends on the structure of the low concentration population subnetwork. The structure that allows the reproduction of immune system behavior in mice is composed of random substructures embedded in one that is very connected: low concentration sites with random low connectivity will allow for the incorporation of new information (new active populations in the network) leading to a different configuration that when compared to the previous one has added relevant new information and discarded the unnecessary.
Another important result that differs from those reported in the literature, is the fact that grows linearly as a function of K. For stable and transition regions the active population subnetworks show a linear growth (positive correlation) with a cut-off , as expected for finite random networks, followed by a negative correlation. However, distributions do not have a cut-off; in the chaotic region, all distributions exhibit a linear growth of as a function of K. Put another way, while the subnetworks exhibit assortative mixing in all phases, there is an interplay between assortative and disassortative mixing in and subnetworks at stable and transition regions, depending on the value of . If there is a positive, linear degree-correlation, but for the correlation is negative. There is a maximum threshold value in the transition region but that disappears in the chaotic region giving way to pure assortative behavior. This suggests that although random, the networks obtained with this model are much more complex than those thus far reported in the literature. We believe that the differences between our results and those from the literature arise from the fact that in this study the networkÂ´s topology reflects the dynamics inspired in the biological process attributed to the immune system, while the majority of biological networks that have been studied neither include nor reflect the dynamics of the process, since in many cases it is unknown.
We thank Silvio Ferreira Jr. for his helpful discussions and suggestions. This study received financial support from the following Brazilian institutions: the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, http://www.cnpq.br), the Coordenação de Aperfeioamento de Pessoal de Nivel Superior (CAPES, http://www.capes.gov.br) and the Fundação de Amparo Ciência e Tecnologia do Estado de Pernambuco (FACEPE, http://www.facepe.br-Pronex-EDT 0012-05.03/04 and Pronex-APQ 0203-1.05/08).
- (1) C. A. Janeway, P. A. Traver, M. Walport, and J. Capra, Immunobiology: The Immune System In Health And Disease (Garland Science Publishing, NY, 2009)
- (2) N. K. Jerne, Ann. Immuno. (Inst. Pasteur) 125 C, 372 (1974)
- (3) A. Coutinho, Immunol. Rev. 110, 63 (1989)
- (4) D. Holmberg, Anderson, L. Carlsson, and S. Forsgren, Immunol. Rev. 110, 89 (1989)
- (5) I. Lundkvist, A. Coutinho, F. Varela, and D. Holmberg, Proc. Natl. Acad. Sci. USA 86, 5074 (1989)
- (6) D. Stauffer and G. Weisbuch, Physica A 180, 42 (1992)
- (7) R. J. D. Boer, L. A. Segel, and A. S. Perelson, J. Theor. Biol. 155, 295 (1992)
- (8) R. M. Z. dos Santos and A. T. Bernades, Physica A 219, 1 (1995)
- (9) R. M. Z. dos Santos, Physica A 196, 12 (1993)
- (10) A. T. Bernardes and R. M. Z. dos Santos, J. Theor. Biol. 186, 173 (1997)
- (11) R. M. Z. dos Santos and A. T. Bernades, Phys. Rev. Let. 81, 3034 (1998)
- (12) R. M. Z. dos Santos and M. Copelli, Brazilian Journal of Physics 33, 628 (2003)
- (13) M. Copelli, R. M. Z. dos Santos, and D. A. Stariolo, Eur. Phys. J. B 34, 119 (2003)
- (14) R. Albert and A.-L. Barabási, Rev. Mod. Phys. 74, 47 (2002)
- (15) H. Jeong, S. P. Mason, Z. N. Oltvai, and A. L. Barabási, Nature 411, 41 (2001)
- (16) A. S. Perelson and G. F. Oster, J. Theor. Biol. 81, 645 (1979)
- (17) A. T. Bernardes and R. M. Z. dos Santos, Int. J. Mod. Phys. C 12, 1 (2001)
- (18) E. Ravasz and A.-L. Barabási, Phys. Rev. E 67 (2003)
- (19) A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge Univ. Press, 2009)
- (20) P. Erdős and A. Rényi, Publ. Math. 6 (1959)
- (21) B. A-L, Z. N. Oltvai, and S. Wuchty, Lect. Notes. Phys. 650, 443 (2004)
- (22) H. J. Ruskin and J. Burns, Physica A 365, 549 (2006)
- (23) D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998)
- (24) A. Vázquez, R. Pastor-Satorras, and A. Vespignani, Phys. Rev. E 65, 066130 (2002)