Critical behaviors in contagion dynamics

TL;DR

This paper analyzes the critical behavior of a general contagion model with nodes switching states due to spontaneous, neighbor-induced, and reverse transitions, revealing three universal regimes and deepening the mathematical understanding of contagion dynamics.

Contribution

It provides a unifying mean-field theory that classifies the critical regimes of complex contagion dynamics, resolving a long-standing debate.

Findings

Identifies three universal regimes: uncorrelated, contact process, and cusp catastrophes.

Derives a mean-field theory that unifies different contagion behaviors.

Deepens the mathematical understanding of contagion phase transitions.

Abstract

We study the critical behavior of a general contagion model where nodes are either active (e.g. with opinion A, or functioning) or inactive (e.g. with opinion B, or damaged). The transitions between these two states are determined by (i) spontaneous transitions independent of the neighborhood, (ii) transitions induced by neighboring nodes and (iii) spontaneous reverse transitions. The resulting dynamics is extremely rich including limit cycles and random phase switching. We derive a unifying 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 spontaneous transition dynamics (b) contact process dynamics and (c) cusp catastrophes. This ends a long-standing debate on the universality classes of complex contagion dynamics in mean-field and substantially deepens its mathematical understanding.

Tables

Table 1. Table ​​ 1: Examples of models and processes that are related to the classes (a-c).

(a) $m=k$	(b) $m=k-1$	(c)
Exogenous factors influencing adoption of innovations Ghanbarnejad et al. (2014)∗ Social response to exogenous factors Crane and Sornette (2008)†	Schlögl I Schlögl (1972); Grassberger (1982) Contact process Harris (1974); Marro and Dickman (2005); Henkel et al. (2008)∗ SIS model Keeling and Rohani (2008); Pastor-Satorras et al. (2015)∗ Reggeon field theory Grassberger and de la Torre (1979)∗ Directed percolation Cardy and Sugar (1980)∗ Bass model Bass (1969); Ghanbarnejad et al. (2014)∗	Schlögl II Schlögl (1972); Grassberger (1982)∗ Quadratic contact process Durrett ∗ General contact process Tomé and de Oliveira (2015)∗ Behavioral adoption Centola (2010)† Threshold models of complex contagions Granovetter (1978); Watts (2002); Leskovec et al. (2007); Centola and Macy (2007); López-Pintado (2008); Rogers (2010); Gleeson (2013)∗† or coordination games Easley and Kleinberg (2010)†
$\ast$exact mean-field correspondence
$†$phenomenological correspondence

Equations

$${\dot{u}}_{\text{spon}}=p\left(1-a\right)-q{u}_{\text{spon}},$$

$${\dot{u}}_{\text{ind}}=r\sum _k{f}_k{E}_k\left(1-a\right)-q^{\prime }{u}_{\text{ind}},$$

$$\dot{a}=f(a,r,p)=rS(a)+p\left(1-a\right)-a,$$

$$\dot{a}=(r+p)\left(1-a\right)-a$$

$$\dot{a}=rka\left(1-a\right)-a.$$

$$r_0(k,m)=\frac{1}{S(a_0)+S^{\prime }(a_0)(1-a_0)},$$

$$p_0(k,m)=\frac{S^{\prime }(a_0)a_0-S(a_0)}{S(a_0)+S^{\prime }(a_0)(1-a_0)}.$$

$$\dot{a}=rS(a_0)+p(1-a_0)+a(r{S}^{\prime }(a_0)-p)+\frac{r_0}{6}{S}^{\prime \prime \prime }(a_0)a^3.$$

$$a(r)={\left[-\frac{6S(a_0)}{r_0{S}^{\prime \prime \prime }(a_0)}\right]}^{1/3}{r}^{1/3}+\mathcal{O}\left(r^{2/3}\right),$$

$$a(p)={\left[-\frac{6\left(1-a_0\right)}{r_0{S}^{\prime \prime \prime }(a_0)}\right]}^{1/3}{p}^{1/3}+\mathcal{O}\left(p^{2/3}\right).$$

$$\dot{a}=f(a,r,p)=rS(a)+p\left(1-a\right)-a,$$

$$\dot{a}=r\left[\sum _{j=1}^k\left(\genfrac{}{}{0ex}{}{k}{j}\right)(-1{)}^{j+1}{a}^j\right]\left(1-a\right)-a$$

$$\dot{a}=rka\left(1-a\right)-a.$$

$$\frac{\partial^{2}f(a,r,p)}{\partial a^{2}}\biggr{\rvert}_{(a,r,p)=(a_{0},r_{0},p_{0})}=0~{}\text{and}~{}\frac{\partial^{3}f(a,r,p)}{\partial a^{3}}\biggr{\rvert}_{(a,r,p)=(a_{0},r_{0},p_{0})}\neq 0.$$

$$\mathbf{v}=\left(\begin{array}[]{c}\frac{\partial f}{\partial r}\\ \frac{\partial f}{\partial p}\\ \end{array}\right)~{}\text{and}~{}\mathbf{w}=\left(\begin{array}[]{c}\frac{\partial^{2}f}{\partial r\partial a}\\ \frac{\partial^{2}f}{\partial p\partial a}\\ \end{array}\right)$$

$$\begin{array}{rl}& \sum _{j=0}^m\left(\genfrac{}{}{0ex}{}{k}{k-j}\right){\left(1-a_0\right)}^j{a}_0^{-j}\left[j\left(j+1\right)+\left(a_0-1\right)\left(1+a_0+2j\right)k+{\left(a_0-1\right)}^2{k}^2\right]\\ & =\left(\genfrac{}{}{0ex}{}{k}{k-1-m}\right){\left(1-a_0\right)}^{m+1}{a}_0^{-m}\left(1+m\right)\left[k(1-a_0)-1-a_0-m\right]=0.\end{array}$$

$$\frac{{\partial }^{3}f(a,r,p)}{\partial a^3}{|}_{(a,r,p)=(a_0,r_0,p_0)}=-\frac{(1+k{)}^3(1+m){\left(\frac{k-1-m}{k+1}\right)}^{k-m}{\left(\frac{2+m}{1+k}\right)}^m(\genfrac{}{}{0ex}{}{k}{k-1-m})}{(k-1-m{)}^2}<0,$$

$$r_0(k,m)=\frac{1}{S(a_0)+S^{\prime }(a_0)(1-a_0)},$$

$$p_0(k,m)=\frac{S^{\prime }(a_0)a_0-S(a_0)}{S(a_0)+S^{\prime }(a_0)(1-a_0)}.$$

$$S^{\prime }(a_0)a_0=\sum _{j=0}^m\left(\genfrac{}{}{0ex}{}{k}{k-j}\right)a_0^{k-j}{\left(1-a_0\right)}^{j+1}\left[\left(k-j\right)-\frac{a_0}{1-a_0}\left(j+1\right)\right]$$

$$>\sum _{j=0}^m\left(\genfrac{}{}{0ex}{}{k}{k-j}\right)a_0^{k-j}{\left(1-a_0\right)}^{j+1}=S(a_0),$$

$$\mathbf{v}=\left(\begin{array}[]{c}S(a_{0})\\ 1-a_{0}\\ \end{array}\right)~{}\text{and}~{}\mathbf{w}=\left(\begin{array}[]{c}S^{\prime}(a_{0})\\ -1\\ \end{array}\right)$$

$$\dot{a}=f(a_0,r_0,p_0)+\mathbf{v}\cdot \lambda +\mathbf{w}\cdot \lambda a+\frac{1}{3!}{r}_0{S}^{\prime \prime \prime }(a_0)a^3+\text{h.o.t.},$$

$$\dot{a}=rS(a_0)+p(1-a_0)+a(r{S}^{\prime }(a_0)-p)+\frac{r_0}{6}{S}^{\prime \prime \prime }(a_0)a^3.$$

$$a(r)={\left[-\frac{6S(a_0)}{r_0{S}^{\prime \prime \prime }(a_0)}\right]}^{1/3}{r}^{1/3}+\mathcal{O}\left(r^{2/3}\right),$$

$$a(p)={\left[-\frac{6\left(1-a_0\right)}{r_0{S}^{\prime \prime \prime }(a_0)}\right]}^{1/3}{p}^{1/3}+\mathcal{O}\left(p^{2/3}\right).$$

$$\gamma(k,m)=\left(\begin{array}[]{c}r_{0}(k,m)\\ p_{0}(k,m)\\ \end{array}\right).$$

$$\dot{a}=$$

$$P_{k_3}\left[r{S}_{k_3}(a)+p(1-a)-a\right]+\cdots +P_{k_n}\left[r{S}_{k_n}(a)+p(1-a)-a\right].$$

$$0=P_{k_i}\left[r\left(1-(1-a{)}^{k_i}\right)(1-a)+p(1-a)-a\right]+(1-P_{k_i})\left[(r+p)(1-a)-a\right].$$

$$a_{st}\left(r,p=0,k_i,P_{k_i}\right)=\frac{(1-P_{k_i})r}{(1-P_{k_i})(1+r)+P_{k_i}(1-k_{i}r)}+\text{h.o.t.},$$

Full text

In 1972 Schlögl proposed two models describing autocatalytic chemical reactions Schlögl (1972) that are commonly known today as Schlögl’s first and Schlögl’s second model henceforth referred to as Schlögl I and Schlögl II). Schlögl I, also known as contact process (Harris 1974), comprises the important case of simple contagion, i.e. the susceptible-infected-susceptible (SIS) model where healthy individuals can be infected due to the exposure to a single infectious source, eventually leading to the spread of an epidemic disease Harris (1974); Marro and Dickman (2005); Keeling and Rohani (2008); Pastor-Satorras et al. (2015). In contrast, Schlögl II that is also known as quadratic contact process Durrett requires contact to two sources. Later studies on Schlögl II sparked a debate on its critical behavior and Grassberger noticed in 1982 that a relation to the Ising universality class ‘would be a most remarkable extension of the universality hypothesis, from models with detailed balance to models without it’ to conclude that Schlögl II ‘is not an example of universality between models with and without detailed balance’ Grassberger (1982).

Closely related to this debate, but more recently, a generalized model of Schlögl II has been proposed where an arbitrary number of sources is necessary to induce a transition Tomé and de Oliveira (2015). The study of the model’s mean-field critical behavior led the authors to conjecture that such general failure-recovery dynamics belong to the Ising universality class Majdandzic et al. (2014). This model is of particular interest since it not only includes simple contagions but also complex contagion phenomena such as the diffusion of innovations Coleman et al. (1957); Rogers (2010), political mobilization Chwe (1999) and viral marketing Leskovec et al. (2007) that require social reinforcement, i.e. the connection to multiple sources Granovetter (1978); Centola and Macy (2007). The model displays an intricate and very rich dynamics including hysteresis effects, limit cycles and cusp catastrophes Ludwig et al. (1978); Zeeman (1979); Strogatz (2014); Majdandzic et al. (2014); Valdez et al. (2016); Böttcher et al. (2016a). Thus, a unifying mean-field theory of the critical behavior is essential for a broad range of dynamical systems..

However, the relation to contact process dynamics and cusp catastrophes has only been shown for specific values of the model’s parameters Böttcher et al. (2016a). But, given the model’s parameter regime, can we generally predict the dynamics type? And does the model’s mean-field critical behavior belong to the Ising universality class or not? Here we answer these questions and analytically demonstrate that the mean-field critical behavior of the model is restricted to only three possible regimes: (a) uncorrelated spontaneous transition dynamics, (b) contact process dynamics and (c) cusp catastrophes. Cusp catastrophes can display abrupt transitions and hysteresis effects — phenomena that can harm the proper functioning of real-world networked systems since small variations in the system’s control parameters may cause catastrophic transitions from a seemingly well-functioning state to global malfunction or severe outages Achlioptas et al. (2009); Araújo and Herrmann (2010); Nagler et al. (2011, 2012); Schröder et al. (2013); Cho et al. (2013); Helbing (2013); Böttcher et al. (2015); D’Souza and Nagler (2015); Böttcher et al. (2016b).

Model

The general contagion dynamics is defined in a network whose constituents (i.e. nodes) are regarded as either active (e.g. not damaged) or inactive (e.g. failed). Three fundamental processes define the transitions between these two states Majdandzic et al. (2014); Böttcher et al. (2016a): (i) nodes undergo a spontaneous transition $A\rightarrow X$ from an active ($A$) to an inactive state ($X$) in a time interval $dt$ with probability $p\text{d}t$, (ii) if fewer than or equal to $m$ nearest neighbors of a node are active, the node becomes inactive ($Y$) due to an induced transition, i.e. $A\rightarrow Y$, with probability $r\text{d}t$ and (iii) a spontaneous reverse transition with probability $q\text{d}t$ if $X\rightarrow A$ or probability $q^{\prime}\text{d}t$ if $Y\rightarrow A$. The inactive states $X$ and $Y$ only differ in their reverse transitions and are equivalent if $q=q^{\prime}$. Process (ii) describes that a node with degree $k$ can become inactive if its number of inactive neighbors is larger or equal to $k-m$. Similar to threshold models describing complex contagion phenomena, the threshold $m$ defines the number of contacts to inactive nodes that is necessary to induce a transition as defined by process (ii) Granovetter (1978); Watts (2002); López-Pintado (2008); Gleeson (2013). A low value of $m$ corresponds to the situation where many inactive neighbors are required to sustain spreading. In contrast, for a large value of $m$ only a few inactive neighbors can sustain the spreading process.

Processes (i-iii) are illustrated in Fig. 1.

Let $a(t)\in[0,1]$ denote the total fraction of inactive nodes. Thus, $a(t)=u_{\text{spon}}(t)+u_{\text{ind}}(t)$ with $u_{\text{spon}}(t)$ and $u_{\text{ind}}(t)$ being the fractions of nodes that are inactive due to spontaneous and induced transitions respectively. The total fraction of inactive nodes in the stationary state is referred to as $a_{st}$. In accordance with Ref. Böttcher et al. (2016a) we derive the mean-field rate equations by assuming a system with homogeneous degrees in the thermodynamic limit that exhibits perfect mixing. Here, perfect mixing either refers to a network of randomly connected nodes with a sufficiently large mean degree or dynamical rewiring Buckee et al. (2004); Marro and Dickman (2005). For the fraction of nodes that spontaneously became inactive we find:

[TABLE]

where the first term accounts for the fact that active nodes spontaneously become inactive with rate $p$ (process (i)) and the second term corresponds to the spontaneous reverse transition with rate $q$ (process (iii)). Eq. (1) is exact since the network structure is not influencing these spontaneous transitions.

Induced transitions (process (ii)) can only occur for nodes whose number of active neighbors is smaller than or equal to $m$. Under the assumption of a perfectly mixed population, the probability that a node of degree $k$ is located in such a neighborhood is $E_{k}=\sum_{j=0}^{m}\binom{k}{k-j}a^{k-j}(1-a)^{j}$ Majdandzic et al. (2014); Böttcher et al. (2016a). The time evolution of the fraction of nodes that are inactive due to induced transitions is therefore given by:

[TABLE]

with $f_{k}$ being the degree distribution. The first term describes the occurrence of induced transitions (process (ii)) with rate $r$ of active nodes in a neighborhood where the number of active neighbors is smaller than or equal to $m$ whereas the second term accounts for the spontaneous reverse transition to an active state with rate $q^{\prime}$ (process (iii)). In order to study the influence of different threshold values $m$ on the mean-field critical behavior of Eqs. (1) and (2), we consider a regular network with degree $k$, i.e. the degree distribution $f_{k^{\prime}}=\delta_{kk^{\prime}}$. We will demonstrate below that the model defined by processes (i-iii) can only exhibit three different regimes depending 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^{\prime}>q$ Böttcher et al. (2016a) as a dynamical feature that does not belong to the critical behavior Strogatz (2014). The equations describing the critical behavior, i.e. $\dot{u}_{spon}=0$ and $\dot{u}_{ind}=0$, can be decoupled by multiplying one of them with an appropriate constant excluding limit cycles Strogatz (2014) — tantamount to setting $q=q^{\prime}=1$. This yields

[TABLE]

with $S(a)=\sum_{j=0}^{m}\binom{k}{k-j}a^{k-j}(1-a)^{j+1}=(1-a)E_{k}$. We use $S(a)$ as shorthand notation for the probability that an active node is located in a neighborhood that is able to induce a transition. Thus, differences in the inactive states $X$ and $Y$, i.e. different $q$ and $q^{\prime}$, do not influence the critical behavior of Eq. (3) but only rescale $r$ and $p$. In the following we analyze the stationary states of Eq. (3) that will be reached in the long-time limit.

Class (a): Uncorrelated spontaneous transitions

We start with the case $m=k$ where the number of active nodes necessary to sustain spreading has to be smaller or equal to the node’s degree $k$ according to the definition of process (ii). This describes the regime where spreading occurs independently of the neighborhood’s state such as in exogenously driven adoption dynamics Crane and Sornette (2008); Ghanbarnejad et al. (2014),

[TABLE]

since $E_{k}=1$ and $S(a)=1-a$. Eq. (4) has only one stationary state, i.e. $a_{st}(r,p)=(r+p)/(1+r+p)$, see Fig. 2 (left).

Class (b): Contact dynamics

By definition $m=k-1$ implies that $k-1$ or less neighbors of a node have to be active to induce a transition. This case describes a contact process where one inactive neighbor is sufficient to sustain spreading Marro and Dickman (2005). As demonstrated in the Appendix, we find for $p=0$ that there exists a critical $r_{c}=k^{-1}$ separating an absorbing and an active phase, i.e. $a_{st}^{(1)}(r,0)=0$ as $r\leq r_{c}$ and $a_{st}^{(2)}(r,0)>0$ as $r>r_{c}$. In the limit of $r\rightarrow r_{c}$ Eq. (3) takes the form:

[TABLE]

Equation (5) describes the mean-field contact process, SIS dynamics or Schlögl I Marro and Dickman (2005); Grassberger (1982). In the limit of $r\rightarrow r_{c}$ the order parameter scales as $a_{st}(r,0)\sim|r-r_{c}|^{\beta}$ with $\beta=1$ and adding the field-like contribution $p(1-a)$ to Eq. (5) yields $a_{st}(r_{c},p)\sim p^{1/\delta_{h}}$ as $p\rightarrow 0$ with the field exponent $\delta_{h}=2$. We illustrate in Fig. 2 (center) the occurrence of only one stable fixed point for $p>0$. This also results in a smeared-out transition for $p>0$ instead of a second-order phase transition for $p=0$ as shown in Fig. 3 (top panel). For $p=0$ and $r>r_{c}$ one clearly sees that $a_{st}^{(1)}(r,0)$ is an unstable but $a_{st}^{(2)}(r,0)$ a stable fixed point.

Class (c): Cusp catastrophes

For some values of $m$, we find a metastable region as illustrated in Fig. 3 (right). Inside this hysteresis region two stable fixed points coexist. Phase-switching is observed when fluctuations in systems of finite size push the dynamics close to the unstable fixed point, cf. Fig. 2 (right). Between the switching events the dynamics remains in one of the two phases for some time. The waiting times thus depend on the fluctuation strength and the distance from one phase to the unstable state in the phase portrait, cf. Fig. 2 (right). For $m<k-1$ where either two or more inactive neighbors are necessary to induce a transition, we now show that the corresponding metastable regions always exist due to the relation to cusp catastrophes Zeeman (1979). For a detailed analytical treatment we refer to the Appendix. The cusp point where the two bifurcation lines intersect (cf. Fig. 3 (right)) is given by $a_{0}(k,m)=(k-1-m)/(k+1)$ together with the corresponding control parameters

[TABLE]

and

[TABLE]

We illustrate the influence of different values of $m$ on $(r_{0},p_{0})$ and on the extent of the hysteresis area in the Appendix. Studying Eq. (3) in the vicinity of $(a_{0},r_{0},p_{0})$, i.e. setting $a=a_{0}+\tilde{a}$, $r=r_{0}+\tilde{r}$, $p=p_{0}+\tilde{p}$, yields for the Taylor expansion (omitted tilde)

[TABLE]

We thus find by setting $p$ or $r$ to zero respectively and solving for the fixed point of Eq. (8):

[TABLE]

[TABLE]

In previous work, the critical behavior at the cusp point of a regular random network ($k=10$ and $m=4$) has been conjectured to belong to the Ising universality class although by definition the dynamics corresponds to a general contact process Böttcher et al. (2016a).

Final remarks

We find that the critical behavior of the general contagion model as formulated in Eq. (3) does not belong to the Ising universality class but to exactly three regimes. The first regime, $m=k$, corresponds to purely spontaneous failure and recovery dynamics. For $m=k-1$ the model recovers the critical behavior of the contact process. A cusp catastrophe is found for all $m<k-1$ with the typical critical behavior at the cusp point (Eqs. (9) and (10)). This sheds analytical insight into a broad range of spreading processes that are determined by the network’s connectivity $k$ and the threshold parameter $m$, cf. examples in Tab. 1.

We have demonstrated that the phase diagram corresponds to a cusp catastrophe, when two or more inactive nodes are needed to trigger induced node-transitions. This scenario typically implies dramatic and uncontrollable global transitions in the network for many systems involving complex contagion dynamics. One could naively expect that it could be beneficial for failure control to design systems such that a component only fails if many of its neighbors already failed, i.e. delaying the failure dynamics. Our results suggest, however, that this delaying procedure might facilitate uncontrollable transitions, 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 Gross et al. (2006); Scarpino et al. (2016). For low spatial dimensions or highly structured networks, the assumptions of perfect mixing or independent node-to-node interactions are not guaranteed. Still mean-field approximations qualitatively describe a given dynamics Gleeson et al. (2012); Keeling and Rohani (2008); Marro and Dickman (2005), see given examples in the Appendix.

Future work should establish the behavior of transients as a function of threshold parameter $m$ and the topology of the network. It has been demonstrated that opinions as well as coinfections may spread faster in clustered networks compared to random ones Centola (2010); Hébert-Dufresne and Althouse (2015). This links our result to the multiple exposure condition in complex contagion phenomena.

In the study of collective behaviors, such as the adoption of innovations, the distinction between exogenous and endogenous factors is of great interest but often solely based on a contact process-like adoption model Crane and Sornette (2008); Ghanbarnejad et al. (2014). Our results suggest studying these processes within our more general framework that incorporates contact process-like adoption as one special case and can account for spreading that relies on multiple contacts.

We acknowledge financial support from the ETH Risk Center (grant RC SP 08-15) and ERC Advanced grant number FP7-319968 FlowCCS.

Appendix A Appendix

A.1 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. $\dot{u}_{spon}=0$ and $\dot{u}_{ind}=0$. Both 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^{\prime}=1$ yielding:

[TABLE]

where $S(a)=\sum_{j=0}^{m}\binom{k}{k-j}a^{k-j}\left(1-a\right)^{j+1}$.

Class (a): Uncorrelated spontaneous transitions

In the case of $m=k$, i.e. spontaneous transitions that occur independent of the neighborhood, the governing rate equation is $\dot{a}=(r+p)(1-a)-a$ since $\sum_{j=0}^{k}\binom{k}{k-j}a^{k-j}\left(1-a\right)^{j}=1$. Therefore, the corresponding stationary state is given by $a_{st}(r,p)=(r+p)/(1+r+p)$.

Class (b): Contact dynamics

For $m=k-1$ we first set $p=0$ and find

[TABLE]

since $\sum_{j=0}^{k-1}\binom{k}{k-j}a^{k-j}\left(1-a\right)^{j}=1-\left(1-a\right)^{k}$. We directly observe that there is a fixed point $a_{st}^{(1)}(r,0)=0$ at the origin of the phase portrait, cf. Fig. 2 (center). Furthermore, the function $f(a,r,p)$ of Eq. (11) is given by $f(a,r)=r\left[\left(1-a\right)-\left(1-a\right)^{k+1}\right]-a$ and we find $f^{\prime}(0,r)=rk-1$. Consequently, $f^{\prime}(0,r)>0$ if $r>k^{-1}=r_{c}$ and there exists one maximum $a_{max}(r,k)=1-\left[\frac{1}{k+1}\left(r^{-1}+1\right)\right]^{1/k}>0$ if $r>r_{c}$ ($f^{\prime\prime}(a_{max})<0$). Thus, $a_{st}^{(1)}(r,0)=0$ is the stable fixed point if $r\leq r_{c}$. For $r>r_{c}$, $a_{st}^{(1)}(r,0)=0$ is unstable with $f^{\prime}(0,r)>0$ and a maximum $a_{max}(r,k)$ exists obeying $0<a_{max}(r,k)<1$ ($f(a_{max})>0$). As a consequence of the latter fact and $f(1,r)=-1$, a second second fixed point $a_{st}^{(2)}(r,0)\in(a_{max},1)$ exists, that is stable for $r>r_{c}$ as illustrated in Fig. 2. In the vicinity of the transition point ($r\rightarrow r_{c}$), we only consider the dominant linear term in the sum of Eq. (12) and obtain:

[TABLE]

As $r\rightarrow r_{c}$ the order parameter scales as $a_{st}(r,0)\sim|r-r_{c}|^{\beta}$ with $\beta=1$ and assuming a non-zero field like term $p(1-a)$ yields $a_{st}(r_{c},p)\sim p^{1/\delta_{h}}$ as $p\rightarrow 0$ with $\delta_{h}=2$.

Class (c): Cusp catastrophes

If $m<k-1$, the critical behavior is described by a cusp catastrophe. The equilibrium point $a_{0}$ of our dynamical system $\dot{a}=f(a,r,p)$ is said to correspond to a cusp catastrophe Zeeman (1979) for the parameters $r_{0},p_{0}$ if it satisfies the following conditions Hoppensteadt and Izhikevich (2012):

1. (i)

The point $a_{0}$ is a non-hyperbolic equilibrium. 2. (ii)

The quadratic term of the function $f$ vanishes but not the cubic one, i.e.

[TABLE] 3. (iii)

The two vectors

[TABLE]

are linearly independent.

We start with condition (ii) and find that the vanishing second derivative at $a_{0}$ implies:

[TABLE]

The cases $a_{0}=0$ or $a_{0}=1$ do not correspond to an equilibrium point of Eq. 3. We therefore focus on the remaining solution $a_{0}(k,m)=(k-1-m)/(k+1)$. The third derivative at this point is given by:

[TABLE]

and clearly non-vanishing since $m<k-1$. In order to satisfy condition (i), i.e. $f(a_{0},r_{0},p_{0})=f^{\prime}(a_{0},r_{0},p_{0})=0$, we find:

[TABLE]

[TABLE]

We need to demonstrate that $r_{0}$ and $p_{0}$ are well behaved, i.e. $0<r_{0}<\infty$ and $0<p_{0}<\infty$. Therefore, we need to show (1) $S^{\prime}(a_{0})a_{0}-S(a_{0})>0$ and (2) $S(a_{0})+S^{\prime}(a_{0})(1-a_{0})>0$. The first condition leads to

[TABLE]

what is fulfilled if $\left(k-j\right)-\frac{a_{0}}{1-a_{0}}\left(j+1\right)>1$. Since $a_{0}=(k-1-m)/(k+1)$, we obtain the condition $(k+1)(m+1-j)>m+2$ which is always satisfied for $m\in\{0,\dots,k-2\}$ and $j\in\{0,\dots,m\}$. For condition (2) it is sufficient to show that $S^{\prime}(a_{0})>0$ — a consequence of the previous statement. The remaining condition (iii) for $\mathbf{v},\mathbf{w}\in\mathbb{R}^{2}$ is clearly met since

[TABLE]

are linearly independent. We have thus shown that Eq. (11) has an equilibrium point $a_{0}$ for the parameters $r_{0},~{}p_{0}$ that corresponds to a cusp catastrophe if $m<k-1$. We study Eq. (11) in the vicinity of $(a_{0},r_{0},p_{0})$ by setting $a=a_{0}+\tilde{a}$, $r=r_{0}+\tilde{r}$, $p=p_{0}+\tilde{p}$ what yields for its Taylor expansion (omitted tilde sign) Hoppensteadt and Izhikevich (2012):

[TABLE]

where $\lambda=(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:

[TABLE]

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):

[TABLE]

[TABLE]

A.2 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) Böttcher et al. (2016a). Here, we study the influence of different values of the threshold $m$ keeping $k$ fixed (Fig. S4). For $k=10$ and $m\in\{0,\dots,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:

[TABLE]

The hysteresis area decreases with $m$.

This means, the more contacts to inactive nodes necessary, the larger the metastable domain and the larger threshold values, the control parameters $r$ and $p$ for which such behavior occurs. We illustrate this effect for hysteresis areas $A_{H|r<r_{c}}$ restricted to $r<r_{c}$ in Fig. S4 (right). The black dots represent the actual values of the corresponding area and the black line is an approximation based on the assumption that $A_{H|r<r_{c}}(r_{0},p_{0})\sim(r_{c}-r_{0})p_{0}$. The value of $r_{c}$, as defined in Fig. 3 (right), is determined by Eq. (11) in the limit $p\rightarrow 0$.

A.3 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 Zeeman (1979). Therefore, our mathematical framework is valid for any choice of $f_{k}$. However, as expected, 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_{k^{\prime}}=\delta_{kk^{\prime}}$) implies multiple terms in the sum $\sum_{k^{\prime}}f_{k^{\prime}}$. A degree distribution with finite cut-off $k_{n}$ is of the form $f_{k^{\prime}}=P_{k_{1}}\delta_{k_{1}k^{\prime}}+P_{k_{2}}\delta_{k_{2}k^{\prime}}+\dots+P_{k_{n}}\delta_{k_{n}k^{\prime}}$ with normalization $\sum_{k^{\prime}\in\{k_{1},k_{2},\dots,k_{n}\}}P_{k^{\prime}}=1$. Here $P_{k_{i}}$ denotes the probability of degree $k_{i}$. We now define $S_{k_{i}}(a)=\sum_{j=0}^{m}\binom{k_{i}}{k_{i}-j}a^{k_{i}-j}(1-a)^{j+1}=(1-a)E_{k_{i}}$. Instead of Eq. 11, a general degree distribution as specified above will lead to:

[TABLE]

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_{k_{i}}\left[rS_{k_{i}}(a)+p(1-a)-a\right]$ that are similar to the ones of a regular degree 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_{i}$, different 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\geq k_{max}$, the overall dynamics is necessarily spontaneous (as $m\geq k_{i}$ for all $i$).

• •

Spontaneous critical dynamics is observed as well if, for example, the $i$-th term leads to contact dynamics ($m=k_{i}-1$) and the remaining terms lead to spontaneous ones. Since $m$ is fixed, only one $k_{i}$ is equal to $m-1$ and with $\sum_{k^{\prime}\in\{k_{1},k_{2},\dots,k_{n}\}\setminus\{k_{i}\}}P_{k^{\prime}}=1-P_{k_{i}}$ (see Sec. Critical behavior in a regular network):

[TABLE]

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:

[TABLE]

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_{k_{i}}=0$ (degrees $k_{i}$ are absent), i.e. $a_{st}(r,p=0)=r/(1+r)$ for spontaneous dynamics in a regular network.

• •

Another scenario is the occurrence of a cusp catastrophe for the $i$-th term ($m<k_{i}-1$) and spontaneous dynamics otherwise. The spontaneous terms just lead to additional constant and linear contributions in the polynomial describing the steady states:

[TABLE]

The cusp point $a_{0}(k,m)=(k-1-m)/(k+1)$ for regular degree distributions has been derived via a vanishing second order derivative, cf. Sec. Critical behavior in a regular network. Therefore, the additional constant 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:

[TABLE]

[TABLE]

We note that for $P_{k_{i}}=1$ (degrees $k_{i}$ occur P-almost surely) Eqs. (28) and (29) correspond to Eqs. (15) and (16) for a regular graph. However, unlike Eq. (16), Eq. (29) might be a negative, unphysical rate $p_{0}$. Therefore, the proportion of spontaneous terms determines if the cusp point exists for physical, i.e. positive rates $r$ and $p$ in the phase-space or not. Consequently, one might not observe a hysteresis region if spontaneous dynamics has too much influence. We illustrate a phase-space configuration for a mixture of spontaneous and cusp dynamics in Fig. 5 (right) with $P_{k_{i}}=0.95$ (degrees $k_{i}$ occur with probability $0.95$) and $m=0$.

In Fig. 6, we study the influence of a power-law degree distribution $P(k)\propto k^{-1}$ with $k_{min}=4$ and $k_{max}=100$ 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.

A.4 Stochastic simulations

In this section, we present stochastic simulations of spreading processes on real networks compared to the corresponding mean-field predictions. In Fig. S7, we considered a friendship network 111The network is based on the public-use dataset from Add Health, a program project designed by J. Richard Udry, Peter S. Bearman, and Kathleen Mullan Harris, and funded by a grant from the National Institute of Child Health and Human Development (P01-HD31921). For data files from Add Health contact Add Health, Carolina Population Center, 123 W. Franklin Street, Chapel Hill, NC 27516-2524, http://www.cpc.unc.edu/addhealth. with 2539 nodes and 20910 edges since for example epidemics or opinions spread through such networks. The corresponding degree distribution is illustrated in the inset of Fig. S7. Based on this degree distribution, we numerically computed a 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ős-Rényi network.