Influence of Network Heterogeneity on Chaotic Dynamics of Infectious Diseases

Influence of Network Heterogeneity on Chaotic Dynamics of Infectious Diseases Carlo Piccardi and Renato Casagrandi Politecnico di Milano, Dipartimento di Elettronica e Informazione, Milano, Italy (e-mail: [piccardi,casagran]@elet.polimi.it). Abstract: It is well known that the topology of contacts among individuals can sharply influence the persistence of a disease. Much less is known on how the network structure influences the dynamical properties in the important case of seasonal diseases. Aim of this work is to study a periodically forced SIR contact process with demography in host populations modeled through complex networks (of either Erd˝os-R´enyi or scale-free type). We systematically perform oneparameter bifurcation analyses with respect to the strength of seasonality, and find that the epidemiological regime is largely independent of the network class. However, the structure of the host networks does matter. In fact, the heterogeneity of the network degree distribution emerges as a key element in determining the epidemic temporal patterns. In particular, we find that the dynamical complexity is maximal (i.e., chaotic) at intermediate values of the heterogeneity. Keywords: Chaos, networks, bifurcations, epidemic models, Lyapunov exponent. 1. INTRODUCTION With no doubts, the discipline of science where the complex network concept has most favorably spread is epidemiology. It is in fact quite spontaneous and does not require much empirical evidence to recognize that individuals cannot be treated “as average” in terms of pathogen transmissions. For example, in sexually transmitted or in childhood diseases – just to limit our focus on human populations – the social behavior of each individual varies her/his risk of infection, thus it can enhance or reduce its role of potential spreader in the population. This is why the extensive use of ODE approaches to epidemics, rooted in the pioneering work by Kermack and McKendrick (1927), is currently facing a deep revision in the light of the emerging complex network paradigm. Particularly studied in this new context is the problem of disease persistence (Pastor-Satorras and Vespignani, 2001a, 2001b), because of its crucial importance for public health policies. However, another undesirable characteristic of epidemiological time-series is their apparently erratic nature. Seasonality is a key component of many human (Grassly and Fraser, 2006) and animal infectious diseases (Altizer et al., 2006) and it is one of the many possible causes of irregular outbreaks observed in field data (see for example Earn et al., 2000). We do not enter the debate about the extent at which the observed dynamics can be sufficiently well explained by deterministic models evolving on a strange attractor. Some researchers in fact propose an interplay of such models with stochasticity (Rohani et al., 2002) while others opted for transient dynamics Bauch and Earn (2003). Rather, here we focus on the effects played by the network structure on the emergence of chaos in simple, yet nonlinear, epidemiological models. We both have analyzed via bifurcation analysis how seasonality can affect the

dynamics of infectious diseases in the SIR (Kuznetsov and Piccardi, 1994) and the SIRC model (Casagrandi et al., 2006), both models being based on the hypotheses of wellmixed host populations. In the present study we concentrate on a single contact process – the SIR model with demography – over three different network topologies. 2. EPIDEMIC DYNAMICS ON NETWORKS We model the population of individuals as a network composed of N nodes. The i-th node has degree ki , i.e., it is connected by 0 < kmin ≤ ki ≤ kmax links to other nodes. The network is characterized by its degree distribution pk (k = kmin , . . . , kmax ), where 0 ≤ pk ≤ 1 is the fraction of nodes having exactly degree k. We denote by µ and σ, respectively, the mean value and the standard deviation of the degree distribution. Thus µ is the average degree of the network. The epidemiological problem studied here is the simple and traditional, yet seasonally forced, SIR contact process (Anderson and May, 1992). At any time instant t, each node i is in one of the three possible states: Susceptible (or infectable), Infective (and infected), or Recovered (i.e., permanently immune). During a short time interval ∆, an infective node can recover with probability γ∆ while a susceptible node (say i) can become infective with probability ρ∆ni , where ni is the number of infectives among its neighbors. It is particularly important to notice that we also account for hosts’ demography. In fact, we assume that each node can die with probability η∆, irrespective of its degree and its current state, while it can give birth to another node with the same probability (η∆). The newborn nodes are attached to the network with rules that are detailed below, but always in a susceptible state. The birth and death rates are assumed to be equal, in

order to simulate populations that exhibit only stochastic fluctuations around the constant value N . As will be discussed in the next section, although the birth-death process underlying the epidemics does not drastically alter the initial number of nodes in the network, it has the power of strongly modifying its degree distribution. To obtain a model that describes the disease dynamics over the entire network, we follow the approach by PastorSatorras and Vespignani (2001a, 2001b) and subdivide the nodes into classes that depend on their current state and their degree. We therefore introduce 2d state variables Sk (t) and Ik (t) (where d = kmax − kmin + 1) that describe the fraction of nodes with degree k that are, respectively, in susceptible and infective state at time t. Note that the fraction of nodes in recovered state are simply Rk (t) = 1 − Sk (t) − Ik (t). Letting ∆ → 0, we derive the following set of 2d differential equations: ˜ S˙ k = η − ηSk − βSk (k/µ)I, (1) ˙ ˜ Ik = βSk (k/µ)I − (η + γ)Ik , where β = ρµ is the so-called transmission rate, while I˜ =

kX max k=kmin

kpk Ik µ

(2)

is the expected proportion of infectives at time t among the neighbors of a node, because qk = kpk/µ is the degree distribution of the neighbors (e.g., Boccaletti et al., 2006). The global disease prevalence at time t is therefore given by kX max I(t) = pk Ik (t). (3) k=kmin

In the simple case of strictly homogeneous network (pµ = 1, pk = 0 for all k 6= µ), it is straightforward to verify that model (1) reduces to the standard SIR model: S˙ = η − ηS − βSI, (4) I˙ = βSI − (η + γ)I, with R(t) = 1 − S(t) − I(t). Provided that β is not too small (β > η+γ), system (4) has a unique, strictly positive equilibrium which is globally asymptotically stable. The long-term behavior of the epidemic is therefore trivially constant. To account for seasonality, we adopt the typical hypothesis formulated in mathematical biosciences (Dietz, 1976) that the transmission rate β = β(t) of the disease is periodically varying through time, with period equal to 1 year: β(t) = β0 [1 + ε cos(2πt)] .

(5)

The two coefficients β0 > 0 and 0 ≤ ε ≤ 1 are the the so-called average transmission rate and strength of seasonality. The properties of epidemiological ODE models with seasonal transmission rates of type (5) have been extensively studied in the past, particularly the SIR model (Smith, 1983; Aron and Schwartz, 1984; Kuznetsov and Piccardi, 1994), but also variations (e.g., Casagrandi et al., 2006). Bifurcation analyses have shown that, by varying β0 and ε, periodic solutions with periods multiple of 1 year, as well as chaotic behavior, can easily be obtained. Both regular and irregular oscillations are suited to explain evidences

emerging from large sets of available data (Schaffer and Kot, 1985; Olsen et al., 1988). The basic scenario giving rise to the emergence of chaos in (4)-(5) is the cascade of period-doubling bifurcations. By keeping β0 fixed, a stable periodic solution with period of 1 year is obtained for sufficiently small ε. By increasing ε, a sequence of bifurcations leads to cycles with period 2, 4, etc., until chaos is reached at a threshold εchaos . Although different phenomena can be detected at different values of β0 , the above scenario is the most common in a range of values for β0 which is reasonable for a number of diseases (Olsen et al., 1988). 3. METHOD OF ANALYSIS Aim of this work is to study how the appearance of chaos – more technically, the minimum strength of seasonality εchaos leading to chaos – depend on the heterogeneity of network topology. To this end, we therefore need to (i) introduce the different network topologies we work with, (ii) define indices of heterogeneity that can serve to synthesize our results, and (iii) identify a method for inspecting εchaos for various levels of network heterogeneity. 3.1 Networks topologies In the present study we deal with three different networks topologies. Two of them are the most widely used types of networks in epidemiology, namely (a) the Erd˝ os-R´enyi and (b) the scale-free networks. The third network type (c), evolved scale-free network, is uncommon only because the discovery that scale-free networks lose their structure while evolving under birth-death processes is very recent (Moore et al., 2006). The degree distribution of such networks and the method used to construct them are as follows: (a) The Erd˝ os-R´enyi (ER) network (1959) is obtained by randomly connecting N nodes with a pre-specified number of links. For this reason, such networks are often, yet improperly, denoted as random networks. The degree distribution of an ER network with large N ’s and average degree µ is given by a Poisson distribution (e.g., Newman, 2003): e−µ µk . (6) k! It can be shown that, under suitable assumptions, degree distributions of type (6) remain unaltered if the population evolves through a birth-death process. A formal proof is available in Moore et al. (2006). Here it suffices to notice that the number of links is preserved. On the one hand, in fact, if a randomly selected node (say i) dies, all its ki links are removed thus, on average, η∆µ links are lost per unit time. On the other hand, however, if the degree of each newborn node is extracted by the same Poisson distribution (6) which characterizes the original network prior the birth event, then on average η∆µ links are added per unit time. (b) Scale-free networks (SFn) are highly heterogeneous networks because they have very few nodes (the hubs) linked to many others, but a large number of nodes with scarce connections. In recent years, SFn have pk =

Technically speaking, we obtain the degree distributions pk ’s that are used with model (1) as follows. First, we fix the values for the minimum and maximum degree of the networks under study (kmin = 1 and kmax = 100). The necessity of truncating the right tail of the degree distribution is due to computational costs (remember that our model has 2[kmax − kmin + 1] equations). Building the degree distribution of ER networks is simple, because the pk ’s can directly be computed via formula (6). The only parameter that must be specified is the average degree µ. As already pointed out above, in fact, ERn are unaltered by birth-death processes, thus no further calculus is needed.

degree distribution pk

0.1

(c) 0.01

0.001

(a)

(b)

0.0001 1

10

100

node degree k

Fig. 1. Plots in log-log scales of the degree distributions considered in this study: (a) Erd˝os-R´enyi, (b) scalefree with scale parameter c = 3, and (c) evolved scale-free networks. In all cases, the average degree is µ = 10. received great attention because they have been discovered in a wide variety of human and artificial contexts (Barab´asi and Albert, 1999; Newman, 2003; Boccaletti et al., 2006), including epidemiology (e.g., Liljeros et al., 2001). The degree distributions of SFn are power-law functions that, at least for large k’s, take the form: pk ≈ k −c , (7) where c > 0 is the scale parameter. (c) In contrast to results obtained with ERn, the degree distribution of SFn is strongly modified by the birthdeath process. Starting from a network structure of type (b), if we assume that death selects nodes at random, on average η∆µ links are cancelled per unit time. This loss can easily be compensated by designing the birth process in such a way that exactly µ links are attached with each newborn node, thus η∆µ new links per unit time. To preserve the scale-free nature of the network, one might think to connect newborn nodes with the so-called preferential attachment paradigm (Barab´asi and Albert, 1999), which is known to create degree distributions of the form (7). The key point is that, even under such operational conditions, the power-law distribution is lost under the demographic dynamics in favor of a new distribution called “stretched exponential” (Moore et al., 2006) that, for large k’s, takes the form pk ≈ k −3/4 e−2

√ k

.

(8)

A network resulting from the process just described is denoted in this paper as evolved scale-free network (ESF). Some exemplificative degree distributions of ER, SF and ESF networks are reported in Fig. 1 for a direct comparison. Notably, in contrast to what happens with the power-law distribution in a suitable range of values for parameter c, the second moment hk 2 i of distributions (8) remains finite even for kmax → ∞. This means that, under the birth-death process, the degree distribution loses the “fat-tail” feature that makes SF so important.

On the contrary, when the epidemic process is studied on a SFn, the dynamics of the network structure per se must be accounted for. To construct the initial degree distribution, we use the Barab´asi-Albert algorithm (1999), which yields a degree distribution (7) with c = 3, with an a priori specified average degree µ. The evolution of the pk ’s under the birth-death process is computed according to the following discrete-time equation, proposed by Moore et al. (2006): N pk (t + ∆) = N pk (t)+ +η∆ (−pk (t) − kpk (t) − µπk pk (t)) + (9) +η∆ ((k + 1)pk+1 (t)) + +η∆ (µπk−1 pk−1 (t)) + η∆δkµ , where δkµ equals 1 if k = µ or 0 otherwise, and πk = k/µ is the preferential attachment kernel, i.e. it is proportional to the probability that a newborn node attaches one of its links to an already existing node of degree k. Each of the other terms in (9) describe one of the possible mechanisms that can alter the generic degree k of a selected node – thus changing the entire degree distribution of the network. Such mechanisms include the death of the degree-k node, its birth, and the birth or death of one of its neighbors. These last two events increase or decrease of one unit the degree k of the existing node. Starting from a power-law of form (7), we let the network evolve according to (9) until we obtain convergence toward the invariant degree distribution, which turns out to be of type (8). 3.2 Measures of heterogeneity From an intuitive viewpoint, the more diversified is the degree distribution of a network, the higher must be its measure of heterogeneity. Figure 2 shows two simplistic examples of network characterized by very different degree distributions. The top panel is a sketch of a strictly homogeneous network (spiky degree distribution), while the bottom panel represents a highly diversified network where each possible degree is held by the same fraction of nodes (flat degree distribution). No index is required in this case to recognize that the bottom network is more heterogeneous than the top one. However, as soon as the examples are not taken at the extremes of the spectrum (spiky vs. flat degree distributions), some quantitative index is needed to confer an appropriate measure of heterogeneity to a network. Several indicators can be used to this end like, for example, the coefficient of variationPσ/µ or the Shannon entropy of the degree distribution k pk log 1/pk . In this work we prefer to use an alternative approach which is based on

leading Lyapunov exponent Λ of the attractor. This last quantity is computed to precisely detect the threshold of chaos εchaos . An important detail is that, except for the value ε0 for which they are irrelevant, the initial values for the state variables of the z-th simulation are the final states obtained while stopping the (z − 1)-th simulation. Not only this procedure usually reduces the transient time toward the attractor, but it lowers the probability of reaching a different attractor in presence of multistable dynamics, a behavior that is known to exist in periodically forced SIR-like models (Kuznetsov and Piccardi, 1994; Casagrandi et al., 2006).

Fig. 2. Two examples of network with N = 14 and µ = 4. Top: a strictly homogeneous network (p4 = 1, pi = 0 for all i 6= 4) with σ/µ = 0, H = 0. Bottom: a nonhomogeneous network (pi = 1/7 for all 1 ≤ i ≤ 7) with σ/µ = 0.519, H = 0.286 (the darkness of a node is proportional to its degree). the Gini index, widely used in microeconomics to quantify income inequality (Colander, 2001). The use of such an index has been recently extended to complex networks and proved to be particularly effective in this context (Hu and Wang, 2008). Network heterogeneity is quantified by an index H that measures the average “degree inequality” in the network, i.e.: PN PN i=1 j=1 |ki − kj | , (10) H= 2N 2 µ where ki is the degree of node i (an alternative expression directly relating H to the degree distribution pk can be found in Hu and Wang (2008)). Larger values of H imply higher levels of heterogeneity. It can be ascertained that H = 0 for a strictly homogeneous network and H → 1 for a completely heterogeneous network, although that limit can only be approached by infinite networks, i.e. for N → ∞. 3.3 Inspecting chaos in different networks For each network type and parameterization, we perform a one-parameter bifurcation analysis with respect to the strength of seasonality ε – see equation (5) above. All other biological parameters of the SIR model are kept constant and fixed at values consistent with those suggested by the epidemiological literature (Olsen et al., 1988), i.e., η = 0.02, γ = 100, β0 = 500. In particular, the bifurcation scenario is obtained by simulating system (1),(5) for increasing values [ε0 , ε1 , . . . , εz , . . . , εmax ] of the strength of seasonality starting from ε0 = 0, and by monitoring the

The numerical analysis outlined above is rather difficult for two reasons. First it can be verified that, for sufficiently large – still realistic – values of ε, model (1) displays chaotic dynamics. Its attractors are such that the variables Ik ’s undergo very broad oscillations that involve many orders of magnitude and become extremely close to zero for quite long time intervals. A more accurate numerical integration of the model can be obtained by replacing equations (1) with their semi-logarithmic counterpart with state variables Lk = log(Ik ), i.e. ˜ S˙ k = η − ηSk − β(t)Sk (k/µ)I, ˜ (11) β(t)Sk (k/µ)I − (η + γ), L˙ k = exp(Lk ) Pkmax where I˜ = k=k q k exp(Lk ). A second difficulty is that min the computation of the Lyapunov exponent Λ requires the integration, in parallel to (11), of a variational equation of the form · ¸ · ¸ v˙ Sk v Sk =J , (12) v˙ Lk vLk where J is a large, 2d × 2d Jacobian matrix having the following structure   ∂ S˙ k ∂ S˙ k  ∂S ∂L   k k  . J = (13)   ˙ ˙  ∂ Lk ∂ Lk  ∂Sk ∂Lk The entries of J are evaluated along the solution of (11). Since the accurate computation of Λ requires very small and regular integration steps (Ramasubramanian and Sriram, 2000), we found that Euler integration is the best strategy for the numerical analysis of equations (11) and (12). 4. RESULTS Figure 3 shows the typical bifurcation scenario obtained by simulating model (11) under sinusoidal forcing (5) while assigning different values to the strength of seasonality ε. The adjective “typical” used above is justified because the same scenario emerges for a wide range of values attributed to all parameters other than ε and, specifically, for both the ER and ESF network topologies described in the previous section. The top panel of Fig. 3 displays the route to chaos (essentially a Feigenbaum cascade) that has been found while studying the contact process over an Erd˝os-R´enyi network. The bottom panel shows the leading Lyapunov exponent of the corresponding attractor. Of much interest to the present study is the understanding of how the network heterogeneity influences the appear-

chaos threshold Hchaos

log prevalence L

0 -10 -20 -30 -40 -50

0.08

0.07

0.06

ESF

0.05

0.04 0.1

0.2

-60

Lyap. exp. /

ER

0.3

0.4

0.5

heterogeneity H

0.4 0.2 0.0

0.00

0.02

0.04

0.06

0.08

Fig. 4. The chaos threshold εchaos as a function of the heterogeneity H in ER and ESF networks. The dashed line marks the chaos threshold in a strictly homogeneous network (H = 0).

strength of seasonality H

ance of chaos. The system tendency to behave irregularly can be captured by the focal parameter εchaos . In particular, lower values of the bifurcation threshold εchaos indicate more possibilities of obtaining less predictable epidemic outbreaks, both in terms of time of occurrence and fraction of the population involved. By systematically performing analyses such as those reported in Fig. 3 for networks of different types (i.e., ER and ESF) and characterized by different levels of heterogeneity, we calculate εchaos (H) and obtain the patterns shown in Fig. 4. As a reference, the dashed line in the same figure represents the bifurcation strength of seasonality εchaos for a strictly homogeneous network. Interestingly, we first notice that the dependence of εchaos on H is rather independent of the network type, not only from a qualitative but also from a quantitative viewpoint. In a sense, this is not surprising because the degree distributions of the ER and ESF networks (see again Fig. 1) are not very different. We remind the reader that the scalefree networks are transformed into ESF network under the pressure of birth-death processes, therefore it would not make any sense to ask whether the similarity of the two patterns in Fig. 4 would be preserved in SFn. Second, all the values of εchaos computed for heterogeneous networks (H > 0) are placed below the dashed line in Fig. 4. This could be effectively synthesized by saying that the heterogeneity of the network enhances chaotic dynamics. However, such a statement does not entirely capture our findings. In fact, third and more important, εchaos (H) is an U-shaped function that reaches its minimum for intermediate values of H. There is therefore an interval of values for the strength of seasonality ε where chaotic epidemics are obtained only in networks that are neither too regular (homogeneous) nor too diversified (heterogeneous). The epidemiological consequences of this result are straightforward. The strength of seasonality ε is a parameter that increases with latitude, at least for the many pathogens transmitted by droplets, aerosols or fomites that must

log prevalence L

Fig. 3. The bifurcation diagram of the logarithmic prevalence L = log(I) (top) and the Lyapunov exponent Λ (bottom) as functions of the strength of seasonality ε for an Erd˝os-R´enyi network with H = 0.43, σ/µ = 0.80.

0 -5 -10 -15 -20 -25 -30 0.1

0.2

0.3

0.4

0.5

heterogeneity H

Fig. 5. The bifurcation diagram of the logarithmic prevalence L = log(I) as function of the heterogeneity H of the network, for a fixed value ε = 0.045 of the strength of seasonality (ER network). therefore survive outside the host. In fact, variations of environmental conditions such as temperature, humidity or exposure to sunlight are much lower in the tropics than toward poles. In such cases, a bifurcation study performed with a fixed value of ε but for changing values of H has the biological meaning of discussing the spread of the disease in populations (e.g., cities) with different degree distributions but located at similar latitude or with similar environmental variations. The outcomes predicted by our model in such a virtual experiments for a given value of ε are reported in Fig. 5. While increasing H we first have a direct period-doubling cascade to chaos, then a reverse cascade – resulting in a ‘bubble’ of chaos at intermediate levels of heterogeneity. To what extent is this scenario realistic? Validation is quite a chimera in these sort of problems. However, if we assume that larger cities have ceteris paribus more diversified degree distributions than smaller cities, then some comparison with historical data is possible. For example, while contrasting two decennial datasets of measles and chickenpox in a relatively small population (Bornholm, Denmark, 5 · 104 inhabitants) and in the most populous city of the US (New York City, 8 · 106 inhabitants in 2007), Olsen and Schaffer (1990) found that there is a pronounced difference between the temporal patterns of the monthly notifications of infectives in large and small cities. In particular, the power spectra shown in the cited work as well as in previous ones (Kot et al., 1988; Olsen et al., 1988) reveal that the overall patterns of infectives,

although irregular, are less complex in very large cities than in smaller ones. A portrait which is compatible with the results presented here. Further and promising directions of research include, therefore, a more quantitative exploration of empirical datasets for those diseases whose dynamics can be described by the SIR model. REFERENCES Altizer, S., Dobson, A., Hosseini, P., Hudson, P., Pascual, M., and Rohani, P. (2006). Seasonality and the dynamics of infectious diseases. Ecol. Lett., 9(4), 467–484. Anderson, R. and May, R. (1992). Infectious Diseases of Humans. Oxford University Press, Oxford, UK. Aron, J. and Schwartz, I. (1984). Seasonality and perioddoubling bifurcations in an epidemic model. J. Theor. Biol., 110(44), 665–679. Barab´asi, A. and Albert, R. (1999). Emergence of scaling in random networks. Science, 286(5439), 509–512. Bauch, C. and Earn, D. (2003). Transients and attractors in epidemics. Proc. R. Soc. B - Biol. Sci., 270(1524), 1573–1578. Boccaletti, S., Latora, V., Moreno, Y., Chavez, M., and Hwang, D. (2006). Complex networks: Structure and dynamics. Phys. Rep., 424(4-5), 175–308. Casagrandi, R., Bolzoni, L., Levin, S., and Andreasen, V. (2006). The SIRC model and influenza A. Math. Biosci., 200(2), 152–169. Colander, D. (2001). Microeconomics. McGraw-Hill, Boston. Dietz, K. (1976). The incidence of infectious diseases under the influence of seasonal fluctuations. In J. Berger, W. B¨ uhler, R. Repges, and P. Tautu (eds.), Lecture Notes in Biomathematics, volume 11, 1–15. Springer. Earn, D., Rohani, P., Bolker, B., and Grenfell, B. (2000). A simple model for complex dynamical transitions in epidemics. Science, 287(5453), 667–670. Erd˝os, P. and R´enyi, A. (1959). On random graphs. Publ. Math.-Debr., 6, 290–297. Grassly, N. and Fraser, C. (2006). Seasonal infectious disease epidemiology. Proc. R. Soc. B - Biol. Sci., 273(1600), 2541–2550. Hu, H. and Wang, X. (2008). Unified index to quantifying heterogeneity of complex networks. Physica A, 387(14), 3769–3780. Kermack, W. and McKendrick, A. (1927). A contribution to the mathematical theory of epidemics. Proc. R. Soc. A, 115(772), 700–721. Kot, M., Schaffer, W., Truty, G., Graser, D., and Olsen, L. (1988). Changing criteria for imposing order. Ecol. Model., 43(1-2), 75–110. Kuznetsov, Y. and Piccardi, C. (1994). Bifurcation analysis of periodic SEIR and SIR epidemic models. J. Math. Biol., 32(2), 109–121. Liljeros, F., Edling, C., Amaral, L., Stanley, H., and Aberg, Y. (2001). The web of human sexual contacts. Nature, 411(6840), 907–908. Moore, C., Ghoshal, G., and Newman, M. (2006). Exact solutions for models of evolving networks with addition and deletion of nodes. Phys. Rev. E, 74(3), 036121. Newman, M. (2003). The structure and function of complex networks. SIAM Rev., 45(2), 167–256. Olsen, L. and Schaffer, W. (1990). Chaos versus noisy periodicity - Alternative hypotheses for childhood epi-

demics. Science, 249(4968), 499–504. Olsen, L., Truty, G., and Schaffer, W. (1988). Oscillations and chaos in epidemics: a nonlinear dynamic study of six childhood diseases in Copenhagen, Denmark. Theor. Popul. Biol., 33(3), 344–370. Pastor-Satorras, R. and Vespignani, A. (2001a). Epidemic dynamics and endemic states in complex networks. Phys. Rev. E, 63(6), 066117. Pastor-Satorras, R. and Vespignani, A. (2001b). Epidemic spreading in scale-free networks. Phys. Rev. Lett., 86(14), 3200–3203. Ramasubramanian, K. and Sriram, M. (2000). A comparative study of computation of Lyapunov spectra with different algorithms. Physica D, 139(1-2), 72–86. Rohani, P., Keeling, M., and Grenfell, B. (2002). The interplay between determinism and stochasticity in childhood diseases. Am. Nat., 159(5), 469–481. Schaffer, W. and Kot, M. (1985). Nearly one dimensional dynamics in an epidemic. J. Theor. Biol., 112(2), 403– 427. Smith, H. (1983). Subharmonic bifurcation in an S-I-R epidemic model. J. Math. Biol., 17(2), 163–177.