Spreading of Persistent Infections in Heterogeneous Populations
Abstract
Up to now, the effects of having heterogeneous networks of contacts have been studied mostly for diseases which are not persistent in time, i.e., for diseases where the infectious period can be considered very small compared to the lifetime of an individual. Moreover, all these previous results have been obtained for closed populations, where the number of individuals does not change during the whole duration of the epidemics. Here, we go one step further and analyze, both analytically and numerically, a radically different kind of diseases: those that are persistent and can last for an individual’s lifetime. To be more specific, we particularize to the case of Tuberculosis’ (TB) infection dynamics, where the infection remains latent for a period of time before showing up and spreading to other individuals. We introduce an epidemiological model for TBlike persistent infections taking into account the heterogeneity inherent to the population structure. This sort of dynamics introduces new analytical and numerical challenges that we are able to sort out. Our results show that also for persistent diseases the epidemic threshold depends on the ratio of the first two moments of the degree distribution so that it goes to zero in a class of scalefree networks when the system approaches the thermodynamic limit.
pacs:
87.23.Ge, 89.20.a, 89.75.FbI Introduction
Disease spreading has been the subject of intense research since long time ago MayBook (); Cambridge (); Murray (). Our current knowledge comprises mathematical models that have allowed to better understand how an epidemic spreads and to design more efficient immunization and vaccination policies MayBook (); Cambridge (); Murray (). These models have gained in complexity in recent years capitalizing on data collections which have provided information on the local and global patterns of relationships in the population geisel (); guimera (); colizza (). In particular, with the advent of modern computational resources and tracking systems, it is now feasible to contacttrace the way the epidemic spreads or at least to predict the paths that a given pathogen might follow. In this way, some of the assumptions at the basis of the theoretical models that were difficult to test the backbone through which the diseases are transmitted are now more accurately incorporated into epidemiological models egamstw04 (); cpv07 (); glmp08 (); cv08 (); mam09 ().
Strikingly, the systems on top of which diseases spread show common nontrivial topological and statistical properties siam (); PhysRep (). A large number of networks of contacts in realworld social, biological and technological systems have been found to be best described by the socalled scalefree (SF) networks. In SF networks, the number of contacts or connections of a node with other nodes in the system, the degree (or connectivity) , follows a power law distribution, . Recent studies have shown that the SF topology has a great impact on the dynamics and function of the system under study cnsw00 (); ceah01 (); siam (); PhysRep (). The reason is that, at variance with homogeneous or regular networks, SF architectures are a limiting case of heterogeneity where the connectivity fluctuations diverges if as the system size tends to infinity (the thermodynamic limit). This means that there are nodes in the network which has an eventually unbounded number of connections compared to the average degree. Examples of such networks include the Internet falou99 (); calda00 (), the worldwideweb (WWW) www99 (), foodwebs, and metabolic or protein networks strog01 (); PhysRep ().
In the context of disease spreading, SF networks lead to a vanishing epidemic threshold in the limit of infinite population when Vespignani (); MaySci (); mpv02 (); newman02 (). This is because the ratio determines the epidemic threshold above which the outbreak occurs. When , is finite while goes to infinity, that is, the transmission probability required for the infection to spread goes to zero. Moreover, the previous result holds both for the SusceptibleInfectedSusceptible (SIS) and SusceptibleInfectedRemoved (SIR) epidemiological models Vespignani (); MaySci (); mpv02 ().
In this paper, we will deal with a different kind of diseases those that are persistent in time and shows a latent period that can be as large as an individual’s lifetime. Our first aim is to enlarge the epidemiological framework for complex networks reported previously for the SIS model Vespignani () and proposed as well for the SIR model MaySci (); mpv02 () by integrating the spreading dynamics of persistent diseases within it. With this purpose, we consider a variation of the SusceptibleExposedInfectedRemoved model MB () on complex heterogeneous networks. As we will see, this kind of dynamics introduces new challenges essentially different to other successfully treated before. In particular, as the latent period is high enough, we shall work with an open system in which new individuals are born and others dead for causes not directly related to the spreading of the disease. We present novel analytical and numerical methods that allow us to obtain the epidemic threshold for this kind of disease in heterogeneous populations. Moreover, we introduce a numerical method that is wellsuited to deal with the kind of problems we face. Our results point out that also for persistent infections the virtually unbounded fluctuations of the degree distribution have an important enhancing effect on the epidemic spreading.
Ii The model
To be more precise and without loss of generality, we will particularize our model to one case of persistent infection the most threatening one which is tuberculosis (TB). TB is an old disease whose worldwide prevalence had been diminishing even before vaccination and prophylaxis strategies were firstly accomplished Wilson (); Styblo (); Daniel (). Its recent return in developing countries, mainly in Southeast Asia, have attracted renewed interest in it. The current world estimate of prevalence is about 33% while the number of deaths per year that it is causing reaches more than 3 million people Bleed (). Depending of the kind and the intensity of immune response that the host immune system performs after initial infection with M. tuberculosis bacillus, the individual can suffer latent infection, (in which the bacteria are under a growtharrest regime and the individual neither suffer any clinical symptom nor becomes infectious) or active infection, where the host suffers clinical symptoms and can transmit the pathogen by air Comstock (); Styblodos (). Latently infected individuals can, generally after an immunedepression episode, reach the active phase. Estimating the probability of developing direct active infection after a contact, or alternatively, the lifetime’s risk for a latent infected individual to evolve into the active phase, are not easy tasks. However, it is generally accepted that only 510% of the infections directly produce active TB Comstock (); Styblodos (), while the ranges concerning the estimation of typical “halflife” of latent state rounds about 500 years MB ().
The spreading dynamics of TB like diseases has been studied in recent years. However, to the best of our knowledge, these works assume the homogeneous mixing hypothesis, that is, a perfectly homogeneous system in which all individuals are dynamically equivalent. As mentioned above, many of the systems on top of which diseases spread, are better described by scalefree connectivity patterns. Therefore, in what follows, our main objective will be to assemble a basic model fitted for tuberculosis spreading that firstly takes into account the heterogeneity in the distribution of the networks of contacts. We note that the increasingly alarming situation about TB epidemiology evidences the need to increase the effort in TB research in a global way. In the context of the study of its epidemiology, new models must be developed in order to gain predictive skills, incorporating the recent theoretical advances referring to disease spreading on complex heterogeneous substrates as well as metapopulation approaches and new computational tools for numerical analysis and simulation. In this sense, ours is a first contribution that addresses one of the most important parameters in epidemiological description: the epidemic threshold.
Let us then introduce our model. We consider that individuals in the population are compartmentalized into three groups: healthy , infected but not infectious or latently infected  and sick individuals which are infected and are infectious as well. The transition between these subpopulations proceeds in such a way that a healthy individual acquires the bacteria through a contact with an infectious subject with probability . In its turn, this newly infected individual may develop the disease directly with probability . However, the most common case is the establishment of a dynamic equilibrium between the bacillus and the host’s immune system, which allows the survival of the former inside the latter. When this happens, newly infected individuals become latently infected, because despite harboring the bacteria in blood, neither becomes sick nor is able to infect others.
On the other hand, after a certain period of time (which may be several years) and usually following an episode of immunosuppression, the balance between the bacterium and its host can be broken. In this case, the bacteria grow and the individual falls ill beginning to develop the clinical symptoms of the disease. In addition, if the infection attacks the lungs (pulmonary TB), the bacillus is present in the sputum, making the guest infectious.
The dynamics of the disease, in a well mixed population, is then described by the following system of nonlinear differential equations:
(1) 
in which:

represents the total population at time ,

is the number of contacts per time unit,

is the probability that the bacteria is transmitted to a new host after a contact with an infectious subject,

is the birth rate per capita and per unit time,

is the natural death rate per capita and unit time,

is the rate of diseaserelated deaths per capita and unit time,

is the transition frequency of latent infection (i.e., the probability that a latently infected individual becomes infectious),
with the closure relationship:
(2) 
The model above is a variation of the archetypal SIR model to which a fourth class has been added (latency class L). This kind of model has been largely treated in the literature in its wellmixed version, and it is frequently referred as SEIR model. As a first step, in this work, we will identify the removed individuals mostly with dead ones, and therefore we do not consider the possibility of natural or medical recovery (this simplification is in part justified by the large latency period of infected individuals and the constant flow of newborns into the system). A more refined model would consist of introducing such eventual recovery fluxes in the model, as well as the possibility of further relapses (the so called endogenous reactivation). These phenomenologies might be important mainly for diseases (like TB) for which the only feasible treatment in many areas consists of supplying large series of antibiotics. Thinking on the tuberculosis case, further refinements, like the inclusion of varieties of less infectious extrapulmonary diseases natmed95 (), could also have consequences on the disease’s dynamics.
Iii Structured Populations
iii.1 Dynamics
The previous system of differential equations describes the dynamics of the epidemics in the wellmixed case. However, as argued above, the number of contacts of a given individual in a population can vary, which is reflected in an heterogeneous distribution of the number of contacts in the system. To account for this fact, we next consider a structured population described by a connectivity distribution . The system of Eqs (1) has to be modified accordingly. Assuming that all individuals with the same number of contacts, i.e., belonging to the same connectivity class , are dynamically equivalent, the new system of differential equations are formulated for each degree class. Therefore, for a structured population, we have that:
(3) 
with:
(4) 
Moreover, it is convenient to express the previous equations in terms of densities, also defined within each connectivity class:
(5) 
so that the following closure relation for any value of is verified:
(6) 
On the other hand, the probability that any given link points to an infectious individual is given by:
(7) 
which leads to the following set of equations that describes the dynamics within each connectivity class:
(8) 
Finally, the number of individuals with connectivity evolves according to:
(9) 
At this point, and building on the previous equation, it is important to point out a feature of the model: the influence of the infection dynamics on the connectivity distribution . First, if we add the above equation for all , we obtain that the total population evolves as:
(10) 
However, if we substitute directly into Eq. (9) and assume to be constant, we would arrive to:
The last expression is only compatible with Eq. (10) under the unrealistic assumption that all connectivity classes have the same proportion of sick individuals. We must therefore assume that the distribution of connectivity is also a function of time: , and therefore:
(11) 
so, if we substitute Eq. (11) into Eq. (9) we get:
expression from which, if we replace from Eq. (10), we get the temporal evolution of as:
(12) 
where
(13) 
Reformulating the equations in terms of densities using the definitions of the densities given above, Eqs. (8) become
(14) 
iii.2 Evolution of the degree distribution
At this point it is appropriate to point out one aspect that will hinder any numerical characterization of the epidemic threshold. Our main goal will be to calculate both numerically and analytically the critical value beyond which the population presents an endemic proportion of sicks individuals. However, we expect to be dependent on the ratio , which is in turn a function of the connectivity distribution . The degree distribution, as previously shown, changes in time as the dynamics of infection progresses. As we should see, we can handle this time dependence analytically, but we should be forced to design a simulation method to account for the rate of births and deaths and the effects of these two processes on the degree distribution.
The aforementioned features might lead to a situation in which the infection dynamics would modify the underlying structure of the network through which the disease is being spread. Therefore, could also vary as one expects it to be intrinsically related to the first two moments of a seemingly timedependent degree distribution. The reason why we consider the distribution of contacts per unit time as heterogeneous, even for the current airbornetransmitted disease is based on the observation that the number of contacts a person can have per unit of time is subjected to two sources of heterogeneity. Firstly, what we can call geodemographic, macroscopic heterogeneity, in which the number of contacts depends on the population density in the region in which an individual inhabits. Secondly, at a more individual, microscopic level, the heterogeneity arises because the number of contacts depends, in a region of constant population density (i.e., a town or neighborhood in a city), on the daily activity pattern of the individual within that region. These two factors define, ultimately, the function . Having that said, the assumption implicitly incorporated in the first equation of system (8) does not hold. Note that this equation implies that the connectivity of individuals is hereditary and therefore that the number of births within each class equals the birth rate times the number of individuals within each class, .
The above situation would be equivalent to assume that the dynamics of the disease being studied is the only one that influences the demographic structure of a population, which is not true since it is clear that there are countless cultural, economic and social factors that ultimately define the above two levels of heterogeneity. We therefore assume in what follows that the newborns of each generation are distributed among the classes according to an invariant distribution function, which we further assume to be the initial degree distribution of the original network: . As we shall see, this assumption, besides being more plausible, has the advantage that makes the connectivity distribution to be roughly stable, and so will be the critical value .
So, we have the following reformulation of the system of differential equations (8):
(15) 
with the definition of the number of individuals in each class of connectivity:
(16) 
and inside each class:
(17) 
Now the total population within each connectivity class verifies:
(18) 
so that, if we add in , the last modification has no effect on the variation of the total volume of the population. The temporal evolution of the degree distribution is now given as:
(19) 
Finally, writing the equations in terms of the densities we get
(20) 
iii.3 Characterization of the equilibrium points.
The previous set of differential equations tells us how the different densities of interest evolves within each connectivity class. Their corresponding macroscopic quantities are defined as
(21)  
where and are the mean densities of healthy, latent and sick individuals, respectively.
Let us now go one step further and characterize the equilibrium points. The magnitudes of interest are the average densities, so that an equilibrium point must verify by definition:
We also impose a further constraint which is that the degree distribution of the network is stationary, that is:
At this point one must ask whether macroscopic stability also implies stability within each connectivity class. The answer is yes, if we also demand stability of the degree distribution. Admittedly, if we equate expression (19) to zero and solve for the stationary we get:
which shows that this value depends on the microscopic scale . Therefore, the stability of the degree distribution imposes a stationary condition on for all , which in its turns extends to the other densities and . Hence, we have:
The above condition is trivially satisfied for the solution , which leads to a degree distribution exactly as the initial distribution. We next analyze the stability of this solution, which shall allow us to characterize the epidemic threshold.
iii.4 Epidemic threshold
As stated before, in this section we will study the stability of the solution . At this point, as no latent or infected individuals are introduced in the network, the degree distribution does not change in time; so that . This situation allows to work with the system of differential equations given by (14) instead of working with the more general case given by the system (20).
iii.4.1 Case
For simplicity and to gain some preliminary insight into the problem, we first study the case , which means that there is no latent phase (i.e., the latent subpopulation disappears, ). Using we get,
where we have omitted temporal dependences, as we will do from now on. Looking for the stationary solution, we have that the condition implies:
from which the meaningful solution is the one with the negative sign. The previous expression is consistent with the meaning of since we recover the expected result when . Moreover, if we calculate the derivative with respect to we get:
which guarantees that will always be less than unity and therefore is a real, valid solution. The study of the value of in the steady state help us to identify the epidemic threshold. We write:
which, after substituting for its value, leads to:
The graphical interpretation of the above equation indicates that the existence of an equilibrium point in which is equivalent to the existence of a point at which crosses the bisector of the first quadrant. Evaluating the second derivative of one gets:
which ensures that the condition for the existence of such intersection is reduced to:
from which the epidemic threshold is derived as:
Note that apart from the factor , the previous result, formally coincides with the epidemic threshold of the SIR model.
iii.4.2 Case
This is a somewhat more involved case. For structured populations, the resolution of the system of differential equations (20) cannot be done explicitly. We next find the epidemic threshold for the case using two approaches. On one hand, we study the time derivative of . On the other hand, we will also make use of the singularity of the Jacobian at the point to argue that the expression for the critical threshold is given by:
(22) 
Time evolution of in populations with a low number of sicks.
A first approach to characterize the epidemic threshold in heterogeneous networks when is to study the sign of the derivative of at the onset of an epidemic outbreak. We consider an initially healthy population in which a small proportion of infectious individuals is introduced so that . The derivative of is:
which, after substitution of the values of the derivatives of and leads to:
At this point we make two simplifications. The first and most easily justifiable is to neglect the term . The second is related to the presence of in the above equation, that we have to transform in a dependency with respect to . Specifically, we assume to be sufficiently close to the stationary point as to be able to assume that the three derivatives vanish. In other words, and focusing our attention on latent and sick classes, we assume that:
from which:
which allows to express the derivative of as:
In the limit the third term within brackets is the ratio , from which the epidemic threshold condition may be derived as:
finally leading to the expected expression for the threshold:
(23) 
Analysis of the Jacobian
While for wellmixed populations the condition of singularity of the Jacobian allows to get the epidemic threshold in a straightforward way, for heterogeneous populations the analysis of the Jacobian is a difficult task because is a function of each and every one of the ’s. This translates into the need of calculating a determinant whose order is three times the number of connectivity classes. What we can reasonably do is to verify is the threshold condition is verified for systems in which there are two or three different connectivity classes. In the first case in which only two different classes of connectivity exist, the Jacobian is just a quite distasteful 6x6 determinant that, after some cumbersome and lengthy algebra, can be reduced to the expression:
which equated to zero leads again to the previously obtained expression for the epidemic threshold. If we instead consider a population with three degree classes, the algebraic complexity of the problem largely increases as we now have to solve a determinant of size 9x9. However, we can proceed as before getting the following expression for the Jacobian:
This leads us to the sensible conjecture that increasing the number of connectivity classes does not add new roots to the Jacobian, but it only would increase the degeneracy of the noninteresting solutions , and .
iii.5 Numerical simulations
When designing numerical simulations to inspect the dynamics of the system under study, we have two difficulties not previously addressed in the literature. These numerical issues with which we have to deal come from the fact that we have a system that is simultaneously open and structured. As a result of dealing with an open system, new individuals are being added to the population at a rate given by the birth rate. Additionally, these new individuals must enter the network of contacts with a predefined connectivity. While deciding how many nodes our new individuals connect to is not a problem, it certainly is to decide what are those nodes the newcomers will be linked to, as this will impact the degree distribution in a nontrivial way. This is an unavoidable numerical complication that we should face relentlessly if the analytical calculations are to be compared with Monte Carlo simulations.
To this end, we have adapted a simulation method based on transition probabilities first proposed in PRE () for SIR models in complex networks. The numerical approach considers all transitions between states that take place during the dynamical evolution of the subpopulations, defined by the system of differential equations (20). When dealing with structured populations, these transition rates depend, in general, on the connectivity class within which they occur. Moreover, within each class seven transitions are possible (see Fig. 1):
 

Birth of healthy individuals,
 

Natural death of healthy individuals,
 

Natural death of latently infected individuals,
 

Natural or diseaserelated death of sick individuals,
 

Transition from a healthy to the latent state,
 

Transition from a healthy to the sick (infectious) state,
 

Transition from a latent to the sick (infectious) state.
Each of these transitions is characterized by a characteristic transition rate that can be directly derived from the system of equations that characterizes the rate at which they occur within the class as:
Similarly, we define the sum of all these transition rates as the average rate at which one transition (of any kind) occurs:
(24) 
This average transition rate in its turn defines the characteristic or average time elapsed between any two consecutive transitions, the latter being defined as the inverse of :
(25) 
Given the previous definitions, the Monte Carlo algorithm is implemented in such a way that at each MC step (of duration ) one single transition takes place. Finally, the probability that a given transition actually happens, is calculated as:
(26) 
that determines which of all possible transitions is realized at each time step .
We have made extensive numerical simulations of the model starting from an initial population made up of individuals, whose network of contacts follows an initial degree distribution . Moreover, every newborn joins the system with a degree that verifies the same connectivity distribution. As for the values of the parameters of the dynamics, and thinking of typical values for persistent diseases, we have set the following values: , , , , and . Demographical parameters and are roughly those of a country like Spain, while the parameters and are in the range of typical values for the case of tuberculosis. has been chosen attending to numerical convenience (tuberculosis reaches a diseaserelated mortality rate that can be as large as 0.8). On the other hand, we note that a numerical criterion to define stationarity should also be adopted. In our simulations, we first let the system evolve for years and later take averages in a window of Monte Carlo steps (which corresponds, roughly, with a temporal lapse of 100 years), for the mean densities of healthy, latent and sick individuals defined in (21). This is a long enough time and ensures that all of the outputs to the left of the threshold are stabilized to the state (this is achieved almost surely already for years).
With the values for the parameters as specified above, the epidemic threshold is . In Fig. 2), we have plotted the stationary proportion of sick individuals for values of the probability of transmission in the interval . As can be seen from the figure, the numerical and analytical results are in a reasonable agreement, despite the numerical challenges of simulating an open system in which the dynamics evolves on top of a complex topology. This indicates that the numerical method is accurate enough as to be used in situations where analytical predictions are not at hand.
Analytical  Numerical  

1000  0.52903  0.46968 
10000  0.42559  0.37595 
50000  0.37620  0.33088 
100000  0.35854  0.31926 
However, it is possible to carry out numerical simulations using a variation of the algorithm above in order to improve the accuracy in the determination of the epidemic threshold . In this variation, instead of using a criterion for stationarity, we focus our attention in the vicinity of the critical value. More specifically, we first evaluate analytically the value for and start the simulation there (obviously, if we don’t have an analytical hint, the simulation can be started at any value of ). At each realization, we expect 6000 years for an eventual arrival of the system to the state . In the case that this state is not reached in that time, we assume that we are to the right of the critical point and so, we move to the left in just a little quantity . If, on the contrary, a minimum number of realizations (we used 10 in our simulations) stabilize at state , we assume to be to the left of the critical point and, consequently, we perform a switch to the right. Each time that such kind of flipflop algorithm (a sort of bisection method) changes direction, we divide by two the value of until the desired precision in is obtained. Using this numerical approach, we have numerically calculated the values of for different system sizes. The results are reported in Table 1. As expected from finitesize effects, the larger the size of the population, the smaller the absolute error between numerical and analytical thresholds is. Moreover, the larger the system size the smaller the epidemic threshold, which eventually should vanish in the thermodynamic limit.
Iv Conclusions.
We have discussed a model for the spreading of persistent infections in complex heterogeneous populations. The framework extends the epidemiological picture proposed in previous works. Our approach is particularly suited for diseases like Tuberculosis, which shows large latency periods. The latent period results from the dynamical equilibrium that is established between the bacterium and the host’s immune system, so that the host might not become infectious during its lifetime. These particular features makes it compulsory to work with an open system where newborns are continuously introduced in the population and individuals might die due to causes different from the disease itself. By assembling a model with all these ingredients, we have shown analytically that the epidemic threshold is proportional to the ratio between the first and second moments of the degree distribution. Therefore, our results point in the same direction that those obtained for the SIS and SIR model on top of the same topologies the virtually unbounded connectivity fluctuations play a key role in the infection dynamics enhancing the epidemic incidence and lowering the epidemic threshold.
From another point of view, we have developed a method well suited to numerically explore the dynamics of the system under study. In particular, we have been able to deal with the new challenge of having a system where the number of individuals in the population is not constant and, moreover, are connected following a given degree distribution (the wellmixed case is not challenging). Although we have applied the numerical approach to our particular system, it is worth stressing that it is general and can be applied to any problem for which transition rates between different classes and states are known. The results obtained agree well with the analytical estimates with the additional advantage that the whole phase diagram can be explored.
Finally, it is also worth mentioning that the model discussed here is probably the simplest one may devise for the spreading of persistent infections in structured populations. However, despite the recent progresses in modeling disease contagion dynamics and pandemic outbreaks, the kind of spreading phenomena analyzed here is one of the issues that have remained less explored. Our aim is to take a first step towards more realistic modeling of persistent diseases. We have left for future investigation possible extensions of the current model that take into account the influence on the dynamics of vaccination, prophylaxis and recovery rates, as well as the effects of genetic heterogeneity in the pathogen but also in the host MB ().
Acknowledgements.
We acknowledge financial support of the Government of Aragón through Grant PI038/08. This work has also been partially supported by MICINN through Grants FIS200801240 and FIS200913364C0201.References
 (1) R. M. Anderson, R. M. May, & B. Anderson, Infectious diseases of humans: Dynamics and Control (Oxford University Press, UK, Oxford, 1992).
 (2) D. J. Daley, & J. Gani, Epidemic Modelling (Cambridge University Press, UK, Cambridge, 1999).
 (3) J. D. Murray, Mathematical Biology (SpringerVerlag, Germany, Berlin, 2002).
 (4) L. Hufnagel, D. Brockmann, & T. Geisel, Proc. Nat. Acad. Sci. USA, 101, 1512415129 (2004).
 (5) R. Guimera, S. Mossa, A. Turtschi, & L. A. N. Amaral, Proc. Nat. Acad. Sci. USA, 102, 77947799 (2005).
 (6) V. Colizza, A. Barrat, M. Barthélemy, & A. Vespignani, Proc. Nat. Acad. Sci. USA, 103, 20152020 (2006).
 (7) S. Eubank, H. Guclu, V. S. AnilKumar, M. V. Marathe, A. Srinivasan, Z. Toroczkai, & N. Wang, Nature 429, 180 (2004).
 (8) V. Colizza, R. PastorSatorras, & A. Vespignani, Nat. Phys. 3, 276 (2007).
 (9) J. G. Gardeñes, V. Latora, Y. Moreno, & E. Profumo, Proc. Nat. Acad. Sci. USA 105, 1399 (2008).
 (10) V. Colizza, & A. Vespignani, J. Theor. Biol. 251, 450 (2008).
 (11) S. Meloni, A. Arenas, & Y. Moreno, Proc. Natl. Acad. Sci. USA 106, 16897 (2009).
 (12) M. E. J. Newman, SIAM Rev., 45, 167256 (2003).
 (13) S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, & D.U. Hwang, Phys. Rep., 424, 175308 (2006).
 (14) D. S. Callaway, M. E. J. Newman, S. H. Strogatz, & D. J. Watts, Phys. Rev. Lett., 85, 54685471 (2000).
 (15) R. Cohen, K. Erez, D. ben Avraham, & S. Havlin, Phys. Rev. Lett., 86, 36823685 (2001).
 (16) M. Faloutsos, P. Faloutsos, & C. Faloutsos, ACM SIGCOMM ’99, Comput. Commun. Rev. 29, 251 (1999).
 (17) G. Caldarelli, R. Marchetti, & L. Pietronero, Europhys. Lett., 52, 386 (2000).
 (18) R. Albert, H. Jeong, & A. L. Barabási, Nature, 401, 130 (1999).
 (19) S. H. Strogatz, Nature, 410, 268 (2001).
 (20) R. PastorSatorras, & A. Vespignani, Phys. Rev. Lett., 86, 32003203 (2001).
 (21) A. L. LLoyd, & R. M. May, Science, 292, 13161317 (2001).
 (22) Y. Moreno, R. PastorSatorras, & A. Vespignani, Eur. Phys. J. B, 26, 521–529 (2002).
 (23) M. E. J. Newman, Phys. Rev. E., 66, 016128 (2002).
 (24) B. M. Murphy, B.H. Singer, S. Anderson, & S. Kirschner, Mathematical Biosciences 180, 161 (2002).
 (25) L. Wilson, J. Hist. Med. Allied Sci. 45 366 (1990).
 (26) K. Styblo, J. Meijer, I. Sutherland, Bull. Int. Union Tuberc. 24(1), 137 (1969).
 (27) T. Daniel, J. Bates, K. Downes, History of Tuberculosis, in Tuberculosis: Pathogenesis, protection and control, Edited by B. R. Bloom (American Society for Microbiology, pp.1324, USA, Washington, 1994).
 (28) D. Bleed, C. Watt, C. Dye World Health Report 2001: Global Tuberculosis Control, Technical Report, World Health Organization, WHO/CDS/TB/2001.287, 2001. Available from http://www.who.int/gtb/publications/globrep01/index.html.
 (29) D. Bleed, C. Watt, C. Dye, Am. Rev. Respirat. Dis. 125, 8 (1982).
 (30) K. Styblo, Tuberculosis control and surveillance, in Recent Advances in Respiratory Medicine, Edited by D. Flenley, T. Petty (Vol. 4, pp. 77108, UK, London, 1986).
 (31) Blower, S. M., McLean, A.R., Porco, T. C., Small, P. M., Hopewell, P. C., Snchez, M. A., Moss A. R., Nature Medicine, 1(8), 815821 (1995).
 (32) Y. Moreno, J. B. Gómez, A. F. Pacheco, Phys. Rev. E 68 035103(R) (2003).