Mathematical modeling of the spread of infectious diseases

Mathematical modeling of the spread of infectious diseases A series of lectures given at PANDA, UNM Guillermo Abramson November 2001 This are informal notes, mostly based on the bibliography listed at the end and on recent papers in the field. 1 Introduction Outline of introduction: The practical use of these models is based on the fact that they can be kept realistic enough. Of course, this doesn’t mean to include every possible detail, but rather of every major mechanism. Bernoulli 1760 nonlinear ordinary differential equation. Effect of cowpox inoculation on the spread of smallpox. The initial plan is to discuss nonextended models in Section 1, and extended models in Section 2. These constitute the classical theory of epidemic spread. Section 3 introduces complex networks, and epidemics on complex networks in Section 4. 2 Simple epidemic models In this section we analyze a simple model of the temporal behavior of an infectious disease which is not extended in space. We can think of it as a kind of mean field model, or a model valid in a well-mixed population, which is usually far from the reality. However, several very interesting results can be derived. We consider that the total population is constant, and that a small number of infected individuals is introduced into a large population. The basic problem is the description of the time evolution of the infected group. For this we make some general reasonable assumptions. 1 SIR model Consider that the disease, after recovery, confers immunity (which includes the deaths, to have a constant population.) We can divide the population into three classes: S: the susceptibles, which can get the disease. I: the infected, who have the disease and can transmit it. R: the removed, either recovered and immune, or isolated, or sadly dead. It is customary to represent the progress of the disease schematically as: S → I → R, (1) and call these models SIR models. Let us call S(t), I(t) and R(t) the numbers (or the densities) of the members of each class. We now make assumptions concerning the transmission and incubation periods, in much the same way as if it were a chemical reaction: 1. The number of infected increases at a rate proportional to both the number of infected and the susceptible: rSI, with r > 0. The number of susceptibles decreases at this same rate. r is called the infection rate. 2. The rate of removal of infected to the R class is proportional to the number of infected only: aI, with a > 0. The number of removed increases at this same rate. a is called the removal rate. 3. The incubation time is negligible, so that a susceptible that catches the disease becomes infectious immediately. With these rates of change, and the further assumption of a well mixed population, where every pair of individuals has the same probability of getting into contact, we can write the following model, which is a kind of mean field model: dS dt dI dt dR dt = −rSI, (2) = rSI − aI, (3) = aI. (4) 2 We complete the formulation giving appropriate initial conditions: S(0) = S0 > 0, I(0) = I0 > 0, R(0) = 0. (5) The conservation of N = S + I + R is automatically ensured by the equations. This is a very simple model, but some rather important and general conclusions can be derived from it. Which are the important questions to make? First, given r, a, S0 and I0 , whether the infection will spread or not. And if it does, how will it develop in time, and when will it start to decline. From (??) we have at time zero: ¯ dI ¯¯ a (6) = I0 (rS0 − a) ≷ 0 if S0 ≷ = ρ. ¯ dt t=0 r And since form (??) we have that dS/dt ≤ 0, then we have always S < S0 and, if S0 < ρ, dI ≤ 0 for all t ≥ 0. (7) dt In this case, the number of infected remains below I0 and goes to zero as time goes to infinity. No epidemic can occur. On the other hand, if S0 > ρ, then I(t) starts increasing, and there is an epidemic (which is the very definition of epidemic). The picture is that of a threshold phenomenon: if S0 > Sc = ρ there is an epidemic while if S0 < Sc there is not. ρ is the relative removal rate. It reciprocal is the contact rate. R0 = rS0 /a is the reproduction rate of the infection. 1/a is the average infectious period. From (??) and (??) we can derive another important result: dI I(rS − a) ρ =− = −1 + , (I 6= 0). dS rSI S Integrating this equation gives (I, S) phase plane trajectories: I + S − ρ ln S = constant = I0 + S0 − ρ ln S0 . (8) (9) These trajectories need to be bounded by the triangle shown in the figure, due to the positiveness of the variables and the finite size of the system. [FIGURE (I,S)] Another important issue is how severe an epidemic will be. From (??) we can find the maximum of I, Imax , and it lies at S = ρ. Then, from (??): µ ¶ ρ Imax = N − ρ + ρ ln . (10) S0 3 In the figure we see that for any values of I0 and S0 > ρ, I increases from I0 to Imax , and an epidemic takes place. It need not be a severe epidemic, as is the case when I0 is near Imax . It can be also seen that there is no epidemic if S0 < ρ. Still more results can be derived from the model. We still have not used the equation for the removed. From (??) and (??) we have: dS dR S ρ ⇒ S = S0 exp[−R/ρ] ≥ S0 exp[−N/ρ] > 0 = − ⇒ 0 < S(∞) ≤ N. (11) In fact in the figure we see that 0 < S(∞) < ρ. And since I(∞) = 0 we have that R(∞) = N − S(∞). So, from (??): · ¸ · ¸ R(∞) N − S(∞) = S0 exp − (12) S(∞) = S0 exp − ρ ρ which is a transcendental equation from which S(∞) can be found as its positive root. Now, from (??) we can calculate the total number of susceptible individuals who get the disease in the full course of the epidemic: Itotal = I0 + S0 − S(∞). (13) An important consequence of this analysis, namely that I(t) → 0 and that S(t) → S(∞) > 0, is that the disease dies out due to a lack of infectives, and not from a lack of susceptibles. The epidemic does not grow unlimited to infect the whole population. There will always be some susceptibles that did not get the disease. Another general remark is that the threshold behavior is directly related to the relative removal rate ρ. For a given disease, the relative removal rate varies with the community, and it determines why an epidemic of a certain disease can occur in a certain community and not in another. For example, if the density of susceptibles is high (S0 is large) and the removal rate a is small (for ignorance, lack of adequate medical care, etc), then an epidemic is likely to occur. (Other things equal, a can be high if the disease is very serious and kills the infected fast, e.g. Ebola. In this case an epidemic is unlikely occur.) In real life epidemics, it is difficult to know how many new infective there are each day. Only those that are removed can be counted. To apply the 4 model to real situations, we need to know the number of removed per unit time, namely dR/dt. From (??) and (??) we get an equation for R alone: µ · ¸¶ dR R = aI = a(N − R − S) = a N − R − S0 exp − . (14) dt ρ Knowing the parameters, it is easy to compute the solution numerically. Unfortunately, the parameters are rarely known, and a fitting has to be made, assuming that the epidemic is well described by the model. In practice, if the epidemic is not large, R/ρ is small, and if it is very small we can expand the exponential in (??) to find: · µ ¶ ¸ S0 S0 R 2 dR = a N − S0 + −1 R− . (15) dt ρ 2ρ2 This approximate equation can be integrated to give: ·µ ¶ µ ¶¸ ρ2 S0 αat R(t) = − 1 + α tanh −φ S0 ρ 2 where α= "µ #1/2 ¶2 2S0 (N − S0 ) S0 −1 + , ρ ρ2 and tanh−1 ³ S0 ρ ´ −1 . α From this, the removal rate is found to be: ¶ µ dR aα2 ρ2 2 αat = sech −φ , dt 2S0 2 φ= (16) (17) (18) (19) which has only three parameters, aα2 ρ2 /2S0 , aα and φ. Bombay plague epidemic 1905-6 This epidemic lasted for almost a year. Since most of the infected died, we can fit the removal rate (??) to the death rate. The epidemic was not severe (compared to the population size) and the approximation performs very well. This was the original SIR model, performed by Kermak and McKendrick in 1927. [FIGURE BOMBAY] 5 Influenza epidemic in an English school 1978 In an English boarding school with a total of 763 boys, there was an epidemic of influenza from 22nd January to 4th February 1978. A total of 512 boys were put to bed during the epidemic, that seems to have started from a single infected boy. Here the epidemic is severe and we cannot use the approximation. However, we have I(t) directly from the data, and a numerical fit can be made to it. The figure shows the data and the numerical results. [FIGURE INFLUENZA] More complex models If a disease is not of short duration, then several changes need to be made to the SIR model. To start with, birth and death terms have to be included in the equation for the susceptibles. A natural death term has to be added also to the equation for the infected and removed classes. The resulting models can show oscillatory behavior that are called epidemic waves. Spatial epidemic waves can appear as an epidemic spreads geographically, as will be seen in the next section. A disease may have a latent or incubation time, when the susceptible has become infected but it is still not infectious. The incubation period for measles, for example, is 8–13 days. For AIDS, on the other hand, it can be anything from a few months to many years. This can be included in the model as a delay, or by introducing a new class, say E(t), in which the susceptible remains for a given time before moving to the I class. Such models give rise to integro-differential equations that also exhibit oscillatory behavior. Finally, age is often an important factor in a disease susceptibility and infectiousness. The model then becomes a partial differential system with independent variables time and age. There are many more modifications that can be incorporated in epidemic models, and these depend crucially on the disease. We will briefly mention some of these. A very important field is that of sexually transmitted diseases, which have several characteristics that make them different from other infectious diseases. For example, many times the carrier is asymptomatic until very late in the development of the disease (such as AIDS or Chlamydia trachomatis, which can produce sterility in women without showing any other symptom). Also, sexually transmitted diseases induce little or none immunity after infection, and the recovered become susceptible again. 6 For example, a considerably simplified model of gonorrhea, where the population is heterosexual and uniformly promiscuous, can be put schematically as: ✓ ❄ ✏ ✲ I ✟ ❍ ❨ ❍ ❨ ✟ ❍❍ ✟ ✙ ✟ ✙ ✲ I∗ S∗ ✻ ✒ S ✲ R ✲ R∗ ✑ where the starred variables stand for female, the contagion is represented by double arrows, and the recovered of both classes return to the corresponding susceptibles at some rate (this models are called SIRS). Further complications arise from the consideration of multigroup models, where the population is divided into more than three groups, such as: susceptible, symptomatic, treated infective, untreated infective, etc. In gastrointestinal parasitic infections, it is necessary to take into account still another mechanism, which is a very complex immunological response. This results into a delay of the immune response of the host, depending on the parasite density, something like: Z t P (τ )dτ, (20) E(t) = t−T where E is a variable describing the immune system and P is the parasite density. This diseases have a huge incidence in the developing world, with an estimated 1000 million of Ascaris lumbricoides and 500 million with Anclystoma duodenale and Nector americanus (figures of 1979). 3 Spread of epidemics The geographical spread of epidemics is even less well understood than their temporal development. And this kind of models are not only relevant to the spread of infectious diseases, but also other socially important phenomena such as rumors, misinformation, contagious habits such as drug abuse and fashions, and others. Diffusive model We consider a simpler model than the SIR of the previous section. Suppose that the population consists of only two classes: infectives I(x, t) and 7 susceptibles S(x, t) which are now functions of a spatial variable x as well as of time t. We model the spatial spread as a diffusive process, where both classes have the same diffusion coefficient D. We do not need to think that the individuals are actually diffusing; we can imagine them as fixed on a lattice, with contacts to their nearest neighbors, through which the disease propagates. The rates of transition from susceptibles to infectives, and of removal from infectives, are the same as in the mean field model. With this, we can write the following equations for the model: ∂S ∂t ∂I ∂t = −rIS + D∇2 S, (21) = rIS − aI + D∇2 I, (22) where r, a and D are positive constants. These are the same as equations (??) and (??) of the SIR model, with the addition of diffusion. The problem we are now interested in is the spatiotemporal spread of the disease when a number of infectives in introduced into an initially uniform population of susceptibles. In one dimension, and properly rescaling the variables to adimensionalize the equations, the model can be written as: ∂S ∂t ∂I ∂t ∂2S , ∂x2 ∂2I = IS − λI + 2 , ∂x = −IS + (23) (24) that has a single parameter λ = a/rS0 . The reproduction rate of the disease is 1/λ. For this model it is possible to analyze the conditions for the existence of traveling waves, that will represent an epidemic wave propagating into a susceptible population. We look for traveling waves by setting I(x, t) = I(z), S(x, t) = S(z), z = x − ct, (25) where c is the wave speed. Substituting these into (??) we have: I ′′ + cI ′ + I(S − λ) = 0, S ′′ + cS ′ − IS = 0, (26) where the prime is differentiation with respect to z. We will find the range of λ for which a solution exists with positive c and non-negative I and S such that I(−∞) = I(∞) = 0, 0 ≤ S(−∞) < S(∞) = 1. 8 (27) This condition means that a pulse of infective will propagate through the susceptible population, reducing its density. [FIGURE TRAVELING WAVE] The system (??) is fourth order. To go a little further analytically we can linearize the equation for I at the leading edge of the wave, where S → 1 and I → 0: I ′′ + cI ′ + (1 − λ)I ≈ 0, (28) which we can solve as: I(z) ≈ exp h³ ´ i p −c ± c2 − 4(1 − λ) z/2 . (29) This solution cannot oscillate about I = 0, since I(z) → 0. Then, the wave speed and λ must satisfy: √ c ≥ 2 1 − λ, λ < 1. (30) If λ > 1 there is no wave solution, so this is the threshold condition for the propagation of an epidemic wave. Going back to dimensional terms, this means: a < 1. (31) λ= rS0 The threshold result (??) has important consequences. We can see that there is a minimum critical density of susceptibles, Sc = a/r, for an epidemic wave to occur. Correspondingly, for a given population S0 and a given mortality rate a, there is a critical transmission coefficient rc = a/S0 which, if not exceeded, prevents an epidemic. There is also a threshold mortality, ac = rS0 which, if exceeded, prevents an epidemic. We see that the more rapidly fatal a disease is, the less chance there is to have an epidemic wave through a population. All of these observations have consequences for control programs. The susceptible density can be reduced with vaccination. The transmission can be reduced by isolation, medical removal, etc, to stop an epidemic. The speed of the waves in dimensional terms results to be: ¶ µ p a 1/2 V = 2 rS0 D 1 − . rS0 (32) This is the minimum speed which turns out to be the observed speed due to dynamical selection of the marginally stable state. Similarly, the traveling wave S(z) can be calculated. 9 The Black Death in Europe 1347-1350 The catastrophic pandemic of plague in the 14th century is so fascinating that it is worth summarizing some facts about it. The Black Death, as it was known, was mainly bubonic plague, caused by Bacillus pestis and transmitted by fleas from rats to humans. The disease was introduced to Italy in December 1347, brought by a single ship from the East. During the next years, it spread through Europe at a speed of 200–400 miles per year. About a quarter of the population died, and 80% of the infected died within 2 or 3 days. The figure shows the spatio-temporal spread of the wave front. [FIGURE BLACK DEATH] After the Black death had passed, around 1350, there was a second outbreak in Germany in 1356. From then on, periodic outbreaks seem to occur every few years, even though none as severe as the Black Death. There is a belief that plague ceased to be a problem after the Great Plague of London in 16651 . However, the truth is far from this. The last plague pandemic started in Yunnan, China, in 1850 and officially lasted until 1959, after killing more than 13 million and affecting almost any part of the world. Plague was bought by ship from the Far East to the northwest of America around 1900. There was a three year epidemic in San Francisco with 200 deaths, starting just after the earthquake in 1906. Since then, the bacillus has been spreading steadily towards the east, moving at about 35 miles a year, carried by many small mammals (squirrels, rabbits, prairie dogs, coyotes, bobcats,...), not only rats. As a result, the western USA and in particular New Mexico is one of the two largest residual foci of the plague. The number of cases in New Mexico have been 9 in 1998, 0 in 1997, 2 in 1996, 4 in 1995, an average of 7.8 per year in the period 1970–1991. In the USA, there have been 341 cases in the period 1970–1995. What could happen when the front reaches the big urban centers of the east coast, with their huge populations of rats, nobody knows. Returning to the Black Death, it is not easy to estimate the parameters of the model. There were about 85,000,000 people in Europe in 1347, which gives a density S0 ≈ 50/mile2 . The transmission and diffusion coefficients are more difficult to estimate. Suppose that news and gossip, with the limited communications of the time, traveled at 100 miles/year, and that the disease travels with these. If this were a diffusion time, we have that 1 By the way, due to the plague Newton returned from Cambridge to Woolsthorpe after his graduation in 1665, and stayed there until 1667. That was his annus mirabilis during which he made most of his achievements in mathematics, mechanics and optics. 10 the time to cover a distance L is o(L2 /D). This gives D ≈ 104 miles2 /year. For transmission a value of r ≈ 0.4 miles2 /year has been estimated, and for mortality rate, based on an infectious period of two weeks, a ≈ 15/year. These give λ ≈ 0.75. From the expression for the velocity of the front, the estimation is: V ≈ 140 miles/year. The uncertainties are very large, but the comparison with the data is reasonable good for such a simple model. The spatial spread of rabies Rabies is a very widespread disease and epidemics (“epizootics” is the proper term) are common. During the last several hundreds years there has been periodic epidemics in Europe. The current one originated in Poland in 1939, and has been moving westward at about 30–60 miles/year. The red fox is the main carrier and victim, but affects many mammals, including bats. There is also an epidemic wave in the east coast of America, moving northward, where the main carrier is the raccoon. Rabies is a viral infection of the central nervous system, transmitted by direct contact. The incidence on man is now rare, with only a few death a year in Europe and the USA, but many more in undeveloped countries. It is a particularly severe disease, of which not a single case of recovery is known after its clinical manifestation (during the incubation, vaccines are 100% effective). Again, the spatial spread of this epidemic is a very complex problem, but certain assumptions can be made. Foxes carry more than 70% of the disease in Europe, so a single species with subpopulations can be considered. Foxes are territorial, and their territories do not overlap. Rabies induces behavioral changes such as dissorientation and aggressivity, so that diffusion only in the infected might be a good transport model for these two reasons. Rabies is invariably fatal, so that infected are removed from the system at a constant rate per fox. We can write the following model, which is very similar to (and simpler than) the previous SI model: ∂S ∂t ∂I ∂t = −IS, = IS − λI + (33) ∂2I , ∂x2 (34) already in dimensionless variables. Also, λ = a/rS0 measures the mortality rate as compared with the contact rate. A traveling wave front can be found 11 with the boundary conditions: S(∞) = 1, S ′ (−∞) = 0, I(∞) = I(−∞) = 0. (35) The condition on S(−∞) is in its derivative, since we expect a number of survivors, already undetermined. Just as in the previous model, the system that describes the wave front becomes: cS ′ = IS, ′′ (36) ′ I + cI + I(S − λ) = 0. (37) Linearizing around I = 0 and S = 1 just as before we find the same threshold behavior: √ c ≥ 2 1 − λ, λ < 1. (38) But now the equations are simpler and the analysis can proceed a little further. From (??) we have I = cS ′ /S, which can be substituted into (??) giving: cS ′ (S − λ) I ′′ + cI ′ + = 0. (39) S This can be integrated, giving: I ′ + cI + cS − cλ ln S = constant. (40) We can determine the constant by using the boundary condition at +∞, where S = 1 and I = 0 (and I ′ = 0). This gives the constant equal to c. Now, at −∞, we have the following equation that determines the surviving susceptible population S1 after the passage of the epizootic wave: S1 − λ ln S1 = 1, λ < 1, S1 = S(−∞), (41) which turns out to be independent of c. It is easier to analyze if put in the form: S1 − 1 = λ ⇒ 0 < S1 < λ < 1. (42) ln S1 S1 can be easily plotted against λ (indeed the inverse of the plot of λ(S1 )). The parameter λ results in a measure of the severity of the epidemic: the smaller λ, the fewer susceptibles survive. [FIGURE RABIES DATA] The figure shows real epizootic data based on the susceptible fox population in France from the wave front of 1977. It can be seen that there are oscillations after the front. It is clear that the disease frees territories, that 12 become available for new young foxes, since the time scale of the wave front is considerably less than the period of the oscillations. Then, it is reasonable to model this as a population growth of the susceptibles, for example with a term such as bS(1 − S/S0 ). Proper adimensionalization gives the model: ∂S ∂t ∂I ∂t = −IS + bS(1 − S), = I(S − λ) + ∂2I . ∂x2 (43) (44) This model has traveling wave solutions of precisely the form of the epizootic data, with a great initial kill followed by oscillations around a surviving fraction. We can see that very simple models are able to capture the major phenomena observed in the propagation of epidemic waves. A more complicated model with three classes (susceptible, infected and rabid) can be justified because the disease has a long incubation time (150 days) during which the infected fox does not shows symptoms and does consequently does not change behavior and does not infect, compared to the clinical development which takes a few days until the fox dies. With this model a quantitative prediction can be made for the control of the epizootic wave by means of a barrier of reduced S. This has been done with success in Denmark, and with less success in Italy and Switzerland. A two dimensional model with good quantitative predictions has also been carried out by Murray (the prediction usually compared to the data is the speed of the epidemic front). (Murray is concerned with the arrival of rabies to England, which is currently free of it, and which has a larger density of red foxes than the rest of Europe.) 4 Complex networks The models of the previous sections constitute the “classical” theory of epidemic propagation. Even though they are able to describe some aspects of the propagation of infectious disease through populations, it is clear that human societies are much more complex than well mixed populations, or two dimensional lattices, or sets of diffusing individuals. This is also true for other complex systems where a set of individuals form a kind of network through which some field can propagate in the same way as an infectious disease. The Internet and the World Wide Web are examples of such networks. The Internet is the physical network of routers and computers linked by telecommunication means. The WWW is the logical web of webpages 13 connected by hyperlinks. The cell can also be thought as a network of chemicals connected by chemical reactions. These are only a few examples of the diversity of this kinds of structure, and you can surely think of many more: telephone calls networks, citation networks of scientific publications, movie actors collaboration, power grids, air travel lines,. . . One that is particularly fascinating for me is the network formed by the words in English that, according to the Merriam-Webster Dictionary, have synonyms. There are 23 thousand such words, of which 22 thousand form a giant cluster. This is already remarkable, that they are almost all synonyms of one another! But this is not all. The average path length is L = 4.5. So you can go from any word to any other, with any completely different meaning, through a chain of just four synonyms! Traditionally, the study of this structures is the field of graph theory. While regular graphs were its focus for many years, in the ’50s Erdös and Rényi pioneered the study of large scale random graphs. In their model, N nodes are given, and every pair of nodes is connected at random with probability p, giving a graph with approximately pN (N − 1)/2 edges distributed at random. Are the complex networks of human society random? Well, to some extent they are, but our intuition tells us that they also have some organization that the ER model will not be able to capture. In the last few years there have been developments in the computarization of many fields, and also in the computer power, enabling the analysis of many very large complex webs in different fields. And motivated by this, there have been many proposals to their analysis. Three main concepts have emerged as most prominent in the study of complex networks: small world, clusterization, and degree distribution. Small world We all know the expression “it’s a small world,” and in other languages there exist also similar expressions to reflect the fact that in the network of human acquaintances, we seem to be rather closely connected to one another. The psychologist Stanley Milgram in 1967 made probably the first quantitative study of this phenomenon. He prepared a number of letters addressed to a person he knew in Boston, and distributed them to a random set of people in Nebraska (as far away from Boston as you can tell). His instructions were that the letters had to be passed from person to person, only to close acquaintances, to someone that the current holder of the letter believed to be closer to the addressee than himself. It turned out that 14 an average of six steps were enough for a letter to get from Nebraska to Milgram’s friend in Boston. The situation has been labeled “six degrees of separation,” and has find its way into a theater play, a Hollywood movie (featuring Donald Sutherland and Will Smith) and several rather frivolous games about the degree of separation between celebrities such as Kevin Bacon or Monica Lewinsky. You can immediately note that, despite the apparent frivolity of all this, there are some very important consequences to the small world phenomenon. These are relevant to the propagation of communications, rumors, fashions, news and, of course, infectious diseases, which spread through a network of people by contact between a pair. The small world phenomenon can have dramatic consequences, since a disease will spread much faster than in a regular lattice, where the typical distance between a pair is relatively long. Random graphs are “small” in this sense. Their average path distance between any pair of nodes scales approximately with the logarithm of its size: L ∼ ln(N )/ ln(hki) (where hki is the average coordination number, or degree, or connectivity). (This is valid for not small p, specifically pN > 1.) Clusterization A common feature of social networks, as opposed to random ER networks, is that the circle of acquaintances of a person tends to overlap to the circle of his acquaintances. Your friends’ friends tend to be also your friends. This produces a varying degree of clusterization, or “cliquishness” in social networks. This side of the phenomenon, which is not present in random graphs, was recognized first by Watts and Strogatz, who incorporated it in their small world networks, in an article that appeared in Nature in 1998 and that eventually was responsible for bringing the whole field to the attention of the physicists’ community. It is possible to define a clustering coefficient C as the average fraction of pairs of neighbors who are also neighbors of each other. That is, suppose that the node i has ki edges that link it to other nodes. If all the neighbors of i were neighbors to each other, there would be ki (ki − 1)/2 edges between them. Suppose that, instead, there are Ei . Then the clusterization around i is defined as: 2Ei Ci = . (45) ki (ki − 1) The clustering coefficient of the network is then defined as the average of Ci over the network. In a fully connected network, C = 1. If you want to think of it in probabilistic terms, C is the probability that to nodes that are 15 neighbors to a particular node, are also neighbors of each other. In an ER random graph, it is C = p = hki/N , which is very small for large networks. (Consider the graph formed by a site and its neighbors. The probability that any pair of these are connected is just p, as is for any other pair of the network.) In real life networks, it has been found that C is significatively less than 1, but much greater than o(N −1 ). The small world networks of Watts and Strogatz show this behavior. Degree distribution In a real life network not all the nodes have the same number of edges (called the degree or the connectivity of the node). It is reasonable then to define a degree distribution P (k) to characterize the way this magnitude is spread. This gives the probability that a randomly selected node has exactly k edges. Since an ER random graph is built precisely by assigning edges at random, most of the nodes have a degree near the average hki, and P (k) is a Poisson distribution peaked around this value. It has been shown that some real networks have exponential tails at large k, while significantly different from Poisson distributions, and still others show power law tails: P (k) ∼ k −γ and are called scale-free networks. Small worlds of Watts and Strogatz Watts and Strogatz introduced a procedure to interpolate, with a single parameter of disorder, between an ordered lattice and a random graph. It consists of taking, as a starting point, a regular lattice of coordination number k (usually a ring, with k/2 links to each side). The clustering coefficient of this lattice can be calculated exactly, giving C= 3(k − 2) , 4(k − 1) (46) valid for k < 2N/3, that tends to 3/4 in the limit of large k. This structures do not show the small world effect, since the node to node distance increases linearly with N . Then, some degree of randomness is added to this basic low-dimensional lattice. This is done by going through each of the nodes and, for each of its edges that points to the right (say), rewire it at random to any node of the network, with probability p. For small p the new graph is almost regular, but with a few links that reach far away in the system. The average coordination number of the network remains k. The procedure is shown schematically in Figure ??, for a small system with k = 4. 16 Figure 1: The rewiring procedure to build small worlds starting from a regular lattice. It turns out that, for these systems, with small p, the average node to node distance is essentially the same as that of a comparable random graph, while its clustering coefficient (essentially 3/4 at low p) is many orders of magnitude greater. For example, for N = 1000 and k = 10, the average distance is L = 3.2 for a random graph. A small world network with p = 1/4 has L = 3.6, compared to L = 50 for the regular lattice. A second version of the model has been proposed, in which instead of rewiring links you add links (called shortcuts) with probability p to an underlying regular lattice. The reason for the introduction of this is that it cannot become disconnected by the procedure, as is the case with the initial model. Some results have been proved for this networks, and many more have been explored numerically. Lets summarize some of them. Most of the work concentrated at first on the properties of the average path length, which was observed to begin to decrease once p > 2/N k which guarantees the existence of at least one shortcut. This can be interpreted as the existence of a p-dependent transition length ξ such that, for a fixed value of k, if N < ξ then L ∼ N , and if N > ξ then L ∼ ln N . It was proposed that the characteristic path length scales with this as: µ ¶ N , (47) L(N, p) ∼ ξF ξ 17 where F (x) = ½ x if x ≪ 1, ln(x) if x ≫ 1, (48) and, from both numerical and analytical results: ξ= 1 . pk (49) The divergence of ξ as p → 0 is somewhat similar to a critical point. There is a kind of phase transition at p = 0 from the behavior L ∼ N to L ∼ ln N . In other words, for an infinitely large system, an infinitesimal amount of disorder induces the small world behavior. It was finally shown, after some controversies, that the characteristic path length obeys the scaling (in dimension d): L(N, p) ∼ where F (x) = ½ N F (pkN d ), k 1 if x ≪ 1, ln x/x if x ≫ 1. (50) (51) In addition to a short average path length, small world networks have a relatively high clustering coefficient for a wide range of values of p. The clusterization of the model depends only on k, and not on N . The value of C(p) remains close to C(0) up to relatively high values of p. A typical situation is shown in Figure ??, where the average path length and the clusterization are displayed simultaneously. The intermediate range, where the clusterization is large—as in a regular lattice—and the average path length is short—as in random graphs—are the small worlds. Scale free networks Scale free networks, as proposed by Barabási and Albert, put the emphasis on the evolution, the construction, of the networks. Their original model, just as many real networks, does not start with a given set of N nodes and some links between them. Instead, starting from a small nucleus, the network is grown by the addition of new nodes and links according to some rules that, presumably, mimic the mechanisms by which real networks grow. The main rule in this models is that of preferential attachment, such that the likelihood of connecting to a node depends on the node’s degree, for example according to a probability: ki Π(ki ) = P . j kj 18 (52) Figure 2: Characteristic path length and clustering coefficient (both normalized) of the small worlds of Watts and Strogatz. After t time steps of adding one node with m edges at a time, the network has grown to a system with N = t + m0 nodes and mt edges. A very simple continuous model is able to show the behavior of the degree distribution. The degree ki of node i will increase every time a new node enters the system and links to it, something that happens with probability Π(ki ). Suppose that ki is a continuous variable of time. Then we can write the rate equation ki ∂ki = mΠ(ki ) = m PN −1 , (53) ∂t j=1 kj where the sum inPthe denominator runs over all the nodes except the one −1 just introduced: N j=1 kj = 2mt − m, so that: ki ∂ki = . ∂t 2t (54) With the initial condition that node i had ki (ti ) = m at the time of its 19 introduction, the differential equation can be solved as: µ ¶1/2 t . ki (t) = m ti (55) Now, before proceeding to the probability of having exactly k edges, we start with the probability that a node has a degree ki (t) smaller than k. This can be written as: µ ¶ m2 t P (ki (t) < k) = P ti > 2 . (56) k Assuming that we add nodes at equal time intervals, the ti have a uniform probability 1 1 P (ti ) = = . (57) N t + m0 Using this into the previous one gives: ¶ µ ¶ µ m2 t m2 t m2 t . = 1 − P ti ≤ 2 =1− 2 P ti > 2 k k k (t + m0 ) (58) And the degree distribution can be obtained by using: P (k) = 2m2t 1 ∂P (ki (t) < k) = ∂k m0 + t k 2+1 (59) that gives the asymptotic (t → ∞) behavior P (k) ∼ 2m2 k −γ , with γ = 2 + 1 = 3, (60) with exhibits an exponent that is independent of m. This power law is universal, which imposes some limitations to model real scale free networks, which have been shown to have different exponents. In fact the exponent can be made to change by changing the form of the preferential attachment probability Π(ki ), which in fact seems to be the case in some real networks. How important are the roles of growth and preferential attachment in the emergence of a power law degree distribution? Without the preferential attachment (connecting the new node with uniform probability to the existing nodes of the network), the continuous model gives P (k) exponential. Without the growing (starting with N nodes and no edges, and adding these according to the preferential attachment rule) P (k) is not stationary. At short times it is a power law, evolving then to a Gaussian (eventually, all the nodes become connected). 20 Scale-free networks are “small.” It has been shown that the average path length of the model increases logarithmically with N . This result is based on an empirical fit, and no analytical arguments exist for it. The clustering behavior of these models is rather different than that of the small world models. Instead of being independent of N , it decays following approximately a power law C ∼ N −0.75 , which is nevertheless slower than the behavior of random graphs, for which C = hki/N . Also, the linear growth of this model can be modified, with effects on the network topology. For example, an increase in the average degree of networks such as the Internet, the WWW, the network of coauthorship of scientific papers, suggests that in real systems the number of edges increases faster than the number of nodes, something called accelerated growth. 5 Epidemics on complex networks All kinds of networks offer a support for interesting dynamical processes, such as the spread of information, ideas, fashions, energy, disease. Games, cellular automata, synchronization of oscillators, percolation, random walks, error tolerance, would be among the interesting fields. The networks can be much more complex that what we have just see. They can be directed, weighted, evolving. . . From all the many fascinating aspects that surround the concept of complex networks, we will discuss here only the phenomenon of disease spreading. In particualr we will be addressing the following questions: How does the dynamics of an infectious disease depend on the structure of a population? Are there consequences to the strategies of control of an epidemic that depend on the topological features of the network? Epidemics in a small world network We will use here a discrete version of an SIRS model. The disease has three stages: susceptible (S), infected (I), and refractory (R). An element of the population is a node of the network described by a single dynamical variable adopting one of these three values. Susceptible individuals can pass to the infected state by contagion by an infected one. Once infected, the state of the individual follows a deterministic trajectory in discrete time through the stages of the disease: S → I → R → S, (61) spending a time τI in the infected state, and a time τR in the refractory state, after which he returns to the susceptible state. The disease is not supposed 21 to confer a permanent immunity. This is, effectively, an excitable system, which is known to show space-time oscillations in extended systems, similar to those of a chemical system like the Belousov-Zhabotinskii reaction. The contagion of a susceptible element occurs stochastically at a local level. If an element i is susceptible, and if it has ki neighbors, of which kinf are infected, then i will become infected with probability pinf = kinf /ki . Observe that i will become infected with probability 1 if all its neighbors are infected. Besides this parameter-free mechanism, there may be other reasonable choices. For example, if the susceptible had a probability q of contagion with each infected neighbor, we would have a probability of infection [1 − (1 − q)kinf ]. The network is put in an initial state composed almost exclusively of susceptibles, with a few infected randomly distributed. The temporal behavior of this system is then analyzed for different values of p, representing societies with different topologies. Figure ?? shows part of the time evolution, once a stationary state is achieved, of the fraction of infected nodes for three values of p. We can see clearly a transition from an endemic situation to an oscillatory one. At p = 0.01 (top), where the network is nearly a regular lattice, the stationary state is a fixed point, with fluctuations. The situation corresponds to that of an endemic infection, with a low and persistent fraction of infected individuals. At high values of p—like the case with p = 0.9 shown in the figure (bottom)—large amplitude, self-sustained oscillations develop. The situation is almost periodic, with a very well defined period and small fluctuations in amplitude. The period is slightly longer than τ0 , since it includes the average time that a susceptible individual remains at state S, before being infected. A mean field model of the system, expected to resemble the behavior at p = 1, can be easily shown to exhibit these oscillations. The transition between both behaviors is apparent—in this relatively small system—at the intermediate value of disorder p = 0.2, shown in the middle curve. Here a low amplitude periodic pattern can be seen, appearing and disappearing again in a background of strong fluctuations. Moreover, the mean value of infection is seen to grow with p. The formation of persistent oscillations corresponds to the spontaneous synchronization of a significant fraction of the elements in the system. We can describe the cycle of the disease of each node i as a phase variable τi (t) = 0, 1, . . . , τI + τR = τ0 . In the oscillatory phase the phases become synchronized, and they go over the disease process together, becoming ill at the same time, and recovering at the same time. We can characterize this 22 ✁✄ ☎ ✆ ✝✞✝ ✡ n inf (t) ✁✂ ✁ ✁✄ ☎ ✆ ✝✞✠ ✁✂ ✁ ✁✄ ☎ ✆ ✝ ✞✟ ✁✂ ✁ 5200 5300 5400 5500 5600 t Figure 3: Fraction of infected elements as a function of time. Three time series are shown, corresponding to different values of the disorder parameter p, as shown in the legends. Other parameters are: N = 104 , K = 3, τI = 4, τR = 9, Ninf (0) = 0.1. behavior with a synchronization parameter: ¯ ¯ ¯ ¯ N X ¯ ¯1 i φj (t) ¯ ¯ e σ(t) = ¯ ¯, ¯ ¯N (62) j=1 where φj = 2π(τj − 1)/τ0 is a geometrical phase corresponding to τj . When the system is not synchronized, the phases are widely spread in the cycle and the complex numbers eiφ are correspondingly spread in the unit circle. In this situation σ is small. On the other hand, when a significant part of the elements are synchronized in the cycle, σ is large. The synchronization would be strictly σ = 1 if all the elements were at the same state at the same time. However, such a state would end up in Ninf = 0 after τ0 time steps, and the epidemic would end since the system is closed and no spontaneous infection of susceptibles is being taken into account in the model. 23 0.4 0.3 ✁ ✂✄☎ ✆✝✄ ✞✟✠✡☛☞✠✌☛✍✎✏✑ I 0.2 ✁ ✂✄✒ ✆✝✄ ✞✟✠✡☛☞✠✌☛✍✎✏✑ ✁ ✂✄✓ ✆✂✄ ✞✟✠✡☛☞✠✌☛✍✎✏✑ ✁ ✂✄✔ ✆✕ ✞✟✠✡☛☞✠✌☛✍✎✏✑ 0.1 0.0 0.2 0.4 0.6 0.8 1.0 p Figure 4: Synchronization of the system as a function of the disorder parameter p. Each point is a time and realizations average. Other parameters are: K = 3, τI = 4, τR = 9, Ninf (0) = 0.1. The synchronization parameter provides a good order parameter to study the system as we increase the system size. Figure ?? shows the synchronization parameter σ, obtained as a time and realizations average, corresponding to system sizes N = 103 , 104 , 105 and 106 . A transition in the synchronization can be seen as p runs from 0 to 1. The transition becomes sharper for large systems, at a value of the disorder parameter pc ≈ 0.4. It is worthwhile to remark that the transition to synchronization occurs as a function of the structure of the network, contrasting the phenomenon of synchronization as a function of the strength of the interaction, as in other systems of coupled oscillators. This behavior is observed for a wide range of values of τI , τR , K and the initial fraction of infected. Effect of immunization Vaccination is one of the most effective strategies against the spread of 24 Figure 5: Reported measles cases in four countries. From Cliff & Haggett, Sci. Am. 250, May, 118 (1984) epidemic (and still many airborne diseases, such as measles, remain endemic despite vaccination). After the many socioeonomical questions for a specific vaccine have been answered, such as its safety, its cost, its effectiveness against the disease, etc., there is still one important question: what proportion of the population must be immunized in order to control the disease as effectively as possible? Certainly the structure of the social links again must play an important role. In scale free networks, that display strong inhomogeneity in the degree distribution, it has been shown that vaccination of those few individuals with maximal connectivity produces the best results [?]. In small world networks the degree distribution is much more homogeneous, but again targeted vaccination of the best connected individuals produces the best results. Moreover, a threshold of immunization exists above which there is a drastic reduction on the impact of the disease [?]. 25 The disease is modeled as an SIR system, implying that there is no transition R → S, the removed being either permanently immune or dead. Immunization is applied before the evolution starts to a fraction ρ of the population, by moving some nodes to the R state, following one of two possible strategies. In random immunization these are chosen uniformly at random in the population. In targeted immunization those individuals with the greatest number of links are chosen. The initial condition contains only one infected node. The evolution of the disease consists of the spread of an infected phase until a fraction of the population has been infected, and some fraction remains susceptible (just as in the SIR model of Section 1). Figure ?? (top) shows the fraction of infected to non-immunized people during the whole course of the epidemic for different levels of immunization. In the figure a transition at a finite value of p can be seen at each level of immunization, above which a finite fraction of the population is affected. For targeted immunization a similar picture arises, albeit quantitatively different, as we shall see. We see that there is a level of immunization that almost suppresses the disease in any network. It is worth observing also that the critical p grows with the immunization level ρ. The meaning of this is that, for small values of the disorder in the network, a little immunization has very strong effects in the control of the epidemic. Figure ?? (bottom) shows a comparison between the random and the targeted strategies (at p = 1 as representative of the maximum infection level). Targeted vaccination appears as a more effective strategy. It has to be kept in mind that the targeted immunization requires the identification of the nodes with the greatest connectivity, which is a time consuming procedure not only in the simulations but also, presumably, in a real situation as well. The important practical consequences of this, however, is that for both strategies there exists a critical value of the immunization level that guarantees that an outbreak of infection will remain localized and will not spread into the population. This can be translated for example as: for a given population with a fixed small world structure, near the critical situation a small increase in the budget for immunization can result in a drastic reduction of the impact of the disease. It is interesting to note a similitude with another propagation phenomenon in a social network, namely the propagation of a rumor. The fraction that remains susceptible is equivalent to a fraction of the population that never hears the rumor, a well known fact in the theory of rumor. Even in every day life we most certainly know some people that complain 26 1.0 random immunization 0.8 ρ = 0.0 ρ = 0.1 ρ = 0.2 0.6 ρ = 0.3 ρ = 0.4 r 0.4 0.2 0.0 0.01 0.1 1 p 1.0 random immunization targeted immunization r (p =1) 0.8 0.6 0.4 0.2 0.0 0.0 0.1 0.2 ρ 0.3 0.4 0.5 Figure 6: Top: Fraction r = NR /(1 − ρ)N of the non-vaccinated population that becomes infected during the disease propagation, as a function of the small-world randomness p, for various levels ρ of random immunization. Bottom: Fraction r for maximal randomness, p = 1, as a function of the immunization level ρ, for random and targeted immunization. 27 about never hearing of anything! This problem has also been studied [?]. Epidemics in a scale free network How is the situation in scale-free networks? There is a recent study [?] on the evolution of an SIS model on a scale-free network that shows the very interesting and relevant phenomenon of the absence of an epidemic threshold. According to this, scale-free networks such as the WWW and perhaps also some social networks might be particularly prone to the propagation of infectious agents. SIS is the standard model for the analysis of computer virus infections, which is a field that has borrowed a lot from the epidemiology of the life sciences, as it is natural. At each time step a susceptible node is infected with probability ν if it is connected to one or more infected nodes. Infected nodes are cured—and become susceptible again—at a rate δ. A relative parameter is defined as λ = ν/δ and δ set to 1 for good. In regular lattices, as well as in random graphs, there is an epidemic threshold λc , much in the spirit of the models we have already discussed. If λ ≥ λc the infection spreads and becomes persistent. If λ < λc , instead, the infection disappears exponentially. The analysis of real informatic virus epidemics shows that the viral density saturates to a very low level, affecting a tiny fraction of the computers, something that could happen if all the different computer virus had their λ tuned just an infinitesimal above λc . (Oh, yes, it could just be a selforganized-criticality phenomenon.) To model this phenomenon, the SIS mechanism was simulated on a scale-free network of the kind discussed in Section 3. Initially, half of the network was infected, and a steady state of density ρ of infected nodes was eventually reached. Figure ?? shows the average asymptotic persistence of the infection as λ decreases. There is no threshold behavior: for any finite value of λ the virus survives in a sufficiently large network. The reason for this is that the unbounded fluctuations in the connectivity of a scalefree network (hk 2 i = ∞) play the role of an infinite connectivity (while in a “normal” network, the higher the average connectivity, the lower the threshold). A mean field calculation is able to describe this result. The rate equation for the density of infected nodes with connectivity k can be written as: ∂ρk (t) = −ρk + λk[1 − ρk (t)]Θ(λ). (63) ∂t Here, the creation term considers the probability that a node with k links is 28 −2 ρ 10 10 10 −3 −4 7 12 17 22 1/λ Figure 7: Persistence ρ as a function of 1/λ for different network sizes, up to several million nodes. The full line is a fit to the form ρ ∼ exp(−C/λ). healthy [1−ρk (t)] get infected form an infected node. The probability of this is proportional to the infection rate λ, the number of its own connections k, and the probability that it points to an infected node Θ. The stationary value of the infected density is found setting the time derivative equal to zero, and: kλΘ(λ) . (64) ρk = 1 + kλΘ(λ) An expression for Θ can be written down. The probability that a link of our node points to a node with s links is sP (s) (from the preferential attachment rule), so that: X kP (k)ρk P Θ(λ) = . (65) s sP (s) k Equations (??) and (??) enables the calculation of both ρk and Θ. Then the average density of infection is calculated as X ρ= P (k)ρ(k), (66) k 29 using the value of P (k) for the scale-free network (P (k) = 2m2 /k −3 ). Finally, this gives ρ ≃ 2 e−1/mλ (67) which is the observed behavior, showing the lack of an epidemic threshold, i.e. λc = 0. (It is said that for scale-free networks with γ > 4 a finite threshold is recovered, γ being the power in the degree distribution P (k) ∼ k −γ .) These are very bad news. As has been pointed [?, ?], the network of human sexual contacts seems to be scale-free with an exponent between 2 and 3 (from a survey of several thousand people in Sweden). Both facts together could imply that reducing the transmission probability of a sexually transmitted disease might not eradicate it. In particular, lacking a vaccine for AIDS, the HIV will eventually reach all the exposed network. Both works offer also some optimistic news: targeted immunization can restore a threshold. It is clear that treating the few nodes with extraordinary large connectivity that these networks have is the best strategy. These, unfortunately, are almost impossible to find in a real-life situation. However, it is shown that any departure from uniform immunization, infinitesimal as it were, would restore an epidemic threshold λc 6= 0. References [1] Murray, J. D., Mathematical biology (Springer, Berlin, 1993). [2] H. W. Hethcote, The mathematics of infectious diseases, SIAM Review 42, 599 (2000). [3] F. C. Hoppensteadt, Mathematical methods in population biology (Cambridge University Press, Cambridge, 1982). [4] R. M. Anderson (ed.), Population dynamics of infectious diseases: Theory and applications (Chapman and Hall, London, 1982). [5] D. J. Watts, Small Worlds (Princeton University Press, Princeton, 1999). [6] R. Albert and A.-L. Barabási, Statistical mechanics of complex networks, cond-mat/0106096 (2001). [7] M. E. J. Newman, Models of the small world: mat/0001118 v2 (2000). 30 A review, cond- [8] M. Kuperman and G. Abramson, Small world effect in an epidemiological model, Phys. Rev. Lett. 86, 2909 (2001). [9] R. Pastor-Satorras and A. Vespignani, Optimal immunisation of complex networks, cond-mat/0107066 (2001). [10] D. H. Zanette and M. Kuperman, Effects of immunization in smallworld epidemics, cond-mat/0109273 (2001). [11] D. H. Zanette, Dynamics of rumor propagation on small-world networks, cond-mat/0110324 (2001). D. H. Zanette, Critical behavior of propagation on small-world networks, cond-mat/0105596 (2001). [12] R. Pastor-Satorras and A. Vespignani, Epidemic spreading in scale-free networks, Phys, Rev. Lett. 86, 3200 (2001). [13] Z. Dezsö and A.-L. Barabási, Can we stop the AIDS epidemic?, condmat/0107420 (2001). 31