Fluctuating epidemics on adaptive networks

0.41 MB·English

by  Leah B. Shaw

Fluctuating epidemics on adaptive networks
Leah B. Shaw
Department of Applied Science, College of William and Mary, Williamsburg, VA 23187
Ira B. Schwartz
US Naval Research Laboratory, Code 6792, Nonlinear Systems Dynamics Section,
Plasma Physics Division, Washington, DC 20375
A model forepidemics on an adaptivenetworkis considered. Nodesfollow an SIRS(susceptible-
infective-recovered-susceptible) pattern. Connections are rewired to break links from non-infected
nodestoinfectednodesandarereformed toconnecttoothernon-infectednodes,asthenodesthat
arenotinfectedtrytoavoidtheinfection. MonteCarlosimulationandnumericalsolutionofamean
8 field model are employed. The introduction of rewiring affects both the network structure and the
0 epidemic dynamics. Degree distributions are altered, and the average distance from a node to the
0 nearest infective increases. The rewiring leads to regions of bistability where either an endemic or
2 adisease-free steady state can exist. Fluctuations around theendemicstate andthelifetime of the
endemicstate are considered. The fluctuationsare found to exhibit power law behavior.
n
a
J
I. INTRODUCTION demic spread and in turn affect future spreading of the
3
epidemic. Here, the new parameter is one that describes

the rewiring rate of the network, which is controlled by
] The study of recurrent epidemics has a long history
E thefractionofsusceptible(S)-infective(I)links. Thenet-
[1],andmanymodels,bothdeterministic andstochastic,
P workaltersdynamicallywhentherearecontactsbetween
have been considered. Deterministic models have been
. S and I, and social pressures (the desire to avoid illness)
o usedsincethetimeofBernoulliandhaveexplainedsome
bi of the mechanisms in the spread of the infectious dis- rewirethe contactstobe insteadbetweenSandS.Infec-
tionsarereducedduetoisolation,andanewphenomenon
- eases. However, deterministic models are not sufficient
q occurs: for appropriatechoicesofparameters,bistability
to account for some of the important stochastic dynam-
[ between the disease free equilibrium and endemic state
ics, such as extinction [2, 3] and sustained fluctuations
has been observed. This is in contrastto static networks
1 [4]. From the general theory of finite Markov chains
ina largepopulation,where there is typically only a sin-
v [5], it was shown that in stochastic models the proba-
6 bilityofextinctionisequaltooneintheasymptotictime gle attracting endemic or disease free state.
0 In this paper, we introduce a recovered, immune class
limit. Numerical [6, 7, 8] and analytic [9] comparisons
6 and consider this slight generalization of the SIS model
of stochastic and deterministic models have been per-
0 on an adaptive network. We examine the structure of
. formed. Thenumericalresultsholdforverysmallampli-
1 the network and the dynamics of the fluctuations of the
tude noise as well as real finite noise. Deterministic SIS
0 epidemic. Our approachis to combine Monte Carlosim-
orSIRSmodelsresultinanequilibriumendemicpresence
8 ulations and stochastic mean field models for epidemic
ofinfectivesforanappropriatechoiceofparameters. Itis
0
evolution on evolving networks. The layout of the paper
: clear that stochastic effects may result in very different
v is as follows. We introduce the model in Section II and
dynamics from deterministic models, particularly when
i present its bifurcation structure in Section III. Proper-
X fluctuations and/or extinction occur.
ties ofthe networkstructure arediscussedin SectionIV.
r Morerecently,thestudyoffluctuationsandthespread
a Wediscussdynamicalpropertiesofthesystem,including
of simple models of epidemics have been simulated on
fluctuations and lifetimes of the states, in Section V.
large networks [10, 11, 12, 13, 14]. In almost all of these
networkmodels, the epidemic propagatesona fixed net-
work. The epidemic dynamics is typically studied as an
SIS or SIR model, in which the population is large and II. MODEL
isolated. In addition to the dynamics on such fixed ar-
chitectures,controlsbasedonvaccinationhavebeencon- We study an SIRS (susceptible-infective-recovered-
sidered as well [15, 16]. Several recent models have con- susceptible) model on an adaptive network. Epidemic
sidered epidemics on a network that changes structure dynamics on the nodes is as follows. The rate for a sus-
dynamically according to rules that do not depend on ceptiblenodetobecomeinfectedispN ,whereN
I,nbr I,nbr
the nodes’ epidemic status [17, 18]. is the number of infected neighbors the node has. The
Incontrasttothemodelsofastaticnetworkormodels recoveryratefor aninfectednode is r. A recoverednode
with externally applied changesin structure,a new class becomes susceptible again with rate q, which we define
ofmodelsbasedonendemicSISpopulationsonanadap- as the re-susceptibility rate.
tive network has been recently introduced [19]. Changes While the epidemic spreads, the network is also be-
tothenetworkstructurearemadeinresponsetotheepi- ing rewired adaptively. If a link connects a non-infected
2
node to an infected node, that link is rewired with rate −pKPSIPSR +w PS P (7)
N PS PS+PR IR
w to connectthe non-infectednode to anotherrandomly
P˙ = 2rP +pKPSIPSR −qP −rP
selected non-infected node. Self links and multiple links IR II N PS IR IR
between nodes are disallowed. −wPIR (8)
In examining steady state solutions, it is sufficient to P
P˙ = rP −2qP +w R P (9)
fix one of the rates,as time may be rescaledaccordingly. RR IR RR IR
P +P
S R
For this reason, we fix r =0.002 throughout this paper.
We performed Monte Carlo simulations of this model The ordinary differential equations of the mean field
on a system with N = 104 nodes and K = 105 links. model can be integrated easily with any well known nu-
(Larger system sizes with the same node to link ratio merical integration technique. We choose initial condi-
were also considered. The major results of this paper tionsso thatwearenearanendemic state. We alsonote
do not depend strongly on system size.) In each Monte that the model does support solutionswith negative val-
Carlo step (MCS), we randomly select N nodes and M ues, but these areunphysical, andso we ignorethem. In
links, where M is the number of links that may poten- thecaseofstochasticsimulations,wehaveconsideredthe
tially rewire (susceptible-infected and recovered-infected effects of both internal fluctuations, modeled as multi-
links), and the links are selected from the pool of links plicative noise, as well as external fluctuations, modeled
thatmayrewire. Initialconditionsareconstructedinone as additive noise. We use a fourth order Runge-Kutta
oftwoways. Weeithergeneratearandom(Erd¨os-R´enyi) solverforeachofthese casesto generatestochasticstim-
graphofsusceptibles andconvertafractionf ofthem to ulations of the mean field. We also tracked the steady
infectives, or we use the final state of a previous run as states as a function of parameters using a continuation
an initial condition. Transients are discarded and simu- package [20].
lationsrunlongenoughthattheinitialconditionsdonot
affect the results.
Following [19], we also developed a corresponding III. BIFURCATION STRUCTURE
mean field model for the system. The mean field model
tracksthedynamicsofbothnodesandlinks. P denotes
A We first consider the steady state bifurcation struc-
the probability of a node to be in state A, where A is ei-
ture of the model. In Figure 1, we show examples of
ther S (susceptible), I (infected), or R (recovered). P
AB theaverageinfectedfractionversusthetransmissionrate
denotestheprobabilitythatarandomlyselectedlinkcon-
p. Two steady states can occur, a disease-free state and
nects a node in state A to a node in state B. We obtain
an endemic state. In the absence of rewiring (Fig. 1a),
the following mean field equations for the evolution of
the disease-free state loses stability for very small trans-
the nodes:
mission rates, and only the endemic state is observed at
P˙ = qP −pKP (1) larger p values. When rewiring is introduced (Fig. 1b),
S R N SI
the disease-free state is stabilized for larger p values. A
P˙ = pKP −rP (2)
I N SI I region of bistability, in which both endemic and disease-
P˙ = rP −qP (3) free states are observed, now occurs. Bistability was ob-
R I R
served previously for the SIS model [19]. Monte Carlo
For example, in the first equation, recovereds are con- simulation points in Fig. 1 were computed as follows.
verted to susceptibles with rate q, and infection spreads To locate the upper branches (endemic state), we swept
withratepalongeachsusceptible-infectedlink. Rewiring p from larger to smaller values, using the final state of
does not appear directly in the node equations, since each run as the initial state for the next run. 5×103
rewiring operates on links, but it affects the system im- MCS of transients were discarded, and the steady state
plicitly through the number of susceptible-infected links was averaged over 104 MCS. To locate the lower branch
(KPSI). We next write a system of mean field equa- (diseasefreestate),webeganwitharandomlyconnected
tions for the links. To close the system, we follow [19] networkinwhichafractionf ofthe nodeswereinfected,
and make the assumption for three point terms that while the rest were susceptible. The system was simu-
PABC ≈ PABPBC/PB. This assumption leads to the lated for 2×104 MCS. f values between 0.025 and 0.9
following system of equations for links: weretried,andatleast5runsweredoneforeachpvalue.
If the epidemic died out in any of the runs, the disease
P
P˙ = qP +w S P −2pKPSSPSI (4) free state was considered stable. Due to the stochastic
SS SR P +P SI N PS
S R nature of the Monte Carlo simulation, and because the
P˙SI = 2pKN PSPSSPSI +qPIR−rPSI −wPSI disease-free state is absorbing, stability designations are
−p P + KPS2I (5) uncertain. It is difficult to distinguish a weakly unstable
(cid:16) SI N PS (cid:17) state from a weakly stable state with a short lifetime.
P˙ = p P + KPS2I −2rP (6) Lifetimes of the endemic state are considered in more
II (cid:16) SI N PS (cid:17) II detail in Section V.
P TheresultsinFigure1showfairlygoodagreementbe-
P˙ = rP +w R P +2qP −qP
SR SI P +P SI RR SR tweenthemeanfieldapproximationandtheMonteCarlo
S R
3
(a) 0.4
0.5
LP(cid:9) (a)
0.4 BP
0.3 HB II
0.3
I
I 0.2
w0.2
0.1
III
0
0.1
−0.1
0 0.002 0.004 0.006 0.008 IV
p
(b) 0
0 2 4 6 8
p −3
x 10
0.4
0.3
1
I 0.2
(b)
0.1
0.8
0
−0.1
0 0.002 0.004 0.006 0.008 0.6
p
I
0.4
FIG. 1: Average infected fraction vs. transmission rate p.
(a) Static network, w = 0. (b) Rewired network, w = 0.04.
Black dots: Monte Carlo simulations; solid gray lines: mean 0.2
fieldsolution (stablebranches);dashedgraylines: meanfield
solution (unstablebranches). q=0.0016, r=0.002.
0
0 2 4 6 8
p −3
x 10
simulation of the full system, and we typically see this
level of agreement in the steady state values, although
FIG.2: (a)Twoparameterplotforwandpofthebifurcation
as we discuss later, the stability and type of bifurcation points for steady states when q =0.0064. The heavy dashed
sometimes differs. Using the mean field model, we next lineisthelineofsaddle-nodepoints(limit points). Thesolid
explorethebifurcationstructureofthesystemforawider line denotes the Hopf bifurcation points, and light dashed
range of parameters. linedenotesthetranscriticalbifurcationpoints(BP).(b)Aa
bifurcation diagram of the infective fraction as a function of
An interesting property of the steady state instabili-
p,wherew=0.04. Thesquaresdenotethesaddle-nodepoint
ties appears when one considers each of the steady state
and transcritical point. Dashed lines are unstable branches.
bifurcation points of the mean field equations. (There
do exist branches of periodic orbits, but since they oc-
cur within a very small range of parameters, we ignore
them in this paper. They will be treated elsewhere.) If whilethe solidlinesdepicttheattractors. Forlowvalues
the re-susceptibilityrateq is heldfixedandabifurcation ofpthediseasefreesteadystateisstable. Trackingalong
diagram constructed, we find the existence of at least the diseasefree branch,ata criticalvalue ofp, anunsta-
two distinct regimes for different q values, illustrated in blebranch(subcritical)ofendemicsteadystatesappears.
Figs. 2 and 3. The instabilities appear as a transcriti- The endemic branch becomes stable at the saddle-node
cal bifurcation from the disease free state steady state, point. There exists a clear region of bistability with co-
a saddle-node bifurcation of endemic steady states, and existing endemic and disease free states for a range of p.
a Hopf bifurcation, from which a branch of subcritical If we now vary the parameter w and track each bifurca-
unstable periodic orbits emanates (not shown). tion curve, we obtain the result in Fig. 2a. We describe
In the case where q = 0.0064, depicted in Fig. 2b for the bifurcationregionsindetail for w largerthan 0.1. In
w =0.04,weshowatypicalbifurcationplotofstableand region I, we have only a stable disease free equilibrium.
unstablesteadystatesastheinfectionratepisincreased. As we cross into region II for large w, the disease free
The dashed lines represent the unstable steady states, equilibrium is stable, and there exists an unstable en-
4
demic state. Region III exhibits bistability between the
disease free equilibrium and endemic state, and region 0.2 (a)
IV has just a stable endemic equilibrium. We note that LP
at w = 0.04, we have the simple saddle-node transition 0.15 HB I II
depicted in Fig. 2b, since there is no Hopf bifurcation to BP
periodic cycles at that particular w value.
The previous discussion presented a case for re- w 0.1 III
susceptibility rate q where the limit (saddle-node) and
IV
Hopf bifurcation branches are close to each other. The
distinction between the saddle-node and Hopf branches 0.05
can be seen more easily if q is lowered to 0.0016, as
shown in Fig. 3a. In Fig. 3a, for sufficiently large w,
as p is increased in region I the system first undergoes 00 2 4 6 8
a limit point bifurcation, and then a Hopf bifurcation as p x 10−3
it passes through region II. The Hopf curve is actually a
closed isola in two parameters. The limit point here is
0.5
a saddle-saddlepoint, where a steady state having a two (b)
dimensionalunstablemanifoldconnectstoasteadystate
withaonedimensionalunstablemanifold. Inbothcases, 0.4
wehavebistablebehaviorforw sufficientlylarge,butthe
region of bistability is much smaller since the Hopf and
0.3
transcriticalbranchesareclosertogetherforthisvalueof
q (region III). I
Forw =0.04,themeanfieldendemicsteadystateloses 0.2
stability in a saddle-node bifurcation for q =0.0064 and
in a Hopf bifurcation for q = 0.0016. In our discussion
offluctuations in the endemic state in SectionV, we will 0.1
refer to q = 0.0064 because the saddle-node bifurcation
structure best corresponds to the scaling of fluctuations
0
that we observe in Monte Carlo simulations of the full 0 2 4 6 8
system. p x 10−3
FIG. 3: (a) Two parameter plot of the steady state bifurca-
tionpointswhenq=0.0016. Theheavydashedlineistheline
IV. NETWORK GEOMETRY
of limit points. The solid line denotes the Hopf bifurcation
points, and the light dashed line denotes the transcritical bi-
A. Degree distributions
furcations. (b)Abifurcationdiagramoftheinfectivefraction
asafunctionofp,wherew=0.04(samebifurcationdiagram
Rewiringleadstosignificantalterationsinthenetwork in Fig. 1).
structure. We firstconsiderthe degreedistribution. Fig-
ure 4 shows degree distributions for each type of node
increase,asothersusceptibleandrecoverednodeswireto
in the absence (4a) and presence (4b) of rewiring. In
connect to it. This leads to the following equations for
Fig. 4b, we averaged over 3×104 MCS. In Fig. 4a, the
d :
network is static, so we averaged over 10 separate runs R,n
to obtain better statistics. d
(d ) = rd −qd −kd
Using mean field ideas, we can understand the recov- dt R,0 I,0 R,0 R,0
ered and susceptible degree distributions as nodes flow d
from infected to recovered to susceptible. We outline (dR,n) = rdI,n−qdR,n−kdR,n (10)
dt
thecalculationbrieflyhere. However,thisapproachdoes
+kdR,n−1 for n>0
notaccuratelypredictthedegreedistributionforinfected
nodesbecausecorrelationsplayamoreimportantrolefor wherek isthe averageratefornodestorewiretoagiven
thesenodes,aswewillexplainbelow. Forthisreason,we non-infected node. k is given by the ratio of the total
cannotwrite a self-consistentsetofequations for the de- rewiring rate to the number of potential target nodes:
greedistributionsthatcouldbesolvedwithoutinputting
wK(P +P )
SI IR
simulation results. k= . (11)
N(P +P )
Let d be the number of nodes in state X (either S R
X,n
S, I, or R) with degree n. Recovered nodes originate Given the degree distribution of the infectives d and
I,n
when infectives recover and are lost when they become the probabilities appearing in k, Eqs. 10 can be solved
susectibleagain. Thedegreeofarecoverednodecanonly for the degree distribution of the recovereds.
5
wherethefractionistheratioofthenumberoflinksthat
(a) cantransmitinfectiontothe totalnumberoflinksinto a
0.1 susceptible. As with the recovereds,the steady state de-
gree distribution for susceptibles can be computed from
Eq. 12. The predicted degree distributions for suscepti-
0.08
bles and recovereds are overlaid on the actual distribu-
y tions in Fig. 4b, using Monte Carlo simulation averages
bilit 0.06 for the infective degree distribution and the node and
a
b link probabilities. (Note that the node and link prob-
o 0.04
pr abilities could instead be obtained from the mean field
system.) Deviations between the prediction and simula-
0.02 tion are smaller than the width of the curves in Fig. 4b,
so they are indistinguishable.
0 The degree of infected nodes, however,cannot be pre-
0 20 40 60 80
dicted by this approximate procedure. We might expect
that
(b) d ′′
(d ) = k d −rd
I,0 I,1 I,0
dt
0.2 d ′′ ′′
(d ) = k (n+1)d −k nd (14)
I,n I,n+1 I,n
dt
′
0.15 +k ndS,n−rdI,n for n>0
y
bilit where
a
b 0.1
pro k′′ =w PSI +PIR (15)
P +P +2P
SI IR II
0.05
istheperlinkrewiringrateforlinksconnectingtoanin-
fective (i.e., the ratio between links that can potentially
0 rewire and total links into infectives). Figure 5a com-
0 20 40 60 80
degree pares the actual degree distribution for infectives with
the distributionpredicted by Eq.14. The number of low
degree infectives is significantly over-predicted. This oc-
curs because the mean field approximation in Eq. 15 is
FIG.4: DegreedistributionsfromMonteCarlosimulationfor not accurate for infectives. Figure 5b shows the frac-
p = 0.002, q = 0.0016, r = 0.002. (a) Static network, w = tion of infected neighbors that a node has, depending
0. All node types have Poissonian degree distribution. (b)
on its degree and disease status. Results are averaged
Rewired network, w = 0.04. Solid gray: infectives; dashed:
over 3 × 104 MCS. Low degree infected nodes tend to
recovereds; black: susceptibles.
haveamuchhigherfractionofneighborsthatarealsoin-
fected,duetotransmissionofthedisease. Onceinfected,
Susceptible nodes originate when recovereds become these neighbors will not rewire away until recovered, so
susceptible again, and they may be lost when they be- the rewiring rate per link is smaller than one would ex-
′′
comeinfectedbyaneighbor. Aswiththerecovereds,the pect from the mean field k in Eq. 15. Predicting the
degree of a susceptible increases due to rewiring. Thus infective degree distribution accurately would require a
the time evolution of the degree distribution for the sus- theory that accounts for these correlations in infection
ceptibles can be written as status of neighboring nodes, which is beyond the scope
of the present work.
d
(d ) = qd −kd
S,0 R,0 S,0
dt
d ′ B. Distance from an infective
(d ) = qd −k nd −kd (12)
S,n R,n S,n S,n
dt
+kdS,n−1 for n>0 We next consider the distribution of distances from a
′ given node to the nearest infective. These distances are
where k is the infection rate per link into a suscepti-
′ ofinterestbecause they relateto the numberofhops the
ble node. We assume that k is independent of degree,
disease must make in order to reach an uninfected indi-
whichweknowfromsimulationsisapproximatelycorrect
vidual. The disease cannot propagate through recovered
(cf. Fig. 5), and write
nodes until they become susceptible again, so the dis-
P tancefromthe nearestinfective does notnecessarilycor-
′ SI
k =p (13)
P +P +2P respond to a path for disease propagation. However, we
SI SR SS
6
(a)
(a)
8000
0.6
6000
0.5 er
billity 0.4 numb 4000
a
b 0.3
o 2000
pr
0.2
0
0.1 0 2 4 6 8 ... ∞
distance
0
0 20 40 60 80
(b)
degree 104
(b)
102
1
er
n 0.8 mb 100
o u
cti n
d fra 0.6 10−2
e
ct
br. infe 0.4 10−40 2 4 distanc6e 8 ... ∞
n 0.2
0
0 20 40 60 80 FIG. 6: Distribution of distances from the nearest infected
degree node. ∞ indicates nodes that are completely disconnected
from an infective. (a) Solid black: with rewiring (w = 0.04,
p = 0.002, q = 0.0016); dashed gray: no rewiring (w = 0,
p = 0.002, q = 0.009). (b) Black: rewired case (same as in
FIG.5: (a)Actual(black)andpredicted(dashedgray)degree (a)); gray: distribution for random graphs, as described in
distributions for infecteds. (b) Average fraction of neighbors text.
that are infected vs. degree. Solid gray: infecteds; dashed:
recovereds; black: susceptibles. p = 0.002, q = 0.0016, r =
0.002, w=0.04.
Given how far from Poisson the degree distributions are
when the network is rewired, it is somewhat surprising
thattheformofthedecayinthedistributionofdistances
note that rewiring acts only on links to infectives, and canbepredictedfromrandomnetworks. Beginningwith
thus the chains of susceptible and recovered links that S, I, and R nodes in numbers matching that observed
this metric identifies will persist until the infection leads in the averageof Monte Carlo simulations, we generated
totheirinterruption. Figure6ashowsthe distributionof 1000randomnetworksandaddedrandomlyselectedlinks
distances from the nearest infective in the presence and untilthenumberofS-S,S-I,S-R,etc.,linksalsomatched
absenceofrewiring. Todisplaythe effectonthe network that from the average of Monte Carlo simulations. As
geometry alone, rather than on the steady state number Figure 6b shows, the distribution of distances for these
of infections as well, we have used a smaller q value for randomnetworksdecaysinthesamewayasitdoesinan
the w = 0 case so that the total number of infectives is adaptive network. Thus the form of the distribution of
approximatelythesameinbothcurves. Resultswereav- distances depends mainly on local dynamics (node and
eraged over 3×104 MCS, sampled every 100 MCS after linkdynamics)ratherthanonthe detailsoflongerrange
removingtransients. Rewiringsignificantlydecreasesthe correlations. The main difference is that the adaptive
number of nodes that are directly connected to an infec- networkhassomenodesfully disconnectedfrominfected
tive. However,despitetherewiring,onlyasmallfraction components,while therandomnetworksdonot. Mostof
of nodes are fully disconnected from the infection. these nodes are recovereds with degree 0, which appear
Figure 6b shows a semilog plot of the w = 0.04 case wheninfectivesofdegree0recover. Eachtypeofnodein
to display the tail of the curve at larger distances. An the random networks has a Poisson degree distribution,
approximationbased on random networks is also shown. so they do not generally have nodes of degree 0.
7
V. FLUCTUATIONS AND OTHER DYNAMICS meanfieldhasdifferentq value,thebifurcationstructure
isequivalenttothatobservedinMonteCarlosimulation,
so we use q =0.0064 in the fluctuation study.
In the previous sections, we have considered steady
The computations reveal that the fluctuations exhibit
states and long time averagesof network properties. We
powerlawscaling,asshowninthelog-logplotsofFig.7b
next consider fluctuations and dynamical properties of
and7d. Onthehorizontalaxis,weplotln(p−p ), where
the endemic state. c
p is the critical point at which the endemic state loses
c
stability. For the mean field, the bifurcation point is
knownexactlybyexaminingtheeigenvaluesoflinearized
A. Fluctuations near bifurcation point
vector field at the steady state. We approximate the
Monte Carlo bifurcation point p as the value that pro-
c
Near the bifurcation point where the endemic state
duces the most linear plot. Although both cases have
loses stability, the number of infectives has larger fluc-
powerlaw scaling,the exponents aredifferent: −0.59for
tuations due to noise overcoming weak attracting forces
MonteCarloand−0.27formeanfield. Thescalingexpo-
[21]. Fluctuations in the SIRS model are significantly
nentforthefullsystemdependsonthenumberofnodes.
larger than those in the previously studied SIS model.
Whetheritwillapproachthemeanfieldvalueinthelimit
We quantify the fluctuations by computing the standard
of infinite system size is a subject for further study. We
deviationdividedbythemeanforlongtimeseriesinboth
note that fluctuations near a Hopf bifurcation point, as
the Monte Carlo and mean field simulations. In Fig. 7,
occurs in the mean field for q =0.0016,wouldproduce a
we plot the fluctuations as the infection rate p is swept
very different form of scaling from the saddle-node case.
towards the bifurcation point. Monte Carlo results were
One of the reasons for the difference is that the instabil-
computed from 5×105 MCS time series sampled every
ity is then two dimensional and underdamped [23]. This
10 MCS, except for the two smallest p values, for which
is known to cause very different scaling laws in generic
shortertime series wereused due to the shorterlifetimes
problems, which can be much slower [24].
ofthesestates. Alltimeserieswerelongerthan105MCS.
To motivate the power law scaling of the fluctuations,
For comparison, the mean field equations can also be
we consider scaling near a generic saddle-node bifurca-
considered near equilibrium in stochastic form. In gen-
tion. The simplest generic case of a saddle-node bifur-
eral, near equilibrium fluctuations can be modeled as
cation for equilibrium points comes from solving for ze-
additive noise [22], and we do so here. (We do note
roes of the vector field at a parameter value where one
that multiplicative noise effects generate similar results
eigenvalueoftheJacobianpassesthroughzero. Standard
to those reported for additive noise.) We assume the
normal form analysis allows one to consider the generic
mean field is of the following form:
problem of a saddle-node bifurcation. In one dimension,
the stochastic differential equation of a saddle-node bi-
X′ =F(X)+ǫη(t), (16)
furcation may be modeled as
where F(X) is the mean field system in Eqs. 1-9, and dx =(a−x2)dt+σ∗dW . (17)
hη(t)η(t′)i=δ(t−t′). ǫ is the noise strength, or ampli- t t t
The parameter is a, and we suppose noise is additive.
tude. We haveconsideredbothadditivenoiseandmulti-
Since noise in general may cause a “shift” in parame-
plicative noise cases in the simulations of the stochastic
ter values where the saddle-node point disappears, we
attractors near the endemic state and have computed
assume that the noise near the bifucation is sufficiently
the standarddeviationdividedbythe meanasdescribed
small, where dW/dt is a white noise term, and dW is a
above for 10 random initial conditions near equilibrium
Brownian increment.
and 10 realizations.
Wefurtherassumethatwearealwaysneartheattract-
We recall from Figs. 2 and 3 that depending on the
ing branch of the saddle-node, so we are in a near equi-
valueofthere-susceptibilityrateq,thebistabilityregions
librium setting. Such an assumption allows us to exam-
III are vastly different. Specifically, for q = 0.0016, we
ine the stationary probability density function (PDF) of
saw that for sufficiently large wiring rates, the saddle-
thestochasticdynamicsbyemployingtheFokker-Planck
saddle and Hopf bifurcation branches were well sepa-
equationnearsteadystate. Forthestochasticdifferential
rated, whereas for q = 0.0064, the branches were very
equation,Eq.17,thePDFiswellknown[21]andisgiven
closeforsmallvaluesofre-wiringratew. InMonteCarlo
by
simulation,wehavenotobservedtheHopfbifurcationsor
stableperiodicoscillationsseeninthemeanfield,evenfor p(a,x,σ)=Ne2(ax−x3/3)/σ2. (18)
systemsizesaslargeas4×105nodes. Althoughthevalue
of q =0.0064 we use in the mean field fluctuation study Here N is a normalization constant. We compute the
is different from that used in the Monte Carlo, the local firstandsecondordermomentsdirectlyusingEq.18and
bifurcation structure is similar when w = 0.04, in that thentaketheratioofthestandarddeviationtothemean.
it is a true saddle-node bifurcation point. This has been Since a=0 is the value of the saddle-node point, we ex-
checked by examining the local linear vector field at the aminethefluctuationsintheneighborhoodofthatvalue.
saddle-node point in question. Therefore, although the The results, shown in Fig. 8, display power law scaling.
8
(a) (b)
0.4 −1
−1.5
0.3
)
n −2
n a
a e
e m
m 0.2 d/ −2.5
d/ st
st g( −3
o
0.1 l
−3.5
0 −4
0.002 0.003 0.004 −10 −9 −8 −7 −6
p log(p−p )
c
(c) (d)
0.01 −4
−4.5
0.008
)
n −5
n a
a 0.006 e
e m
m d/ −5.5
d/ 0.004 st
st g( −6
o
l
0.002
−6.5
0 −7
0 0.005 0.01 −12 −10 −8 −6 −4
p log(p−p )
c
FIG. 7: Fluctuations in infectives (standard deviation divided by mean) vs. infection rate p near the bifurcation point: (a)
MonteCarlo,(b)meanfield. Curvesaretoguidetheeye. Log-logplots(datapointswithbestfitlines)showpowerlawscaling
forbothMonteCarlo(b)andmeanfield(d). MonteCarloparameters: q=0.0016,w=0.04,r=0.002. Meanfieldparameters:
q=0.0064,w=0.04,r=0.002,ǫ=0.0001.
B. Delayed outbreaks Cross correlations between the infectives and the infec-
tive neighbors were computed for varying shifts between
the time series, and the lag maximizing the cross corre-
Wenextconsiderphaserelationshipsbetweenthefluc-
lationwasidentified. AsshowninFig.9b,rewiringleads
tuating variables. We tracked the number of infectives
to increasing lag times and delayed outbreaks. We also
in the system at each time point as well as the number
computedtime lagsforthe meanfieldmodelwith multi-
of nodes that neighbored an infective (i.e., the number
plicative noise and found the same trend of increasingly
of SI and IR links). In the rewired system, fluctuations
delayed outbreaks with larger rewiring.
in the number of infectives laggedbehind fluctuations in
thenumberofinfectiveneighborsthatarenotthemselves
infected, as shown in Fig. 9a.
C. Lifetime of endemic steady state
Both mean field and Monte Carlo simulations of the
full system displayed this effect, and we studied its de-
pendence on the rewiring rate. Monte Carlo simulations The final dynamic effect we consider is the lifetime of
weresampledevery1MCSfor3×104MCSafterdiscard- theendemicsteadystate. Becausethesystemisstochas-
ing transients. The meanfield wassampled every1 time tic and the disease free state is absorbing, all parameter
unit for 5×104 units. Additive noise was included in valueswillleadtoeventualdieoutofthediseaseinthein-
the mean field equations with noise strength ǫ=0.0001. finite time limit. These lifetimes become shorter and die
9
−0.5
PDF model 13
linear fit
−1
12
n)
a −1.5
e
m T
d/ g 11
g(st −2 lo
o
l
10
−2.5
9
−−37 −6.5 −6 −5.5 −5 −4.5 0 2 4 6
log(a) (p−p )3/2 [× 10−8]
0
FIG. 10: Dependence of endemic state average lifetimes on
FIG. 8: Fluctuation size of a generic saddle-node bifurcation infection rate p. Points: Monte Carlo simulations; line: best
as a function of bifuration parameter a near the bifurcation fit line. q=0.0016, r=0.002, w=0.04. See text for details.
point usingEq. 18. The noise level used is 0.005.
out is more easily observed in the bistable regime, near
the bifurcation point where the endemic state has weak
(a)
stability. We measured the dependence of the lifetimes
on the infection rate p. For each p value, we prepared
s. a steady state initial condition and computed multiple
nbr duplicate runs to obtain a distribution of lifetimes. (We
nd computed 100 duplicate runs for all but the two high-
a
s est p values and over 25 runs for the two highest.) We
e
v then calculated the average lifetime T for each p. For
ecti a generic saddle-node bifurcation in one dimension, it
inf is expected that lnT is a linear function of (p−p0)3/2,
where p is the bifurcation point [25, 26]. We obtained
0
0 200 400 600 800 1000 spurious results for our system with 104 nodes, possibly
MCS becausethesmallsystemsizeledtorapiddieoutandwe
wereunabletoobtaingoodstatisticsnearthebifurcation
(b)
point. Awayfromthebifurcationpoint,thesystemhasa
100
weakly damped oscillatory component and behaves like
80 a focus. By switching to a system with 4×104 nodes
and4×105 links,we wereable to runlongersimulations
S) 60 closer to the bifurcation point and operate in a regime
C
M where only one dimension mattered and the oscillations
g ( 40 could be ignored. Preliminary scaling results are shown
a
l in Figure10. We usedthe bifurcationpoint p that gave
0
20
the best fit line, which led to an R value of 0.99. The
scaling results appear consistent with expectations, but
0
0 0.02 0.04 0.06 0.08 0.1 more study and better statistics are needed.
w
VI. CONCLUSIONS AND DISCUSSION
FIG. 9: Delayed outbreaks due to rewiring. (a) Monte Carlo
time series. Black: infectives; gray: neighbors of infectives. We have explored the stable states, network proper-
Curves are scaled in arbitrary units for comparison of peak ties, and dynamics of an SIRS model on a network with
times. p=0.0065,w=0.09,q =0.0016,r =0.002. (b) Time
adaptive rewiring. As with the SIS model studied previ-
in MCS by which infectives lag behind infective neighbors
ouslybyGrossetal.[19],therewiringleadstobistability
vs. rewiring rate. Solid black: Monte Carlo; dashed gray:
of the endemic and disease free states. A mean field ver-
mean field. p = 0.0065,q = 0.0016,r = 0.002, mean field
sion of the model predicted the steady states with good
noise strength ǫ=0.0001.
numericalaccuracyandwasalsovaluableinstudyingthe
10
fluctuationsofthesystem,withthecaveatthatonemust Monte Carlo simulation agrees qualitatively with the es-
be near the appropriate type of bifurcation in the mean capetimes andyields awell-knownpowerlaw. However,
field to obtain corresponding results. With the addition this is for parameters in which endemic and extinction
of the recovered class and re-susceptibility rate, we can states are not too far apart. Such extinction regimes
controlthe width of the bistability regionby manipulat- can be analyzed using a Fokker-Planck approach [27].
ing the location of the bifurcation points. Although preliminary results are in agreement with the
The fluctuations in the infectives near the bifurca- expected scaling, further study of the lifetime scaling is
tion point showed power law scaling. This agreed with needed,includinginregimesthatarephysicallyrealistic,
mean field results and our expectations for scaling near and where the usual extinction rates cannot be modeled
a saddle-node bifurcation. with a Fokker-Planck approach.
We studied the effects of the rewiring on the network In addition to the directions for further researchmen-
geometry. Degree distributions were altered, and mean tioned above, a major challenge is to develop network
fieldargumentswereable to predictthe distributions for geometries and rewiring rules that are more consistent
susceptibles and recovereds. However, a new analytical with human social networks. Real social networks are
approachthatincludescorrelationsisneededtofullyun- expected to have community structure [28]. The net-
derstand the degree distributions. works studied here did not. In our rewired networks,
The other network property we considered was the the number of connections from non-infected nodes to
distribution of distances from non-infected nodes to the infected nodes was reduced in comparison to a random
nearest infective, a quantity that may be important in network, while the number of connections between non-
disease spreading and its control. This distribution de- infectednodeswasincreased. However,thisprocessdoes
pendedprimarilyonnodeandlinkdynamicsratherthan not induce a community structure on the network in the
onhigherordercorrelations,soitcouldbepredictedfrom Newman-Girvan Q-modularity sense [29]. Rewiring was
random graphs. It is possible to compute the distribu- a non-local effect; new neighbors were chosen at random
tion of distances from an infective analytically for ran- amongst all non-infected nodes in the network rather
dom graphs, but this calculation is awkward for a three thanintroducingalocalstructure. Ithasbeenshownfor
speciessystem(S,I,R)anddoesnothaveasimplefunc- static networks that community structure affects disease
tional form, so we have omitted discussion of analytical dynamics [30, 31], and we expect an impact in adaptive
results here. networks as well.
Delayed outbreaks were observed in the rewired sys- In the current work, an ideal setting was proposed
tem. Peaks in the infective fraction laggedbehind peaks where non-infected nodes were assumed to behave ra-
in the number of nodes that neighbor an infective. For tionallyandhaveperfectknowledgeofthe diseasestatus
the parametervalues studied here,the lag time is onthe of their current neighbors and potential new neighbors.
order of 10% of the mean infectious period, which might It would be of interest to consider a situation in which
be considered a very short lag time. However, the pa- not all contagious individuals appear ill or know that
rameters in this study were not selected to correspond they are contagious, as might be the case for a sexu-
to any specific disease. For most real diseases, we would ally transmitted disease possessing asymptomatic indi-
expect a much slower rate for immunity to wear off and viduals.. This effect might be modeled by simultaneous
recovered individuals to become susceptible again, com- spreading through the network of both the disease and
pared to the mean infectious period. This regime would informationaboutthedisease. Ifthecurrentmodelisex-
be more difficult to study in Monte Carlo simulations, tendedwithinformationandcommunitystructure,social
since the average number of infected nodes would be dynamicscouldbe extrapolatedtoimprovecontacttrac-
much smaller than seen here. Further work is needed to ing and epidemic control in organized populations with
determine whether the observation of delayed outbreaks local structure.
due to rewiring would persist or perhaps become more
significant in a physically realistic system.
Finally, we considered lifetimes of the endemic state
Acknowledgments
near the saddle-node bifurcation where it loses stability.
Inordertoachieveextinctionfromasteadystate,thedis-
ease must first overcome the attractive forces, which are This work was supported by the Office of Naval Re-
weaknearthebifurcationpoint. Duetothegenericlocal search, Center for Army Analysis, and Armed Forces
topology of the saddle-node structure, the escape rate is Medical Intelligence Center. The authors wish to thank
well characterized analytically [24]. We found that the Thilo Gross and Bernd Blasius for helpful discussions.
[1] R. M. Anderson and R. M. May, Infectious Diseases of Bernardino, J. M. Pacheco, and M. C. Gomes, Journal
Humans. (Oxford UniversityPress, 1991). of Theoretical Biology 233, 553 (2005).
[2] J. Verdasca, M. M. Telo Da Gama, A. Nunes, N. R. [3] M. J. Keeling, Ecology, Genetics, and Evolution. (Else-