Critical behaviors in contagion dynamics

1.9 MB·

by  Lucas Böttcher

Critical behaviors in contagion dynamics
L. B¨ottcher,∗ J. Nagler,† and H. J. Herrmann
ETH Zurich, Wolfgang-Pauli-Strasse 27, CH-8093 Zurich, Switzerland
(Dated: February 1, 2017)
Westudythecriticalbehaviorofageneralcontagionmodelwherenodesareeitheractive(e.g.with
opinion A, or functioning) or inactive (e.g. with opinion B, or damaged). The transitions between
thesetwostatesaredeterminedby(i)spontaneoustransitionsindependentoftheneighborhood,(ii)
transitions induced by neighboring nodes and (iii) spontaneous reverse transitions. The resulting
dynamicsisextremelyrichincludinglimitcyclesandrandomphaseswitching. Wederiveaunifying
mean-field theory. Specifically, we analytically show that the critical behavior of systems whose
dynamics is governed by processes (i-iii) can only exhibit three distinct regimes: (a) uncorrelated
7 spontaneoustransitiondynamics(b)contactprocessdynamicsand(c)cuspcatastrophes. Thisends
1 along-standingdebateontheuniversalityclassesofcomplexcontagiondynamicsinmean-fieldand
0 substantially deepens its mathematical understanding.
2

n
In 1972 Schl¨ogl proposed two models describing au- ior belong to the Ising universality class or not? Here
a
J tocatalytic chemical reactions [1] that are commonly we answer these questions and analytically demonstrate

1 known today as Schl¨ogl’s first and Schl¨ogl’s second that the mean-field critical behavior of the model is re-
3 model henceforth referred to as Schl¨ogl I and Schl¨ogl stricted to only three possible regimes: (a) uncorrelated

II). Schl¨ogl I, also known as contact process (Harris spontaneoustransitiondynamics,(b)contactprocessdy-
]
h 1974), comprisestheimportantcaseofsimple contagion, namicsand(c)cuspcatastrophes. Cuspcatastrophescan
p i.e.thesusceptible-infected-susceptible (SIS)modelwhere display abrupt transitions and hysteresis effects — phe-
- healthy individuals can be infected due to the exposure nomena that can harm the proper functioning of real-
c
o to a single infectious source, eventually leading to the world networked systems since small variations in the
s spread of an epidemic disease [2–5]. In contrast, Schl¨ogl system’scontrolparametersmaycausecatastrophictran-
.
s II that is also known as quadratic contact process [6] re- sitions from a seemingly well-functioning state to global
c
quirescontacttotwosources. LaterstudiesonSchl¨oglII malfunction or severe outages [21–30].
i
s sparkedadebateonitscriticalbehaviorandGrassberger Model The general contagion dynamics is defined in
y
noticed in 1982 that a relation to the Ising universality anetworkwhoseconstituents(i.e.nodes)areregardedas
h
p class ‘would be a most remarkable extension of the uni- either active (e.g. not damaged) or inactive (e.g. failed).
[ versality hypothesis, from models with detailed balance to Three fundamental processes define the transitions be-

models without it’ to conclude that Schl¨ogl II ‘is not an tween these two states [9, 20]: (i) nodes undergo a spon-
1
example of universality between models with and without taneoustransitionA X fromanactive(A)toaninac-
v
→
8 detailed balance’ [7]. tive state (X) in a time interval dt with probability pdt,
8 Closelyrelatedtothisdebate,butmorerecently,agen- (ii) if fewer than or equal to m nearest neighbors of a
9 eralized model of Schl¨ogl II has been proposed where an node are active, the node becomes inactive (Y) due to
8
arbitrarynumberofsourcesisnecessarytoinduceatran- an induced transition, i.e. A Y, with probability rdt
0 →
sition [8]. The study of the model’s mean-field critical and(iii)aspontaneousreversetransitionwithprobability
.
1 behavior led the authors to conjecture that such general qdt if X A or probability q(cid:48)dt if Y A. The inactive
0 → →
failure-recoverydynamicsbelongtotheIsinguniversality states X and Y only differ in their reverse transitions
7
1 class [9]. This model is of particular interest since it not and are equivalent if q = q(cid:48). Process (ii) describes that
: only includes simple contagions but also complex con- a node with degree k can become inactive if its number
v
tagion phenomena such as the diffusion of innovations of inactive neighbors is larger or equal to k m. Similar
i
X [10, 11], political mobilization [12] and viral marketing to threshold models describing complex co−ntagion phe-
r [13] that require social reinforcement, i.e. the connection nomena, the threshold m defines the number of contacts
a
tomultiplesources[14,15]. Themodeldisplaysanintri- to inactive nodes that is necessary to induce a transition
cate and very rich dynamics including hysteresis effects, as defined by process (ii) [14, 31–33]. A low value of m
limit cycles and cusp catastrophes [9, 16–20]. Thus, a corresponds to the situation where many inactive neigh-
unifying mean-field theory of the critical behavior is es- bors are required to sustain spreading. In contrast, for a
sential for a broad range of dynamical systems.. largevalueofmonlyafewinactiveneighborscansustain
However,therelationtocontactprocessdynamicsand the spreading process. Processes (i-iii) are illustrated in
cuspcatastropheshasonlybeenshownforspecificvalues Fig. 1.
of the model’s parameters [20]. But, given the model’s Let a(t) [0,1] denote the total fraction of inactive
∈
parameterregime,canwegenerallypredictthedynamics nodes. Thus, a(t) = u (t)+u (t) with u (t) and
spon ind spon
type? And does the model’s mean-field critical behav- u (t) being the fractions of nodes that are inactive due
ind
2
on the choice of m. It is important to notice that for
more general degree distributions the mean-field critical
behavior still falls into these classes, see Supplementary
Material.
The coupled equations (1) and (2) admit oscillatory
behavior for q(cid:48) >q [20] as a dynamical feature that does
Figure 1. Model. Spontaneous failure (A → X) and spon- not belong to the critical behavior [18]. The equations
taneousrecovery(X →A)occurswithratespandq, respec-
describingthecriticalbehavior,i.e.u˙ =0andu˙ =
spon ind
tively. A node may also fail (become inactive) dependent on
0, can be decoupled by multiplying one of them with
its neighborhood, if too few active nearest neighbors n ≤ m
an appropriate constant excluding limit cycles [18] —
sustainthenode’sactivity(A→Y withrater). Inaddition,
a failed node Y recovers (Y →A) with rate q(cid:48). Active nodes tantamount to setting q =q(cid:48) =1. This yields
(A) are purple while failed ones (X and Y) are grey.
a˙ =f(a,r,p)=rS(a)+p(1 a) a, (3)
− −
with S(a)=(cid:80)m (cid:0) k (cid:1)ak−j(1 a)j+1 =(1 a)E . We
j=0 k−j − − k
tospontaneousandinducedtransitionsrespectively. The use S(a) as shorthand notation for the probability that
total fraction of inactive nodes in the stationary state is an active node is located in a neighborhood that is able
referred to as ast. In accordance with Ref. [20] we de- to induce a transition. Thus, differences in the inactive
rive the mean-field rate equations by assuming a system states X and Y, i.e. different q and q(cid:48), do not influence
with homogeneous degrees in the thermodynamic limit the critical behavior of Eq. (3) but only rescale r and
that exhibits perfect mixing. Here, perfect mixing either p. In the following we analyze the stationary states of
refers to a network of randomly connected nodes with Eq. (3) that will be reached in the long-time limit.
a sufficiently large mean degree or dynamical rewiring Class (a): Uncorrelated spontaneous transitions We
[3, 34]. For the fraction of nodes that spontaneously be- start with the case m = k where the number of active
came inactive we find: nodesnecessarytosustainspreadinghastobesmalleror
equaltothenode’sdegreekaccordingtothedefinitionof
u˙spon =p(1 a) quspon, (1) process (ii). This describes the regime where spreading
− −
occursindependentlyoftheneighborhood’sstatesuchas
where the first term accounts for the fact that active
in exogenously driven adoption dynamics [35, 36],
nodesspontaneouslybecomeinactivewithratep(process
(i)) and the second term corresponds to the spontaneous a˙ =(r+p)(1 a) a (4)
− −
reverse transition with rate q (process (iii)). Eq. (1) is
since E = 1 and S(a) = 1 a. Eq. (4) has only one
k
exactsincethenetworkstructureisnotinfluencingthese −
stationary state, i.e.a (r,p) = (r+p)/(1+r+p), see
st
spontaneous transitions.
Fig. 2 (left).
Induced transitions (process (ii)) can only occur for
Class (b): Contact dynamics Bydefinitionm=k 1
nodes whose number of active neighbors is smaller than −
implies that k 1 or less neighbors of a node have to be
or equal to m. Under the assumption of a perfectly −
active to induce a transition. This case describes a con-
mixed population, the probability that a node of de-
tact process where one inactive neighbor is sufficient to
gree k is located in such a neighborhood is E =
k sustain spreading [3]. As demonstrated in the Appendix,
(cid:80)mj=0(cid:0)k−kj(cid:1)ak−j(1−a)j [9, 20]. The time evolution of wefindforp=0thatthereexistsacriticalrc =k−1sepa-
the fraction of nodes that are inactive due to induced ratinganabsorbingandanactivephase,i.e.a(1)(r,0)=0
transitions is therefore given by: st
as r r and a(2)(r,0) > 0 as r > r . In the limit of
≤ c st c
(cid:88) r r Eq. (3) takes the form:
u˙ind =r fkEk(1 a) q(cid:48)uind, (2) → c
− −
k a˙ =rka(1 a) a. (5)
− −
with f being the degree distribution. The first term Equation (5) describes the mean-field contact process,
k
describes the occurrence of induced transitions (process SIS dynamics or Schl¨ogl I [3, 7]. In the limit of r rc
→
(ii)) with rate r of active nodes in a neighborhood where the order parameter scales as ast(r,0) r rc β with
∼ | − |
the number of active neighbors is smaller than or equal β = 1 and adding the field-like contribution p(1 a) to
−
to m whereas the second term accounts for the sponta- Eq. (5) yields ast(rc,p) p1/δh as p 0 with the field
∼ →
neous reverse transition to an active state with rate q(cid:48) exponent δh = 2. We illustrate in Fig. 2 (center) the
(process (iii)). In order to study the influence of differ- occurrence of only one stable fixed point for p>0. This
entthresholdvaluesmonthemean-fieldcriticalbehavior also results in a smeared-out transition for p>0 instead
of Eqs. (1) and (2), we consider a regular network with of a second-order phase transition for p=0 as shown in
degree k, i.e. the degree distribution fk(cid:48) =δkk(cid:48). We will Fig. 3 (top panel). For p=0 and r >rc one clearly sees
demonstrate below that the model defined by processes that a(1)(r,0) is an unstable but a(2)(r,0) a stable fixed
st st
(i-iii) can only exhibit three different regimes depending point.
3
Figure2. Phase portraits of the general contagion model’s regimes. (toppanel)Thephaseportraitisshownfork=4
anddifferentvaluesofmwitharrowsindicatingthesignofa˙ (rightarrow: a˙ >0,leftarrow: a˙ <0). Blackcirclescorrespondto
stable fixed points and white circles to unstable ones. (left) For m=k there exists only one stable fixed point, whose position
depends on p. We set r=1 and for the grey solid line p=1. (center) For m=k−1 the dynamics resembles the phase space
of a contact process with a stable non-zero fixed point for r > r , where r defines the critical spreading rate below which a
c c
approaches zero if p = 0. For p > 0 and r > r the second-order phase transition gets smeared out. We set r = 1 and for
c
the grey solid line p = 1. (right) If m < k−1, r exceeding r and 0 < p < p (inside the hysteresis region enclosed by the
0 bif
bifurcation lines) it is possible to find two stable fixed points. For p>p the dynamics exhibits only one stable fixed point.
bif
We set r=10 and m=0. For the grey solid line (p<p ) we set p=0.2 and for silver solid line (p>p ) we used p=0.6.
bif bif
(bottom panel) The time evolution for different parameters in a regular random graph with N =256 nodes and k=4. (left)
Form=k=4,p=0.01,r=1.0thedynamicsgrowsuntilastationarystateisreached. (center)Form=k−1=3,p=0.01,
r=1.0wefindthetypicallogisticgrowthpattern. (right)Form<k−1=1,p=0.24,r=10weencounterphase-switching.
Class (c): Cusp catastrophes For some values of m, and
wefindametastableregionasillustratedinFig.3(right).
S(cid:48)(a )a S(a )
Insidethishysteresisregiontwostablefixedpointscoex- p (k,m)= 0 0− 0 . (7)
0 S(a )+S(cid:48)(a )(1 a )
ist. Phase-switchingisobservedwhenfluctuationsinsys- 0 0 0
−
tems of finite size push the dynamics close to the unsta-
We illustrate the influence of different values of m on
ble fixed point, cf. Fig. 2 (right). Between the switching
(r ,p ) and on the extent of the hysteresis area in the
eventsthedynamicsremainsinoneofthetwophasesfor 0 0
Appendix. Studying Eq. (3) in the vicinity of (a ,r ,p ),
some time. The waiting times thus depend on the fluc- 0 0 0
i.e. setting a = a +a˜, r = r +r˜, p = p +p˜, yields for
tuation strength and the distance from one phase to the 0 0 0
the Taylor expansion (omitted tilde)
unstablestateinthephaseportrait,cf.Fig.2(right). For
m < k 1 where either two or more inactive neighbors r
are nec−essary to induce a transition, we now show that a˙ =rS(a0)+p(1−a0)+a(rS(cid:48)(a0)−p)+ 60S(cid:48)(cid:48)(cid:48)(a0)a3. (8)
thecorrespondingmetastableregionsalwaysexistdueto
We thus find by setting p or r to zero respectively and
the relation to cusp catastrophes [17]. For a detailed an-
solving for the fixed point of Eq. (8):
alytical treatment we refer to the Appendix. The cusp
point where the two bifurcation lines intersect (cf. Fig. 3
(cid:20) 6S(a ) (cid:21)1/3 (cid:16) (cid:17)
(right))isgivenbya0(k,m)=(k 1 m)/(k+1)together a(r)= 0 r1/3+ r2/3 , (9)
with the corresponding control p−ara−meters −r0S(cid:48)(cid:48)(cid:48)(a0) O
1 (cid:20) 6(1 a )(cid:21)1/3 (cid:16) (cid:17)
r (k,m)= , (6) a(p)= − 0 p1/3+ p2/3 . (10)
0 S(a )+S(cid:48)(a )(1 a ) −r S(cid:48)(cid:48)(cid:48)(a ) O
0 0 0 0 0
−
4
tivenodesareneededtotriggerinducednode-transitions.
This scenario typically implies dramatic and uncontrol-
lable global transitions in the network for many sys-
tems involving complex contagion dynamics. One could
naively expect that it could be beneficial for failure con-
trol to design systems such that a component only fails
if many of its neighbors already failed, i.e. delaying the
failuredynamics. Ourresultssuggest, however, thatthis
delaying procedure might facilitate uncontrollable tran-
sitions, hence achieving exactly the opposite as initially
intended. This result agrees well with previous findings
on delaying procedures which have been applied to a
SIS model [38, 39]. For low spatial dimensions or highly
structured networks, the assumptions of perfect mixing
or independent node-to-node interactions are not guar-
anteed. Still mean-field approximations qualitatively de-
scribe a given dynamics [3, 4, 40], see given examples in
the Appendix.
Futureworkshouldestablishthebehavioroftransients
asafunctionofthresholdparametermandthetopology
of the network. It has been demonstrated that opinions
aswellascoinfectionsmayspreadfasterinclusterednet-
works compared to random ones [41, 42]. This links our
resulttothemultipleexposureconditionincomplexcon-
Figure 3. Phase diagrams of contact dynamics and
tagion phenomena.
cusp catastrophes. (top panel) In the case of m = k−1,
In the study of collective behaviors, such as the adop-
thephasespacecorrespondstotheoneofSchlo¨glIorcontact
dynamics. For vanishing p, a second-order phase transition tion of innovations, the distinction between exogenous
from an absorbing to an active phase occurs at r . For a and endogenous factors is of great interest but often
c
positive value of p, the transition gets smeared out [3, 37]. solely based on a contact process-like adoption model
(bottompanel)Forallm<k−1,cuspcatastrophesdefinethe [35, 36]. Our results suggest studying these processes
phase space. Two in the thermodynamic limit stable states
within our more general framework that incorporates
coexistinsidethehysteresisregion(II)thatissurroundedby
contactprocess-likeadoptionasonespecialcaseandcan
bifurcationlines(blacksolidlines)mergingatthecusppoint
account for spreading that relies on multiple contacts.
(r ,p ). Outsidethehysteresisregion(I)onlyonestablestate
0 0
exists. The critical point r is defined as the transition point
c
for vanishing p.
Table I. Examples of models and processes that are related
to the classes (a-c).
(a) m=k (b) m=k−1 (c) m<k−1
In previous work, the critical behavior at the cusp point
Exogenous fac- Schl¨ogl I [1, 7] Schlo¨gl II [1, 7]∗
of a regular random network (k = 10 and m = 4) has
tors influencing Contact process Quadratic con-
beenconjecturedtobelongtotheIsinguniversalityclass
adoption of inno- [2, 3, 37]∗ tact process [6]∗
although by definition the dynamics corresponds to a vations [36]∗ SIS model [4, 5]∗ General contact
general contact process [20]. Social response Reggeon field process [8]∗
Final remarks We find that the critical behavior of to exogenous theory [43]∗ Behavioral adop-
thegeneralcontagionmodelasformulatedinEq.(3)does factors [35]† Directed percola- tion [41]†
not belong to the Ising universality class but to exactly tion [44]∗ Thresholdmodels
Bass model of complex con-
three regimes. The first regime, m = k, corresponds to
[36, 45]∗ tagions [11, 13–
purely spontaneous failure and recovery dynamics. For
15, 31–33]∗† or
m = k 1 the model recovers the critical behavior of
coordination
−
the contact process. A cusp catastrophe is found for all games [46]†
m < k 1 with the typical critical behavior at the cusp ∗exactmean-fieldcorrespondence
point (−Eqs. (9) and (10)). This sheds analytical insight †phenomenologicalcorrespondence
into a broad range of spreading processes that are deter-
mined by the network’s connectivity k and the threshold
We acknowledge financial support from the ETH Risk
parameter m, cf. examples in Tab. I.
Center (grant RC SP 08-15) and ERC Advanced grant
We have demonstrated that the phase diagram cor-
number FP7-319968 FlowCCS.
responds to a cusp catastrophe, when two or more inac-
5
Appendix
Critical behavior in a regular network
This section is dedicated to the analytical derivation of the mean-field critical behavior of Eqs. (1) and (2) in the
main manuscript. Thus, we have to focus on the corresponding stationary states, i.e. u˙ = 0 and u˙ = 0. Both
spon ind
equations can be decoupled by multiplying one of the two with an appropriate constant and adding them up — the
same can be achieved by setting q =q(cid:48) =1 yielding:
a˙ =f(a,r,p)=rS(a)+p(1 a) a, (11)
− −
where S(a)=(cid:80)m (cid:0) k (cid:1)ak−j(1 a)j+1.
j=0 k−j −
Class (a): Uncorrelated spontaneous transitions In the case of m = k, i.e. spontaneous transitions that occur
independentoftheneighborhood,thegoverningrateequationisa˙ =(r+p)(1 a) asince(cid:80)k (cid:0) k (cid:1)ak−j(1 a)j =
− − j=0 k−j −
1. Therefore, the corresponding stationary state is given by a (r,p)=(r+p)/(1+r+p).
st
Class (b): Contact dynamics For m=k 1 we first set p=0 and find
−
 
k (cid:18) (cid:19)
(cid:88) k
a˙ =r ( 1)j+1aj(1 a) a (12)
j − − −
j=1
since (cid:80)k−1(cid:0) k (cid:1)ak−j(1 a)j = 1 (1 a)k. We directly observe that there is a fixed point a(1)(r,0) = 0 at
j=0 k−j − − − st
the origin of the phase portrait, cf. Fig. 2 (center). Furthermore, the function f(a,r,p) of Eq. (11) is given by
(cid:104) (cid:105)
f(a,r) = r (1 a) (1 a)k+1 a and we find f(cid:48)(0,r) = rk 1. Consequently, f(cid:48)(0,r) > 0 if r > k−1 = r and
c
− − − − −
there exists one maximum a (r,k) = 1 (cid:104) 1 (cid:0)r−1+1(cid:1)(cid:105)1/k > 0 if r > r (f(cid:48)(cid:48)(a ) < 0). Thus, a(1)(r,0) = 0 is
max − k+1 c max st
the stable fixed point if r r . For r > r , a(1)(r,0) = 0 is unstable with f(cid:48)(0,r) > 0 and a maximum a (r,k)
≤ c c st max
exists obeying 0 < a (r,k) < 1 (f(a ) > 0). As a consequence of the latter fact and f(1,r) = 1, a second
max max
−
second fixed point a(2)(r,0) (a ,1) exists, that is stable for r >r as illustrated in Fig. 2. In the vicinity of the
st ∈ max c
transition point (r r ), we only consider the dominant linear term in the sum of Eq. (12) and obtain:
c
→
a˙ =rka(1 a) a. (13)
− −
As r r the order parameter scales as a (r,0) r r β with β = 1 and assuming a non-zero field like term
c st c
→ ∼ | − |
p(1 a) yields ast(rc,p) p1/δh as p 0 with δh =2.
− ∼ →
Class(c): Cuspcatastrophes Ifm<k 1,thecriticalbehaviorisdescribedbyacuspcatastrophe. Theequilibrium
−
point a of our dynamical system a˙ = f(a,r,p) is said to correspond to a cusp catastrophe [17] for the parameters
0
r ,p if it satisfies the following conditions [47]:
0 0
(i) The point a is a non-hyperbolic equilibrium.
0
(ii) The quadratic term of the function f vanishes but not the cubic one, i.e.
∂2f(a,r,p)(cid:12)(cid:12) ∂3f(a,r,p)(cid:12)(cid:12)
(cid:12) =0 and (cid:12) =0.
∂a2 (cid:12) ∂a3 (cid:12) (cid:54)
(a,r,p)=(a0,r0,p0) (a,r,p)=(a0,r0,p0)
(iii) The two vectors
(cid:32) ∂f (cid:33) (cid:32) ∂2f (cid:33)
v= ∂r and w= ∂r∂a
∂f ∂2f
∂p ∂p∂a
are linearly independent.
We start with condition (ii) and find that the vanishing second derivative at a implies:
0
(cid:88)m (cid:18) k (cid:19)(1 a )ja−j(cid:104)j(j+1)+(a 1)(1+a +2j)k+(a 1)2k2(cid:105)
k j − 0 0 0− 0 0−
j=0 − (14)
(cid:18) (cid:19)
k
= (1 a )m+1a−m(1+m)[k(1 a ) 1 a m]=0.
k 1 m − 0 0 − 0 − − 0−
− −
6
The cases a =0 or a =1 do not correspond to an equilibrium point of Eq. 3. We therefore focus on the remaining
0 0
solution a (k,m)=(k 1 m)/(k+1). The third derivative at this point is given by:
0
− −
∂3f(a,r,p) (1+k)3(1+m)(cid:16)k−k1+−1m(cid:17)k−m(cid:16)21++mk(cid:17)m(cid:0)k−1k−m(cid:1)
= <0,
∂a3 |(a,r,p)=(a0,r0,p0) − (k 1 m)2
− −
and clearly non-vanishing since m < k 1. In order to satisfy condition (i), i.e. f(a ,r ,p ) = f(cid:48)(a ,r ,p ) = 0, we
0 0 0 0 0 0
−
find:
1
r (k,m)= , (15)
0 S(a )+S(cid:48)(a )(1 a )
0 0 0
−
S(cid:48)(a )a S(a )
p (k,m)= 0 0− 0 . (16)
0 S(a )+S(cid:48)(a )(1 a )
0 0 0
−
We need to demonstrate that r and p are well behaved, i.e. 0 < r < and 0 < p < . Therefore, we need to
0 0 0 0
∞ ∞
show (1) S(cid:48)(a )a S(a )>0 and (2) S(a )+S(cid:48)(a )(1 a )>0. The first condition leads to
0 0 0 0 0 0
− −
m (cid:18) (cid:19) (cid:20) (cid:21)
S(cid:48)(a )a =(cid:88) k ak−j(1 a )j+1 (k j) a0 (j+1) (17)
0 0 k j 0 − 0 − − 1 a
0
j=0 − −
m (cid:18) (cid:19)
>(cid:88) k ak−j(1 a )j+1 =S(a ), (18)
k j 0 − 0 0
j=0 −
whatisfulfilledif(k−j)−1−a0a0 (j+1)>1. Sincea0 =(k−1−m)/(k+1),weobtainthecondition(k+1)(m+1−j)>
m+2 which is always satisfied for m 0,...,k 2 and j 0,...,m . For condition (2) it is sufficient to show
∈ { − } ∈ { }
that S(cid:48)(a ) > 0 — a consequence of the previous statement. The remaining condition (iii) for v,w R2 is clearly
0
∈
met since
(cid:18) S(a ) (cid:19) (cid:18)S(cid:48)(a )(cid:19)
v= 0 and w= 0
1 a 1
0
− −
are linearly independent. We have thus shown that Eq. (11) has an equilibrium point a for the parameters r , p
0 0 0
that corresponds to a cusp catastrophe if m < k 1. We study Eq. (11) in the vicinity of (a ,r ,p ) by setting
0 0 0
−
a=a +a˜, r =r +r˜, p=p +p˜what yields for its Taylor expansion (omitted tilde sign) [47]:
0 0 0
1
a˙ =f(a ,r ,p )+v λ+w λa+ r S(cid:48)(cid:48)(cid:48)(a )a3+h.o.t., (19)
0 0 0 · · 3! 0 0
where λ = (r,p) and h.o.t. describes the higher-order terms in a, r and p. Neglecting the terms of higher-order, the
latter equation can be also written as:
r
a˙ =rS(a )+p(1 a )+a(rS(cid:48)(a ) p)+ 0S(cid:48)(cid:48)(cid:48)(a )a3. (20)
0 − 0 0 − 6 0
We find the critical behavior in the vicinity of the cusp point by setting p or r to zero respectively and solving for the
fixed point of Eq. (20):
(cid:20) 6S(a ) (cid:21)1/3 (cid:16) (cid:17)
a(r)= 0 r1/3+ r2/3 , (21)
−r S(cid:48)(cid:48)(cid:48)(a ) O
0 0
(cid:20) 6(1 a )(cid:21)1/3 (cid:16) (cid:17)
a(p)= − 0 p1/3+ p2/3 . (22)
−r S(cid:48)(cid:48)(cid:48)(a ) O
0 0
7
3 5
A
m=0 H|r<rc 4
)
2 m
A (r r )p
q (c3 H|r<rc ∼ c− 0 0
/ m=1 <r
p
r 2
1 m=2 H|
. A
.. 1
0 0
0 5 10 15 20 25 30 0 1 2 3 4 5 6 7 8
r/q m
0
Figure 4. Extent and position of hysteresis regions. (left) Hysteresis regions for k = 10 and different m. (right)
Hysteresis areas for k = 10 as function of m. The black line is the approximation A (r ,p ) ≈ 1(r −r )p where
H|r<rc 0 0 3 c 0 0
the m-dependence is adopted from the ones of r and p in Eqs.(15) and (16). The parameters r and p are varied whereas
0 0
q=q(cid:48) =1.
Parameter dependence of the hysteresis area
Position and extent of the hysteresis region impact the controllability of failure-recovery dynamics defined by
processes (i-iii) [20]. Here, we study the influence of different values of the threshold m keeping k fixed (Fig. S4). For
k =10 and m 0,...,8 the grey dots in Fig. S4 (left) are the cusp points described by Eqs. (15) and (16) and the
∈{ }
following curve describes their position in phase space:
(cid:18) (cid:19)
r (k,m)
γ(k,m)= 0 . (23)
p (k,m)
0
The hysteresis area decreases with m.
Thismeans,themorecontactstoinactivenodesnecessary,thelargerthemetastabledomainandthelargerthreshold
values, the control parameters r and p for which such behavior occurs. We illustrate this effect for hysteresis areas
A restricted to r <r in Fig. S4 (right). The black dots represent the actual values of the corresponding area
H|r<rc c
and the black line is an approximation based on the assumption that A (r ,p ) (r r )p . The value of r ,
H|r<rc 0 0 ∼ c− 0 0 c
as defined in Fig. 3 (right), is determined by Eq. (11) in the limit p 0.
→
Critical behavior for general degree distributions
Any degree distribution leads to a characteristic polynomial of maximum degree three (Eq. (8) in the manuscript).
This is, in general, a consequence for two-parameter bifurcations of smooth dynamical systems with exactly two
control parameters [17]. Therefore, our mathematical framework is valid for any choice of f . However, as expected,
k
the mean-field theory is not exact for, e.g., networks with low degree, or networks with a broad degree distribution,
but nevertheless, the maximal degree of the characteristic polynomial necessarily remains invariant under the specific
form of the degree distribution and therefore constraints the critical behavior to fall into classes (a-c). Below we
describe the treatment of general degree distributions in our mean-field theory and also discuss specific examples —
bearing in mind that the arguments apply for a perfectly mixed mean-field configuration.
The three classes (a-c) we discussed above in Sec. Critical behavior in a regular network are recurring even for
more general degree distributions and are useful to qualitatively understand the mean-field dynamics in more general
situations. For values of m that depend on the actual node degree, an analogous treatment is possible. In our
mean-field theory a general degree distribution different from a regular one (f =δ ) implies multiple terms in the
k(cid:48) kk(cid:48)
(cid:80)
sum f . A degree distribution with finite cut-off k is of the form f = P δ +P δ + +P δ
k(cid:48) k(cid:48) (cid:80) n k(cid:48) k1 k1k(cid:48) k2 k2k(cid:48) ··· kn knk(cid:48)
with normalization P =1. Here P denotes the probability of degree k . We now define S (a)=
(cid:80)mj=0(cid:0)kik−ij(cid:1)aki−j(1−ak(cid:48))∈j{+k11,=k2,(..1.,k−n}a)Ekk(cid:48)i. Instead ofkEiq. 11, a general degree distribution ais specified above wkiill lead
8
1.0
10% spont. and 90% contact
0.8
50% spont. and 50% contact 1.0
)0.6 90% spont. and 10% contact p) (r ,p )
(r Approximation for small ast(r) r, 0 0
t ( 0.5
s0.4 t
a s
a
0.2
0.0
0.3
0.0
0.2
0.0 0.2 0.4 0.6 0.8 1.0 p 0.1 10
8
r
0.0 6 r
Figure 5. Examples of mean-field critical behavior for general degree distributions. (left) An illustration of the
growth of the order parameter a (r) for different mixtures of spontaneous and contact dynamics. The bold grey lines are an
st
approximation for small values of a (r), cf. Eq. (26) with k = 4. (right) The phase-space of a mixture of spontaneous and
st i
cusp dynamics according to Eq. (27) with P =0.95 and k =4, m=0. The cusp point, as defined by Eqs. (28) and (29), is
ki i
indicated by (r ,p ), The black solid lines enclose the hysteresis region.
0 0
to:
a˙ =P [rS (a)+p(1 a) a]+P [rS (a)+p(1 a) a]+
k1 k1 − − k2 k2 − − (24)
P [rS (a)+p(1 a) a]+ +P [rS (a)+p(1 a) a].
k3 k3 − − ··· kn kn − −
We study the steady states of Eq. (24) to learn about the critical behavior of a general degree distribution. We
are now dealing with a sum of terms P [rS (a)+p(1 a) a] that are similar to the ones of a regular degree
ki ki − −
distribution. Therefore, each term corresponds to cusp, contact or spontaneous dynamics, respectively and the more
general behavior can be reduced to the ones of a regular network. Depending on the difference m k , different
i
−
phase-space configurations are possible. We describe some possible cases below to illustrate that the three classes of
critical behavior of a regular network are recurring and essential to understand this more general configuration.
For m k , the overall dynamics is necessarily spontaneous (as m k for all i).
max i
• ≥ ≥
Spontaneous critical dynamics is observed as well if, for example, the i-th term leads to contact dynamics
•
(m=k 1) and the remaining terms lead to spontaneous ones. Since m is fixed, only one k is equal to m 1
i i
−(cid:80) −
and with P =1 P (see Sec. Critical behavior in a regular network):
k(cid:48)∈{k1,k2,...,kn}\{ki} k(cid:48) − ki
0=P (cid:2)r(cid:0)1 (1 a)ki(cid:1)(1 a)+p(1 a) a(cid:3)+(1 P )[(r+p)(1 a) a]. (25)
ki − − − − − − ki − −
A second-order phase transition characterizing the critical behavior of the contact process exists no longer due
to spontaneous transitions. For vanishing p, the critical behavior of Eq. (25), i.e. in this case an expansion for
small a, yields cf. Fig. S5:
(1 P )r
a (r,p=0,k ,P )= − ki +h.o.t., (26)
st i ki (1 P )(1+r)+P (1 k r)
− ki ki − i
where h.o.t. denotes the higher order terms in the occurring variables. We note that Eq. (26) corresponds to
the steady states described by Eq. (4) for P = 0 (degrees k are absent), i.e. a (r,p = 0) = r/(1+r) for
ki i st
spontaneous dynamics in a regular network.
Another scenario is the occurrence of a cusp catastrophe for the i-th term (m < k 1) and spontaneous
i
• −
dynamics otherwise. The spontaneous terms just lead to additional constant and linear contributions in the
polynomial describing the steady states:
0=P [rS (a)+p(1 a) a]+(1 P )[(r+p)(1 a) a]. (27)
ki ki − − − ki − −
The cusp point a (k,m)=(k 1 m)/(k+1) for regular degree distributions has been derived via a vanishing
0
− −
second order derivative, cf. Sec. Critical behavior in a regular network. Therefore, the additional constant
9
and linear terms are not affecting the existence of the cusp point in Eq. (27). As for the regular graph, the
non-hyperbolic equilibrium condition yields:
1
r (k,m,k ,P )= , (28)
0 i ki P (cid:2)S (a )+S(cid:48) (a )(1 a )(cid:3)
ki ki 0 ki 0 − 0
P 1+P (cid:2)S(cid:48) (a )a S (a )(cid:3)
p (k,m,k ,P )= ki − ki ki 0 0− ki 0 . (29)
0 i ki P (cid:2)S (a )+S(cid:48) (a )(1 a )(cid:3)
ki ki 0 ki 0 − 0
We note that for P = 1 (degrees k occur P-almost surely) Eqs. (28) and (29) correspond to Eqs. (15) and
ki i
(16) for a regular graph. However, unlike Eq. (16), Eq. (29) might be a negative, unphysical rate p . Therefore,
0
the proportion of spontaneous terms determines if the cusp point exists for physical, i.e. positive rates r and p
inthephase-spaceornot. Consequently, onemightnotobserveahysteresisregionifspontaneousdynamicshas
too much influence. We illustrate a phase-space configuration for a mixture of spontaneous and cusp dynamics
in Fig. 5 (right) with P =0.95 (degrees k occur with probability 0.95) and m=0.
ki i
1.2
m=4
1.0
0.8
m=3
0.6
˙a
0.4
0.2 m=2
0.0
0.2
− 0.0 0.2 0.4 0.6 0.8 1.0
a
Figure6. Phase-portraitconfigurationsofapower-lawdegreedistribution. Themean-fieldphase-portraitforadegree
distribution P(k)∝k−1 with k =4, k =100, r =10 and different values of m. For m=2 and m=4 we set p=0.05
min max
and p = 0 for m = 3 to illustrate the existence of the unstable fixed point at the origin. Arrows are indicating the sign of a˙
(right arrow: a˙ >0, left arrow: a˙ <0). Black circles correspond to stable fixed points and white circles to unstable ones.
In Fig. 6, we study the influence of a power-law degree distribution P(k) k−1 with k = 4 and k = 100
min max
∝
on the mean-field critical behavior as described by Eq. (24). We qualitatively recover the three classes describing
spontaneous, contact and cusp dynamics. In particular, for m = 2 we find a polynomial that corresponds to cusp
dynamics and for m = 3 the phase-portrait configuration corresponds to contact dynamics. In the case of m = 4,
the dynamics is spontaneous although the phase-portrait is not linear as in Fig. (2) (in the manuscript) for a regular
network. This is due to the higher order contributions but still there exists only one fixed point and the dynamics is
spontaneously driven.
10
Stochastic simulations
In this section, we present stochastic simulations of spreading processes on real networks compared to the corre-
spondingmean-fieldpredictions. InFig.S7,weconsideredafriendshipnetwork[48]with2539nodesand20910edges
sinceforexampleepidemicsoropinionsspreadthroughsuchnetworks. Thecorrespondingdegreedistributionisillus-
trated inthe inset ofFig. S7. Based onthis degree distribution, we numerically computeda mean-field solution(grey
solid lines). In Fig. S8, we present simulations of the model’s dynamics in a regular random and in an Erd˝os-R´enyi
network.
1.0
m=5
0.8
0.12
0.6
F
) D0.06 MFT
t
( P
a0.4 Simulation
0.00
0 10 20 30
0.2 Degree
m=0
0.0
0 10 20 30 40 50
t
Figure 7. Dynamics in a social network. The time-evolution of the fraction of inactive nodes a(t) for r = 10, p = 0.1,
q = q(cid:48) = 1, m = 0 and m = 5 (red solid lines) in a friendship network with 2539 nodes and 20910 edges — the inset shows
its degree distribution. The grey solid lines are the numerical solution of Eqs. (1) and (2) with the corresponding degree
distribution that can be approximated by a truncated normal distribution.
1.0 1.0
stable stable
0.8 0.8 0.10
F
0.6 unstable 0.6 D0.05
) ) P unstable
t t
a(0.4 Simulation a(0.4 0.00
0 20 40
MFT
Degree MFT
0.2 0.2 Simulation
stable stable
0.0 0.0
0 2 4 6 8 10 0 2 4 6 8 10
t t
Figure 8. Stochastic simulation and mean-field theory. The time-evolution of the fraction of inactive nodes a(t) for
different initial conditions a(t = 0) in a (left) regular random network with k = 20 and (right) an Erdo˝s-R´enyi network with
averagedegree(cid:104)k(cid:105)=20—itsdegreedistributionisshownasaninset. Bothnetworkshave100.000nodes. Theredsolidlines
representthesolutionsofthestochasticsimulation. Thegreysolidlinesarethecorrespondingmean-fieldsolutions. Twostable
states indicate that the dynamics is located in the hysteresis region. The parameters are r=10, p=0.1, q=q(cid:48) =1, m=5.