Robust oscillations in SIS epidemics on adaptive networks: Coarse-graining by automated moment closure

0.33 MB·English

by  Thilo Gross

epl draft
Robust oscillations in SIS epidemics on adaptive networks:
Coarse-graining by automated moment closure
Thilo Gross12 and Ioannis G. Kevrekidis3
1 Dept. of Chemical Engineering, Princeton University, Princeton NJ 08544, USA.
8 2 Max-Planck Institute fu¨r Physik komplexer Systeme, No¨thnitzer Straße 38, 01187 Dresden, Germany. email:
0
[email protected]
0
2 3 Dept. of Chemical Engineering and PACM, Princeton University, Princeton NJ 08544, USA. email: yan-
[email protected]
n
a
J

4

] PACS 89.75.Hc–Networks and genealogical trees
O PACS 87.19.Xx–Diseases
A PACS 89.75.Fb–Structuresand organization in complex systems
n. Abstract. - We investigate the dynamics of an epidemiological susceptible-infected-susceptible
i (SIS)model on an adaptivenetwork. This model combinesepidemic spreading (dynamicsonthe
l
n network)withrewiringofnetworkconnections(topologicalevolutionofthenetwork). Wepropose
[ and implement a computational approach that enables us to study the dynamics of the network

directlyonanemergent,coarse-grainedlevel. Theapproachsidestepsthederivationofclosedlow-
2
dimensionalapproximations. Ourinvestigationsrevealthatglobalcoupling,whichentersthrough
v
theawarenessofthepopulationtothedisease,canresultinrobustlarge-amplitudeoscillations of
7
4 thestate and topology of thenetwork.
0
2
0
7
0
/ Introduction.– Overrecentyearsthephysicsofnet- olddiseases [17]as wellas the emergence of new ones [18]
n
works has received much attention [1–3]. Research has atanacceleratingpace. Conceptualmathematicalmodels
i
nl focused on the one hand on dynamics of networks, i.e., canprovidefundamentalinsightsontheunderlyingmech-
: networkgrowthandrewiring[4,5],andonthe other hand anisms that govern epidemic dynamics.
v
on dynamics on networks; the coupling of individual dy- One of the most simple conceptualmodels ofepidemics
i
X namical systems according to a static network topology isthesusceptible-infected-susceptible(SIS)model[19],de-
r [6,7]. Only recently have these two distinct strands of scribed in more detail below. If considered on a static
a
network research been brought together in the study of network, varying the parameters in this model reveals at
adaptive networks,whichcombinetopologicalevolutionof mostonedynamicaltransition. Thistransitioniscontinu-
the network with dynamics on the network nodes [8–14]. ous and corresponds to the epidemic threshold: the point
Adaptive networks have been shown to exhibit a num- in parameter space beyond which the disease can invade
ber of new phenomena, including robust self-organization the network. If one enables the individuals in the model
towards dynamical criticality [8], the formation of com- to avoid contact with infected ones by altering their lo-
plex global topologies based on simple local rules [13,15], calinteractiontopologythesystemturnsintoanadaptive
the emergence of new bifurcations and phase transitions network. It has recently been shown that this natural ex-
involving local as well as topological degrees of freedom tensionofthespatialSISmodelgivesrisetomorecomplex
[12,13], and finally a spontaneous ‘division of labor’ in dynamics, including bistability, hysteresis and a narrow
which an initially homogeneous population of nodes self- region of oscillatory dynamics [13].
organizes into functionally distinct classes [10]. In order to understand the subtle interplay between
A prominent application of networks is epidemiology. local dynamics and topological evolution that character-
In an ever more populated and tighter connected world izes adaptive networks the application of tools from dy-
epidemicdiseasesareonceagainbecomingamajorthreat namical systems theory is desirable. While these meth-
[16]. In the past decades we have witnessed the return of ods are in general only applicable to relatively low-
p-1
T. Gross and I.G. Kevrekidis
dimensional equation-based models, previous results [13, Coarse graining. – It has been shown in [13] that
20] suggest that often a small number of topological de- the system-level dynamics of sufficiently large networks
grees of freedom suffices to characterize the dynamics of (here N = 105,L = 106) can be captured by the number
the individual-based model on an emergent level. This of SS-links lSS, the number of II-links lII and the number
means that a small number of observables can describe ofinfectedindividualsi. Thenumberofsusceptiblessand
the macroscopic (system-level) state of the network and the number of SI-links lSI then follow from the conserva-
thusalow-dimensionaldescriptionatthislevelis feasible. tion laws s+i = N and lSS +lSI +lII = L. A moment
In [13] a low-dimensional set of system-level equations of expansionrevealsthatthedynamicsofthecollectivestate
motion for the adaptive SIS model was derived by means variables s,lSS,lII can be described by
of a moment expansion and subsequent moment closure
approximation. While this approach has yielded satisfac- ddti = plSI−ri
tory results, it involves a strong homogeneity assumption ddtlII = plSI(cid:16)llISSII +1(cid:17)−ri2liII (1)
and therefore can fail in certain topologies. ddtlSS = rilSiI −plSIllSSSII +wlSI
Inthis Letterweextendthepreviousworkintwoways.
First,weconsidertheeffectofavariablelevelofawareness where lISI and lSSI denote the number of I-S-I and S-S-I
to the disease. This awareness is shared instantaneously chains in the network, respectively.
among the population and therefore acts as a global cou- Note that we have so far only assumed that the net-
pling. Second, we propose and implement an alterna- work is large, so that the state can be described by con-
tive computationalapproach,basedonthe coarse-grained tinuous variables and the effect of individual-based (de-
‘equation-free’ modeling and analysis framework [21]. In mographic) stochasticity can be neglected. In particular
this approach the analytical moment-closure approxima- no strong homogeneity has been used since many effects
tionis replacedby anumericalprocedurewhichautomat- of non-homogeneous topology are captured by taking the
ically extracts the effects of appropriate closure. This ap- densities of links lSS,lSI,lII and chains lISI,lSSI into ac-
proachenablesustoinvestigatetheemergentdynamicsof count. However, Eq. (1) does not yet constitute a closed
complexnetworkswithwell-establishedtoolsofdynamical model, as we haven’tspecified how lISI andlSSI are deter-
systemstheorywithoutexplicitlyderivingaclosedanalyti- mined.
caldescriptiononthislevel. Ourinvestigationrevealsthat In Ref. [13] the moment-closure approximation lISI =
the globalcouplingstronglyextends the parameterregion lS2I/s,lSSI = lSSlSI/2s was used to close the model. This
inwhichoscillatorydynamicscanbe observed. Itthereby approximation introduces a homogeneity assumption as
gives rise to persistent large-amplitude oscillations affect- it requires that the SI- and SS-links are homogeneously
ingthe prevalenceofthediseaseaswellasthetopologyof distributed among the possible three node chains. While
the network. the moment closure approximation has yielded good re-
sultsin[13]ithastofailinnetworkswithslowlydecaying
The adaptive SIS model. – We consideranetwork degree distribution, such as scale free graphs.
with a fixed number of nodes N and undirected links L. In this Letter we avoidthe moment closure approxima-
Everynode representsanindividual,while links represent tion: We do not attempt to derive analytical expressions
socialcontacts between individuals. An individual can be for lISI and lSSI in order to close Eq. (1). Instead, we
either infected (I) or susceptible (S). We denote the links specify a numerical procedure that generates the effect of
betweenindividualsasSS-links,SI-linksorII-linksaccord- appropriateclosuretermson-demand fromshortburstsof
ing to the states of the nodes that they connect. In every numericalsimulation. Thereforethe majortechnicalchal-
timestepinfectedindividualsrecoverwithaprobabilityr, lengeinthisLetteristodeterminethecorrectvalueoflISI
becoming susceptible. For every SI-Link there is a proba- and lSSI for a given set of state variables i,lSS,lII.
bility p that the susceptible individual becomes infected. In general one can image a large number of candidate
InadditiontothesestandardrulesofSISmodels,weallow network topologies that agree with the desired values of
susceptible individuals to avoid contact with infected in- i,lSS,lII. Thus, for a given set of i,lSS,lII, many different
dividuals by altering their local interaction topology: For values of lISI and lSSI are in principle possible. However,
everySI-Linkthereisaprobabilitywthatarewiringevent mostof the candidate topologies will not arise in the long
occursinagiventimestep. Inarewiringevent,thesuscep- term dynamics of the network. In fact one finds that the
tible individualcutsits linkto theinfectedindividualand lISI and lSSI that are consistent with given i,lSS,lII are
establishesanewlinktoanothersusceptibleindividual. In restricted to a small set in the long-term dynamics.
contrastto [13]we do not fix the rate w, but assume that Consider the following: Large adaptive networks can
it changes according to the individual’s awareness of the easilyhavemillionsoftopologicaldegreesoffreedom. The
disease. In human populations with accessto mass media assumption that the state of such a network can be char-
it is reasonable to assume that information on the dis- acterized by a low number of variables implies in general
ease is instantaneously disseminated across the network. that there is a time scale separation between a few slow
We therefore set w =w0ρ, where ρ=i/N, is the infected variables and the remaining fast degrees of freedom. The
fractionofthepopulationandw0 isaconstantparameter. system will then approach a slow manifold, on which the
p-2
Coarse Graining and Oscillations in an adptive SIS model
Fig. 2: A slice of the slow manifold (i=97693, lss =884707).
Shown are theclosure terms lssi and lisi computed by thepro-
posed procedure (solid line) and an analytical approximation
(dashed line). The corresponding results from a long simula-
tionrun(circles)agreewellwiththenumericalresult. N =105,
Fig. 1: (Color online) Schematic of the lifting procedure. A
network with initially arbitrary third moments is settled to L=106,w0 =0.006, p=0.0007, r=0.0002, .
the slow manifold (dotted line) on which the consistent val-
ues encountered in the long-term dynamics are located. For
this purpose we run a short simulation (solid line) and then provide us with a first estimate cISI,1, cSSI,1. To obtain
reset the first and second moments while retaining the third a better estimate, we generate a network at the operat-
∗ ∗ ∗
moments (dashed line). Inset: Illustration of coarse projec- ing point i , lII, lSS; with third moments according to
tiveintegration. Thenecessaryinformationforlongprojective cISI,1, cSSI,1 (explained below). The endpoint of a sec-
leaps (dashed) is extracted from short bursts of individual- ond burst of simulation yields a better estimate cISI,2,
based simulation (solid). cSSI,2. This process is repeated ninit times to converge
on the slow manifold (see Fig. 1). This is followed by
naver = ntot −ninit further repetitions. Averaging over
longtermdynamicstakesplace. Theslowmanifoldcanbe
the endpoints of these repetitions we compute the consis-
ptoapraomlogeticearilzdeedgerenetsiroelfyfrbeyedtohme salroewslvaavreidabtloest,haesslaolwl ovtahreir- tent values c∗x = Pnnt=otninitcx,n/naver, for x=ISI, SSI. The
variance of the approximation and therefore the width of
ables. Thus a set ofslowvariablesdetermines the state of
thebandinwhichtheconsistentvaluesarelocatedcanbe
tchaensbyestwemrit.teTnhiassima pfulinecsttiohnatoafltlhoethselorwdevgarreiaesbloefs.freedom estimated by Pnnt=otninit(c∗x−cx,n)2/(naver2−naver).
Animportantstepinthealgorithmoutlinedaboveisthe
In the system consideredhere the time scale separation
generation of a network with specific i,lSI,lSS,lSSI, and
betweentheslowvariablesi,lSS,lIIandfastvariables,such lISI. Starting from a given network the desired number
aslSSIandlISI,isfiniteandthusthereasoninggivenabove of i,lSI and lII can be reached straightforwardly by first
isonlyapproximatelytrue. Nevertheless,thevaluesoflISI flipping the state of nodes and then selectively rewiring
andlSSI thatcanpossiblyappearinthelongtermdynam- links. Then we set lSSI, lISI: We chose a random link and
ics(whileagivensetofi,lSS,lII isobserved)arerestricted additionallyarandomnodethatis notpartofthe link. If
to a relatively narrow band. This band is simply a noisy
atleastoneofthenodesconnectedtothelinkhasthesame
graphofthe slow manifold. If the time scale separationis
stateastheindividualnodeweselected,wecanrewirethe
sufficiently large, the band of possible lISI or lSSI is nar- link in a way that it leaves i,lSI and lII unchanged but
row enough to treat it effectively as a single value. In
can possibly affect lSSI and lISI. If this rewiring brings us
the following we denote such a value of lISI or lSSI that closertothe desiredvaluesitisaccepted,otherwiseanew
is consistent with a given set of i,lSS,lII in the long term node and a new link is chosen randomly. This process is
dynamics simply as a consistent value.
repeateduntil the desiredlSSI andlISI have been reached.
Let us now describe an algorithm that computes the
consistent values of lISI and lSSI for a given set of state Numerical Methods. – The algorithm described
variablevaluesi∗,l∗ ,l∗. Inordertoaccountforthecom- above allows us to generate consistent states for Eq. (1)
SS II
binatorialeffectoftheslowvariables,itisadvantageousto computationally on-demand. This enables us to use
workwiththe normalizedthirdmomentscISI =slISI/lSI2, Eq. (1) to numerically estimate the temporal derivatives
cSSI = lSSIs/(lSSlSI). We start by randomly generating a ofthesystem-levelstatevariablesatagivenpointinstate
network at the desired operating point i∗, l∗, l∗ . The space. Ittherebyprovidesuswithalltheinformationthat
II SS
third moments lISI and lSSI in a random graph will not, is needed to apply many established tools for the compu-
in general,be consistent, that is, they are not onthe slow tational investigation of dynamical systems.
manifold of the system. In order to “settle” the network Let us consider a simple example. A standard task in
to the slow manifold, where the consistent values are lo- theinvestigationofdynamicalsystemsisthecomputation
cated, we proceed as follows: We run a short simulation ofatrajectory,startingfromagiveninitialpointx0. Per-
of j steps. Over the course of the simulation the sys- haps, the most simple computational tool to perform this
tem approaches the (attracting) slow manifold; the val- taskisforwardEulerintegration: Givenapointx onthe
n
ues of the third moments at the end of the simulation trajectory a subsequent point xn+1 is computed by the
p-3
T. Gross and I.G. Kevrekidis
first order Taylor approximation xn+1 = xn +τ f(x)|x , thus only a small speed-up is realized (a factor of ∼ 2).
n
where f(x) is the right-hand side of the equations of mo- However, a significantly bigger speedup can be expected
tion and τ is a fixed increment in time. By repeated ap- in larger networks. In addition there are other advan-
plication of this procedure a sequence of points is pro- tages: In the derivation of Eq. (1) the individual-based
duced that traces the desiredtrajectory. The approachof based stochasticity of the system is to a large extend fil-
coarse-grained projective integration [22] can be used to teredout. Furthermore,thecoarseintegrationschemecan
apply methods like forward Euler integration to systems runbackwardsintime,whichcanforinstancebeusefulin
in which no closed analytic expression for the right-hand approximatingthe boundaries ofbasins of attraction[23].
sideoftheequationsofmotionisavailable. Followingthis But most importantly the projective Euler integration is
approachthederivatives,suchas dx/dt| = f(x)| ,are a proof of principle; many more powerful methods of dy-
xn xn
not computed analytically, but are estimated from prop- namical systems theory can be applied in essentially the
erly initialized bursts of microscopic simulation. same way.
We have performed coarse projective integration of the ThesimpleexampleofforwardEulerintegrationproves
SIS model considered here. We have also used a variant that it is possible to close the model on-demand by us-
scheme that takes some additional information from the ing only short bursts of microscopic simulation. This al-
moment expansion into account: In order to compute the lows us for instance to use fixed-point algorithms such as
temporalderivative of the state variables at a givenpoint Newton’s method. These methods do not only compute
in state space x = (i,lII,lSS), we first use the procedure stationary states more efficiently than integration to sta-
described in the previous section to compute appropri- tionarity; they can also be used to find dynamically un-
ate lISI and lSSI for Eq. (1). We can then close Eq. (1) stable steady states, such as saddles, that are inaccessible
and use it to estimate the temporal derivatives. Thus, to direct simulation.
in every step of this projective Euler method a new net- Startingfromaninitialestimatex0ofastationarystate,
work is generated in order to compute the derivatives of Newton’s method produces progressivelybetter estimates
the system-levelvariables. However,in orderto construct byxn+1 =xn−J−1f(xn),whereJistheJacobianmatrix
the network and find the desired closure terms only short atx . As the convergenceofourNewtonmethodwasnot
n
burstsofsimulationarenecessary;nolongsimulationruns very sensitive to minor errorsin the Jacobianwe found it
are required. advantageous to estimate f(x) by the automated closure
procedure and approximate the Jacobian by differentiat-
ing the (inaccurate but faster) analytical closure.
The Jacobian is also needed for continuation of a solu-
tion branch and for detection of bifurcations [24]. Since
the eigenvalues can be highly sensitive to errors in the
Jacobianwe approximatethe elements ofJacobianinthis
casebyfinitedifferencesoftemporalderivativescomputed
numerically. Because of the sensitivity of the eigenvalues,
averaging over many (∼ 103) temporal derivative evalu-
ations is necessary to obtain reasonably accurate eigen-
values. A possibly better alternative is the application
of matrix-free iterative algorithms, in which the eigenval-
ues are determined without explicit approximation of the
Jacobian [25].
Fig. 3: (Color Online) Coarse projective Integration. Shown
Results. – By applying the Newton method for sev-
are time series of the infected fraction of the population, the
densityofSS-linkslSS andSI-linkslSI. Thelinknumbershave eralparametervalues,wehavecomputedaone-parameter
been normalized with respect to the total number of links in bifurcationdiagramfortheadaptiveSISnetwork(Fig.4).
the network. The results from a coarse projective integration For p>0.71the networkapproachesa stationarystate in
(black)areingoodagreementwithafullindividual-basedsim- whichtheprevalenceofthediseasesishighρ≈0.97. Ifpis
ulation (orange). Parameters: N = 105, L = 106, p = 0.06, lowered,athresholdiscrossed(pointB)atwhichacoarse
r=0.0002, p=0.0003, ninit =2, naver =1, s=20, τ =100. limit cycle attractor emerges – most likely from a discon-
tinuous non-equilibrium phase transition which appears
A comparison of the coarse projective integration asa foldbifurcationof limit cyclesin the low dimensional
scheme with full individual-based simulation is shown in system. We expect that a coarse saddle cycle is also born
Fig. 3. Coarse projective integration has previously been there. Thestationaryattractorandthelimitcycleattrac-
proposed as a way to speed up lattice Boltzmann, ki- torcoexistuntilthe stabilityofthe stationaryattractoris
netic Monte Carlo as well as molecular dynamics simula- lostinasubcriticalHopfbifurcation(pointA).Insuchbi-
tions[22]. Inthepresentimplementationfornetworks,ini- furcationsanunstablelimitcyclevanishes. Whilewewere
tialization of consistent networks is time consuming, and notable to compute this unstable oscillationdirectly, it is
p-4
Coarse Graining and Oscillations in an adptive SIS model
The dynamics of the adaptive SIS model with variable
awareness reported here differ from the dynamics of the
model with constant rewiring studied in [13]. Rewiring is
in both cases a strong protective mechanism that isolates
infectedindividualsandthusreducestheprevalenceofthe
disease. However, in the present model this mechanism
onlyplaysthe roleifthe prevalenceandthereforealsothe
population’s awareness of the disease is high. In partic-
ular, compared to the non-adaptive SIS model, rewiring
does not increase the epidemic threshold (point F) in the
present model, while a strong increase of the threshold
was found in the model with constant rewiring. Perhaps
the most striking difference between the model variants
is thatthe awarenessmechanismgreatlyincreasesthe pa-
rameterrangeinwhichoscillationscanbeobserved. Also,
Fig. 4: System-level bifurcation diagram. Branches of stable
theamplitude oftheoscillationsisincreasedconsiderably.
steady states (solid lines) and saddle points (dashed/dotted)
Inordertounderstandhowvariableawarenesspromotes
have been computed by a coarse-grained Newton method.
oscillations, the underlying mechanism has to be taken
The white line tentatively marks a branch of saddle limit cy-
into account. While rewiring of network connections iso-
cles. Shaded regions mark ranges of ρ observed during long
individual-based simulations in the neighborhood of the at- lates infected individuals and thus slowly decreases the
tractive limit cycle (light gray) and of stable of stationary so- prevalenceofthedisease,italsoleadstotheformationofa
lutions (dark gray). Computation of the Jacobian eigenvalues highlyconnectedclusterofsusceptibleindividuals(see[13]
revealasubcritical Hopfbifurcation (A),twofold bifurcations fordetails). Ifthediseasemanagestoinvadethisclusterit
(C,E) and a transcritical bifurcation (F). In addition there is canrapidlypropagatecausingalargeoutbreak. Thisleads
a fold bifurcation of cycles (B) and a homoclinic bifurcation backtoastateofhighprevalencecompletingthecycle. In
(D). Two small insets indicate the eigenvalue configuration at
the model with constant rewiring the invasion of the sus-
points A and C. Inset: time series on the limit cycle attractor
ceptible cluster is difficult as it is still protected by rapid
at p = 0.0006. Parameters: w0 = 0.06, r = 0.0002, N = 105,
L=106. rewiring. Iftheamplitudeofthecycleislargethenumber
of infected individuals at the low point is small, which re-
duces the chance of successful invasionfurther. Therefore
consistent with the saddle cycle mentioned above. large amplitude oscillations cannot persist in the model
Atlowerinfectionratesastablesteadystateandasad- with constant rewiring and the parameter range in which
dle point emerge from a fold bifurcation (point C). This oscillations are observedis small. In the model studied in
discontinuoustransitionmarkstheedgeofanarrowinter- this Letter large oscillations also result in small numbers
valinwhichthelimitcycleattractorcoexistswiththenew ofinfected individuals in the low point of the cycle. How-
stationary attractor. At slightly lower infection rates the ever in this case the disease is unlikely to become extinct
limit cycle vanishes in what appears as a homoclinic bi- as the system is still beyond the epidemic threshold and
furcation(pointD).Afterthisbifurcationthesteadystate the protection offered by rewiring decreases proportional
that has emerged from C is the only attractor of the sys- to the number of infected. Since the disease free state is
tem. The twosaddlesvanishinafoldbifurcationatpoint thus dynamically unstable reinfection of the susceptible
E. Finally, at the epidemic threshold (point F), the sta- cluster is bound to occur.
ble equilibrium undergoes a transcritical bifurcation with
the trivial equilibrium (i = lII = lIS = 0) as it leaves the Summary and Conclusions. – In this Letter we
positive cone of the state space. have proposed a computational approach to the investi-
Letus emphasizethat Fig.4 describesthe dynamicson gation of emergent properties of adaptive networks. By
the emergent level: inthenon-trivialstationarystatesthe applying the methods of coarse-grained, ‘equation-free’
networkisnever frozen,butconstantlyundergoeschanges modeling and analysis we have shown that existing nu-
inits topologyandthe statesofits nodes. However,these merical tools of dynamical systems theory can be used to
changesarebalancedinsuchawaythatsystem-levelvari- investigate networks directly on an emergent level. The
ables, such as i, lSS and lII remain effectively stationary. approachenabled us to obtain information that is not di-
Furthermore, note that two of the state variables (lSS rectly accessible by simulation alone, while avoiding the
and lII) describe topological degrees of freedom. With- often prohibitively difficult derivation of closed system-
out adaptive rewiring, it is well known that the emergent level equations of motion. In particular we have avoided
level dynamics of the SIS model can be described by the thestronghomogeneityassumptionthatisinherentinpre-
single variable i and only exhibits transition E [19]. Dy- vious analytical moment closure approximations.
namics suchas the self-organizedcycles for SIS model are Ourapproachdependscriticallyontheassumptionthat
therefore only possible on adaptive networks. bona fide emergent-level variables exist and can be iden-
p-5
T. Gross and I.G. Kevrekidis
tified; else closure (numerical or otherwise) at the desired REFERENCES
level is impossible. The moment-based description used
here is applicable to networks in which the number of [1] Dorogovtsev S. N.andMendesJ. F. F.,Evolution of
different states of the nodes is small and the topological Networks (Oxford University Press, Oxford) 2003.
changecanbeexpressedintermsofthenetworkmoments. [2] Albert R. and Barabasi A., Rev. Mod. Phys. , 74
(2002) 1.
The success of this description depends critically on the
[3] Newman M. E. J., BarabasiA.andWattsD. J.,The
time scale separation between the state variables and the
structureanddynamicsofnetworks(PrincetonUniversity
higher moments of the network. In our model the time
Press, Princeton) 2006.
scaleseparationisfiniteanddecreasesastherewiringrate [4] Baraba`si A. and Albert R., Science , 286 (1999) 509
is increased. At very high rewiring rates the computation
(DOI:10.1126/science.286.5439.509).
of closure terms ceases to converge. [5] Watts D. J. and Strogatz S. J., Nature , 393 (1998)
Inothersystemsitmaybeadvantageoustousedifferent 440.
sets of system-level variables, such as heuristic variables, [6] Pastor-Satorras R. and Vespignani A.,
based on the researcher’s experience or variables identi- Phys. Rev. Lett. , 86 (2001) 3200 (DOI: 10.1103/Phys-
fied by automateddata-mining algorithms [26]. Although RevLett.86.3200).
dynamical equations, such as Eq. (1) are in this case not [7] May R. M. and Lloyd A. L., Phys. Ref. E , 64 (2001)
066112 (DOI:10.1103/PhysRevE.64.066112).
available, one can still estimate temporal derivatives of
[8] Bornholdt S. and Rohlf T., Phys. Rev. Lett. , 84
the chosen set of variables directly from short bursts of
(2000) 6114.
properly initialized simulation.
[9] Zhou C. S.andKurths J.,Phys. Rev. Lett. ,96(2006)
Here we only applied two simple numerical tools: For-
164102 (DOI:10.1103/PhysRevLett.96.164102).
ward Euler integration and the Newton method. These [10] Ito J. and Kaneko K., Phys. Rev. Lett. , 88 (2002)
tools revealed significant information on the coarse bifur- 028701.
cation structure of the system. Demonstrating the ap- [11] Pacheco J. M., Traulsen A. and Nowak M. A.,
plicability of such tools to a network problem constitutes Phys. Rev. Lett. ,97(2006) 258103 (DOI:10.1103/Phys-
a proof of principle. More sophisticated integrators and RevLett.97.258103).
fixed point algorithms, continuation of solution branches, [12] Holme P. and Newman M. E. J., Phys. Rev. E , 74
automated bifurcation detection, computation of normal (2007) 056108 (DOI:10.1103/PhysRevE.74.056108).
form coefficients etc. can in principle be analogously ap- [13] Gross T., Dommar D’Lima C. and Blasius B.,
Phys. Rev. Lett. , 96 (2006) 208701.
plied. These tools are available in the form of free, well-
[14] Gross T. and Blasius B., Adaptive coevolutionary
tested software packages that are commonly used for the
networks – a review j. Roy. Soc. Interface in press
investigationofsystems ofordinarydifferentialequations.
(DOI:10.1098/rsif.2007.1229) (2007).
By writing numerical wrappers that operate along the [15] Rosvall M. and Sneppen K., Euro. Phys. Lett. , 74
lines discussed above, this arsenal of existing, highly effi- (2006) 1109.
cientmethodscanbebroughttobearontheinvestigation [16] Karlen A., Man and microbes: Diseases and plagues in
of adaptive networks. history and modern times (Touchstone, New York) 1995.
In the present Letter we have used the proposed proce- [17] Oldstone M.B.A.,Viruses, Plagues, andHistory(Ox-
dure to investigate an epidemiological SIS model. In our ford University Press, Oxford) 1998.
modelindividualscanprotectthemselvesbyalteringtheir [18] Omi S., SARS – How a global epidemic was stopped
(WHO Press, New York,Geneva) 2006.
topological neighborhood at a rate determined by their
[19] Anderson R. M.andMay R. M.,Infectious diseases of
awareness to the disease. Our analysis shows that in a
humans.DynamicsandControl(OxfordUniversityPress,
considerable parameter range the prevalence of the dis-
Oxford) 2005.
ease and the topology of the network exhibits oscillations
[20] Ehrhardt G. C. M. A., Marsili M. and Vega-
of large amplitude. Let us emphasize that the observed RedondoF.,PhysicalReviewE,74(2006)036106(DOI:
oscillations are not caused by the variable level of aware- 10.1103/PhysRevE.74.036106).
ness, since in our model awareness changes without time [21] KevrekidisI.G.,GearC.W.andHummerG.,AIChE
lag. Instead variable awareness levels extend the parame- Journal ,50 (2004) 1346.
terrangeinwhichoscillationscanbeobserved;decreasing [22] Gear C. W., Hyman J. M., Kevrekidis P. G., Run-
awareness levels at low prevalence stabilize cycles which borg O. and Theodoropoulos K., Comm. Math. Sci.
would otherwise drive the disease to extinction. The sim- , 1 (2003) 715.
[23] Gear C. W. and Kevrekidis I. G.,Physics Letters A ,
ple mechanism studied in our conceptual model may ex-
321 (2005) 335.
plain oscillatory dynamics observed in certain real world
[24] Kuznetsov Y., Elements of applied bifurcation theory
diseases.
(Springer Verlag, Berlin, Germany) 1989.
[25] KelleyC.T.,Iterative methodsforlinearandnonlinear
∗∗∗
equations (SIAM,Philadelphia) 1995.
[26] NadlerB.,Lafon S.,CoifmanR.R.andKevrekidis
This research is partially supported by the Humboldt I. G., Applied and Computational Harmonic Analysis ,
Foundation (TG), DARPA and DTRA (IGK, TG). 21 (2006) 113.
p-6