Physical Review. E
American Physical Society
Critical behavior of a two-step contagion model with multiple seeds
DOI 10.1103/PhysRevE.95.062115, Volume: 95, Issue: 6,
•
•
• Altmetric

### Notes

Abstract

A two-step contagion model with a single seed serves as a cornerstone for understanding the critical behaviors and underlying mechanism of discontinuous percolation transitions induced by cascade dynamics. When the contagion spreads from a single seed, a cluster of infected and recovered nodes grows without any cluster merging process. However, when the contagion starts from multiple seeds of O(N) where N is the system size, a node weakened by a seed can be infected more easily when it is in contact with another node infected by a different pathogen seed. This contagion process can be viewed as a cluster merging process in a percolation model. Here we show analytically and numerically that when the density of infectious seeds is relatively small but O(1), the epidemic transition is hybrid, exhibiting both continuous and discontinuous behavior, whereas when it is sufficiently large and reaches a critical point, the transition becomes continuous. We determine the full set of critical exponents describing the hybrid and the continuous transitions. Their critical behaviors differ from those in the single-seed case.

Choi, Lee, and Kahng: Critical behavior of a two-step contagion model with multiple seeds

## I.. INTRODUCTION

Nonequilibrium dynamic transitions driven by cascade dynamics on complex networks have attracted considerable attention recently [1–3]. The spreading of epidemic disease on complex networks [4–23] is an instance, in which a pathogen is transmitted from an infected node (e.g., a person) to a susceptible neighbor, who then becomes infected with a certain probability. If the transmission probability is sufficiently large (small), the pathogen spreads out to a macroscopic scale (remains local). An epidemic transition occurs between these two limits. The extent of spreading also depends on the structure of an underlying network [1,24]. When the degree distribution of a network is highly heterogeneous, diseases can spread out massively even for a small transmission probability, so that an epidemic transition point can be zero [25]. Information spreading in social media from one page to others may be modeled in a similar manner [5,6].

Among the several epidemic models, one of the simple contagion models is the so-called susceptible-infected-recovered (SIR) model [26,27], in which each node has one of three states, susceptible (denoted as $S$), infected ($I$), or recovered ($R$). Initially, all the nodes are in state $S$ except for one seed node in state $I$. The contagion process starts from a single node in state $I$. Each node in state $I$ transmits pathogens to its neighbors in state $S$ and infects each of them with probability $\kappa$; then it changes its state to $R$ with a unit probability. This contact process is repeated until the system reaches an absorbing state in which no infected node is left in the system. When the probability $\kappa$ is sufficiently small (large), the order parameter defined as the density of nodes in state $R$ after the system falls into the absorbing state becomes ; i.e., the system falls into a subcritical (supercritical) state. In between, an epidemic transition occurs at ${\kappa }_{c}$, and the system exhibits critical behavior. It is known that when the dynamics starts from a single seed on Erdős-Rényi (ER) random networks [28], the SIR model undergoes a continuous percolation transition following the universal behavior of ordinary percolation.

The SIR model with multiple seeds has been considered [29], in which two percolation transitions occur successively at ${\kappa }_{c1}$ and ${\kappa }_{c2}$ as $\kappa$ is increased. The density of nodes in state $R$ is finite for $\kappa >{\kappa }_{c1}$, whereas the density of nodes in state $S$ disappears for $\kappa >{\kappa }_{c2}$. Thus, there exists a state of coexisting nodes in states $R$ and $S$ between ${\kappa }_{c1}$ and ${\kappa }_{c2}$.

The SIR model was extended to a two-step contagion model, in which a weakened state ($W$) can exist between the $S$ and $I$ states. Accordingly, this model is called the SWIR model [7,20]. Nodes in state $W$ are involved in the reactions $S+I\to W+I$ and $W+I\to 2I$, which occur in addition to the reactions $S+I\to 2I$ and $I\to R$ in the SIR model. The properties of the epidemic transition in the SWIR model were extensively investigated for the single-seed case [7,10,15,20–22]. The order parameter, defined as the density of nodes in state $R$ after an absorbing state is reached, displays a discontinuous transition, whereas other physical quantities such as the outbreak size distribution exhibit critical behaviors. Thus, the phase transition occurring in the SWIR model with a single seed is regarded as a mixed-order phase transition [22]. The dynamic rule of the SWIR model is so simple that its underlying mechanism for the discontinuous behavior of the order parameter was disclosed [30]. Moreover, the mechanism turned out to be universal in other models such as $k$-core percolation [31–34], the cascading failure model on interdependent networks [35–39], and epidemic-related models [5,6,8,10,12–14,16,17].

Here we investigate the phase transitions of the SWIR model with multiple seeds. The model with multiple seeds has been investigated in Refs. [15,20,40]: The authors of Refs. [15,40] used the mean-field approach and performed numerical simulations, obtaining the phase diagram as a function of the reaction rates. The order parameter exhibits either a discontinuous or continuous transition depending on the density of the infectious seeds and mean degree of a given network [15,40]. In Ref. [20] the discontinuous transition is regarded as a spinodal transition, because there is no coexistence phase in the system while the order parameter jumps. Even though such results were obtained, the properties of the phase transitions and critical behaviors were not deeply investigated yet.

Here we reveal that the spread of contagion in the SWIR model with multiple seeds proceeds differently from that in the SWIR model with a single seed: in the multiple-seed case, the reactions $W+I\to 2I$ often occur even in early time steps, because nodes in states $W$ and $I$ involved in that reaction can originate from different seeds (see Fig. 1). We note that the number of multiple seeds was taken as $O\left(N\right)$. On the contrary, in the single-seed case, such reactions rarely occur until the system reaches a characteristic dynamic step ${n}_{c}\left(N\right)\sim {N}^{1/3}$: When dynamic step $n$ is less than ${n}_{c}\left(N\right)$, the reactions $S+I\to 2I$ and $I\to R$ are dominant, but the number of nodes in $R$ still remains as $o\left(N\right)$. The contagion spreads in the form of a branching tree. When the dynamics reaches ${n}_{c}\left(N\right)$, the branching tree forms long-range loops due to finite-size effect. Once such loops form, the reaction $W+I\to 2I$ occurs massively, in which the nodes in state $W$ were generated in early time steps. Thus, the density of nodes in state $R$ increases drastically as many as $O\left(N\right)$ in short time steps. Due to these different contagion mechanisms, the properties of epidemic transitions in the multiple seed case become different from those in the single-seed case. We will determine the full set of critical exponents describing the phase transitions in the multiple seed case, and compare them with those obtained in the single-seed case [22].

FIG. 1.
Schematic illustration of epidemic spreading processes in the SWIR model with multiple seeds. (a) Epidemic spreading begins from each infectious node. (b) These nodes can infect susceptible neighbor nodes and change their state to either $I$ or $W$. (c) A node in state $I$ contacts a node in state $W$ from a different root. (d) Then the node in state $I$ infects the node in state $W$ and changes its state to $I$. The two clusters merge.

This paper is organized as follows: In Sec. II we present the rules of the SWIR model in detail. In Sec. III we set up the self-consistency equation to derive the mean-field solution using the local tree approximation of the order parameter for the epidemic transition on the ER networks. We show that, depending on the initial density of infectious nodes, different types of phase transition can occur. In Sec. IV we report numerical results for the epidemic transitions. In the final section, a summary and discussion are presented.

## II.. THE SWIR MODEL

The SWIR model with multiple seeds is simulated on ER networks with $N$ nodes. Initially, $N{\rho }_{0}$ nodes are selected randomly from among those $N$ nodes and assigned to state $I$; the other nodes are assigned to state $S$. At each time step $n$, the following processes are performed: (1) All the nodes in state $I$ are listed in random order. (2) The states of the neighbors of each node in the list are updated sequentially as follows: If a neighbor is in state $S$, it changes its state in one of the two ways: either to $I$ with probability $\kappa$ or to $W$ with probability $\mu$. If a neighbor is in the state $W$, it changes to $I$ with probability $\eta$, where , and $\eta$ are the contagion probabilities for the respective reactions. (3) All nodes in the list change their states to $R$. This completes the time step, and we repeat the above processes until the system reaches an absorbing state in which no infectious node is left in the system. The reactions are summarized as follows:

$\mathrm{S}+\mathrm{I}\genfrac{}{}{0pt}{}{\kappa }{⟶}\mathrm{I}+\mathrm{I},$
$\mathrm{S}+\mathrm{I}\genfrac{}{}{0pt}{}{\mu }{⟶}\mathrm{W}+\mathrm{I},$
$\mathrm{W}+\mathrm{I}\genfrac{}{}{0pt}{}{\eta }{⟶}\mathrm{I}+\mathrm{I},$
$\mathrm{I}\genfrac{}{}{0pt}{}{1}{⟶}\mathrm{R}.$

The order parameter exhibits a discontinuous transition at a transition point ${\kappa }_{c}$ when ${\rho }_{0}$ is less than a critical value ${\rho }_{c}$, and it shows a continuous transition at ${\kappa }_{c}$ when ${\rho }_{0}={\rho }_{c}$ for given parameter values , and $\eta$, where $z$ is the mean degree of a given ER network.

## III.. SELF-CONSISTENCY EQUATION AND PHYSICAL SOLUTIONS

In an absorbing state, each node is in one of three states: the susceptible $S$, weakened $W$, and recovered $R$ states. The order parameter $m\left(\kappa \right)$, the density of nodes in state $R$ in an absorbing state, is written using the local tree approximation as

$m\left(\kappa \right)={\rho }_{0}+\left(1-{\rho }_{0}\right)\sum _{k=1}^{\infty }{P}_{d}\left(k\right)\sum _{\ell =1}^{k}\left(\genfrac{}{}{0pt}{}{k}{\ell }\right){q}^{\ell }{\left(1-q\right)}^{k-\ell }{P}_{R}\left(\ell \right).$
The first term in Eq. (5), ${\rho }_{0}$, is the initial density of infected nodes. In the second term, the factor $\left(1-{\rho }_{0}\right)$ represents the probability that a node is originally in state $S$. ${P}_{d}\left(k\right)$ is the probability that a randomly selected node has degree $k$; $q$ is the probability that an arbitrarily chosen edge leads to a node that is in state $R$ but not infected through the chosen edge in the absorbing state. Thus, ${P}_{d}\left(k\right)\left(\genfrac{}{}{0pt}{}{k}{\ell }\right){q}^{\ell }{\left(1-q\right)}^{k-\ell }$ is the probability that a node has degree $k$ and $\ell$ of them are in state $R$ in the absorbing state. ${P}_{R}\left(\ell \right)$ is the conditional probability that a node is finally in state $R$, provided that it was originally in state $S$ and its $\ell$ neighbors are in state $R$ in the absorbing state.

Similarly to ${P}_{R}\left(\ell \right)$, we define ${P}_{S}\left(\ell \right)$ as the conditional probability that a node remains in state $S$ in the absorbing state, provided that it has $\ell$ neighbors in state $R$ and was originally in state $S$. ${P}_{W}\left(\ell \right)$ is defined similarly. We note that for a certain node to have $\ell$ neighbors in state $R$ in the absorbing state means that the node receives $\ell$ attempts to infect it when the recovered neighbors are in state $I$. Thus, a node still remaining in state $S$ with $\ell$ neighbors in state $R$ has to be unchanged from $\ell$ infection attempts through the entire process. Thus, we obtain

${P}_{S}\left(\ell \right)={\left(1-\kappa -\mu \right)}^{\ell }.$
Next, the probability ${P}_{W}\left(\ell \right)$ is given as
${P}_{W}\left(\ell \right)=\sum _{j=0}^{\ell -1}{\left(1-\kappa -\mu \right)}^{j}\mu {\left(1-\eta \right)}^{\ell -j-1},$
where $j$ denotes the number of attacks that a node sustains before it changes to state $W$. Using the relation ${P}_{S}\left(\ell \right)+{P}_{W}\left(\ell \right)+{P}_{R}\left(\ell \right)=1$, one can determine ${P}_{R}\left(\ell \right)$ in terms of ${P}_{S}$ and ${P}_{W}$.

The local tree approximation allows us to define ${q}_{n}$ similarly to $q$ but now at the tree level $n$. The probability ${q}_{n+1}$ can be derived from ${q}_{n}$ as follows:

$\begin{array}{ccc}\hfill {q}_{n+1}& =& {\rho }_{0}+\left(1-{\rho }_{0}\right)\sum _{k=1}^{\infty }\frac{k{P}_{d}\left(k\right)}{z}\sum _{\ell =1}^{k-1}\left(\genfrac{}{}{0pt}{}{k-1}{\ell }\right)\hfill \\ & & ×\phantom{\rule{0.16em}{0ex}}{q}_{n}^{\ell }{\left(1-{q}_{n}\right)}^{k-1-\ell }{P}_{R}\left(\ell \right)\equiv {\rho }_{0}+\left(1-{\rho }_{0}\right)f\left({q}_{n}\right),\phantom{\rule{1em}{0ex}}\hfill \end{array}$
where the factor $k{P}_{d}\left(k\right)/z$ is the probability that a node connected to a randomly chosen edge has degree $k$. As a particular case, when the network is an ER network with becomes
$\begin{array}{ccc}\hfill f\left({q}_{n}\right)& =& 1-\left(1+\frac{\mu }{\eta -\kappa -\mu }\right)\hfill \\ & & ×\phantom{\rule{0.16em}{0ex}}{e}^{-\left(\kappa +\mu \right){q}_{n}z}+\frac{\mu }{\eta -\kappa -\mu }{e}^{-\eta {q}_{n}z}.\hfill \end{array}$

Equation (8) reduces to a self-consistency equation for $q$ for given epidemic parameter values in the limit $n\to \infty$. Once we obtain the solution of $q$, we can obtain the outbreak size $m\left(\kappa \right)$ using Eq. (5). For ER networks, however, $m\left(\kappa \right)$ becomes equivalent to $q$, thus the solution of the self-consistency equation (8) yields the order parameter. We remark that the methodology we used here is similar to those used in previous studies of epidemic spreading on complex networks [6,10,15,20,21].

Hereafter, we set $\mu =\kappa$ for convenience and define a function:

$F\left(m,{\rho }_{0}\right)\equiv f\left(m\right)-\frac{m}{1-{\rho }_{0}}+\frac{{\rho }_{0}}{1-{\rho }_{0}}.$
Using formula (9), we approximate $F\left(m,{\rho }_{0}\right)$ in the limit $m\to 0$ as
$F\left(m,{\rho }_{0}\right)=\frac{{\rho }_{0}}{1-{\rho }_{0}}+am+b{m}^{2}+c{m}^{3}+O\left({m}^{4}\right),$
where
$\begin{array}{ccc}\hfill a& =& \kappa z-\left[1/\left(1-{\rho }_{0}\right)\right],\hfill \end{array}$
$\begin{array}{ccc}\hfill b& =& \frac{1}{2}\kappa \left(\eta -2\kappa \right){z}^{2},\hfill \end{array}$
$\begin{array}{ccc}\hfill c& =& \frac{1}{6}\kappa \left(4{\kappa }^{2}-2\eta \kappa -{\eta }^{2}\right){z}^{3}.\hfill \end{array}$
For convenience, we neglect the higher-order terms and redefine $F\left(m,{\rho }_{0}\right)$ as
$F\left(m,{\rho }_{0}\right)\equiv \frac{{\rho }_{0}}{1-{\rho }_{0}}+am+b{m}^{2}+c{m}^{3}.$

Depending on the relative magnitudes of $a$ and $b$, various solutions of the self-consistency equation $F\left(m,{\rho }_{0}\right)=0$ can be obtained. However, we need to check whether these solutions are indeed physically relevant in the steady state when we start epidemic dynamics from a certain initial condition. The stability criterion was established in a previous work [22]: The solution $F\left({m}_{0},{\rho }_{0}\right)=0$ is stable if and only if $\partial F\left(m,{\rho }_{0}\right){/\partial m|}_{m={m}_{0}}<0$.

The equation of state in the steady state can be obtained using $F\left(m,{\rho }_{0}\right)=0$. We notice that $F\left(m=0,{\rho }_{0}\right)={\rho }_{0}/\left(1-{\rho }_{0}\right)>0$ and $F\left(m=\infty ,{\rho }_{0}\right)=-\infty$ because $c<0$, as shown in Fig. 2. We examine the solutions of $dF\left(m,{\rho }_{0}\right)/dm=0$, which are obtained as

${m}^{±}=-\frac{b}{3c}±\sqrt{\frac{{b}^{2}}{9{c}^{2}}-\frac{a}{3c}},$
where , and $c$ are given in formulas (12)–(14). Note that $a$ depends on ${\rho }_{0}$. At these extreme points has either a local maximum or a local minimum. For a given , and $\eta$, both ${m}^{±}$ values exist, and they are positive in the range of ${\kappa }_{d}<\kappa <{\kappa }_{a}$, where ${b}^{2}/9{c}^{2}-a/3c=0$ at $\kappa ={\kappa }_{d}$, and $a=0$ at $\kappa ={\kappa }_{a}$. For a given $z$ and $\eta$, diverse types of phase transitions occur depending on ${\rho }_{0}$. When ${\rho }_{0}$ is less than a certain value ${\rho }_{c}$, the order parameter jumps at a transition point. On the other hand, when ${\rho }_{0}\ge {\rho }_{c}$, the order parameter increases continuously with $\kappa$. At and $F\left({m}_{0},{\rho }_{c}\right)=0$ at $\kappa ={\kappa }_{d}={\kappa }_{c}$, as schematically shown in the blue (lower) curve in Fig. 3.
FIG. 2.
Schematic plot of $F\left(m,{\rho }_{0}\right)$ versus $m$ for $0<{\rho }_{0}<{\rho }_{c}$. Curves represent $F\left(m,{\rho }_{0}\right)$ for different $\kappa$. , and ${m}_{u}$ are the solutions of $F\left(m,{\rho }_{0}\right)=0$, where ${m}_{d}<{m}_{m}<{m}_{u}$. ${m}_{0}^{±}$ are the solutions of $F\left(m,{\rho }_{0}\right)=dF\left(m,{\rho }_{0}\right)/dm=0$ with ${m}_{0}^{-}<{m}_{0}^{+}$.
FIG. 3.
Schematic plot of $F\left(m,{\rho }_{c}\right)$ versus $m$ for ${\rho }_{0}={\rho }_{c}$. There exists a solution ${m}_{0}$ at which the self-consistency equations $F\left({m}_{0},{\rho }_{c}\right)=0$ and ${d}^{2}F\left(m,{\rho }_{c}\right)/{d}^{2}{m|}_{m={m}_{0}}=0$ hold.

### A.. When ρ0<ρc

When ${\rho }_{0}<{\rho }_{c},$ there exists a range of $\kappa$ in which $F\left(m,{\rho }_{0}\right)=0$ has more than one solution, as shown in Fig. 2. The order parameter $m$ versus $\kappa$ is shown in Fig. 4(a) and 4(b). In particular, when $\kappa$ has a certain value obtained using Eq. (16) satisfies $F\left({m}^{-},{\rho }_{0}\right)=0$. The ${m}^{-}$ value at ${\kappa }_{c}^{-}$ is denoted as ${m}_{0}^{-}$. We also define ${\kappa }_{c}^{+}$ and ${m}_{0}^{+}$ similarly to ${m}^{+}$ in Eq. (16). We note that ${\kappa }_{c}^{+}<{\kappa }_{c}^{-}$. Depending on the magnitude of the reaction probability $\kappa$ relative to ${\kappa }_{c}^{+}$ and ${\kappa }_{c}^{-}$, the order parameter behaves differently, as follows:

• (1) For $\kappa <{\kappa }_{c}^{+}$, there exists one stable solution $m={m}_{d}\left(\kappa \right)$, which increases slowly with $\kappa$. It is obtained that ${m}_{d}\approx {\rho }_{0}/\left(1-\kappa z\right)+O\left({\rho }_{0}^{2}\right)$.
FIG. 4.
(a) Schematic plot of the order parameter $m\left(\kappa \right)$ versus $\kappa$ for ${\rho }_{0}<{\rho }_{c}$. Thick solid (dashed) curves represent stable (unstable) solutions of the self-consistency equation. Dashed-dotted lines represent the pandemic probability ${P}_{\infty }\left(\kappa \right)$. (b) Plot of $m\left(\kappa \right)$ versus $\kappa$ for ER networks with mean degree $z=8$ and ${\rho }_{0}=2×{10}^{-3}$. Solid curve represents analytic solution of the self-consistency equation. Red dots (blue squares) represent averaged values of ${m}_{u}$ (${m}_{d}$) obtained by numerical simulations on ER networks of system size $N=4.096×{10}^{7}$. (c) Plot of ${P}_{\infty ,N}$ versus $\kappa$ for various system sizes $N$. Data are obtained for ER networks with mean degree $z=8$ and ${\rho }_{0}=2×{10}^{-3}$. At , which implies that ${P}_{\infty }\left(\kappa \right)$ behaves like a step function in the limit $N\to \infty$, as depicted schematically in (a) with dashed-dotted lines.
• (2) At $\kappa ={\kappa }_{c}^{+}$, there exist two solutions, ${m}_{d}$ and ${m}_{0}^{+}$ (${m}_{d}<{m}_{0}^{+}$). However, ${m}_{0}^{+}$ is not accessible because ${m}_{d}$ is stable.
• (3) When ${\kappa }_{c}^{+}<\kappa <{\kappa }_{c}^{-}$, there exist three solutions, , and ${m}_{u}$, with relative magnitudes ${m}_{d}<{m}_{m}<{m}_{u}$; however, the solution ${m}_{m}$ is unstable. Thus, only ${m}_{d}$ is accessible from the initial density ${\rho }_{0}<{m}_{d}$. The order parameter behaves as ${m}_{0}^{-}-{m}_{d}\left(\kappa \right)\sim {\left({\kappa }_{c}^{-}-\kappa \right)}^{1/2}$ for $\kappa <{\kappa }_{c}^{-}$. Thus, the critical exponent of the order parameter is obtained as $\beta =1/2$.
• (4) At $\kappa ={\kappa }_{c}^{-}$, there exist two stable solutions, ${m}_{0}^{-}$ and ${m}_{u}$. Thus, the order parameter jumps between the two values, exhibiting discontinuous behavior. Hence, a hybrid phase transition occurs at the point $\left({\kappa }_{c}^{-},{m}_{0}^{-}\right)$.
• (5) For $\kappa >{\kappa }_{c}^{-}$, there exists one solution, denoted as ${m}_{u}\left(\kappa \right)$, which increases with $\kappa$ as ${m}_{u}\left(\kappa \right)-{m}_{u}\left({\kappa }_{c}^{-}\right)\sim \left(\kappa -{\kappa }_{c}^{-}\right)$.

### B.. When ρ0=ρc

When ${\rho }_{0}={\rho }_{c}$, there exists a reaction probability ${\kappa }_{c}$ that satisfies the relation $F\left({m}_{0},{\rho }_{c}\right)=dF\left(m,{\rho }_{c}\right){/dm|}_{m={m}_{0}}=0$, and ${b}^{2}-3ac=0$. Thus, the two solutions, ${m}_{0}^{-}$ and ${m}_{0}^{+}$, reduce to the same one, which is denoted as ${m}_{0}$. The function $F\left(m,{\rho }_{c}\right)$ versus $m$ is shown in Fig. 3, and the order parameter $m$ versus $\kappa$ is shown with the analytic solution and simulation data in Fig. 5. At ${\kappa }_{c}$, singular behavior occurs, and the order parameter $m$ behaves as $m-{m}_{0}\sim {|\kappa -{\kappa }_{c}|}^{1/3}$ on both sides. The derivation of this exponent $1/3$ is presented in the Appendix.

FIG. 5.
Plot of $m\left(\kappa \right)$ versus $\kappa$ for ER networks with mean degree $z=8$ and ${\rho }_{0}={\rho }_{c}\approx 0.00747762$. Solid curve represents analytic solution of the self-consistency equation. Red dots (blue squares) represent average values of ${m}_{u}$ (${m}_{d}$) obtained by numerical simulations on ER networks of the system size $N=1.024×{10}^{7}$.

## IV.. NUMERICAL RESULTS

To estimate various critical exponents, we perform extensive numerical simulations on ER networks with mean degree $z=8$. For simplicity, the reaction probability $\mu$ is set equal to $\kappa$, and $\eta =1/2$. With these parameter values, we numerically solve Eq. (8) and determine ${\rho }_{c}\approx 0.00747762$ and ${\kappa }_{c}\approx 0.108021$. Specifically, we increase the value ${\rho }_{0}$ from zero until the values ${m}_{0}^{-}$ and ${m}_{0}^{+}$ become the same. Then ${\rho }_{c}$ is determined. Moreover, at that point, the values ${\kappa }_{c}^{-}$ and ${\kappa }_{c}^{+}$ also become the same as ${\kappa }_{c}\approx 0.108021$. The values ${\rho }_{c}$ and ${\kappa }_{c}$ will be used in numerical analysis later. For ${\rho }_{0}<{\rho }_{c}$, we take ${\rho }_{0}=2×{10}^{-3}$ in the simulations. We take the average over 50 different dynamics samples for each of 1600–4000 network configurations. Thus, 80 000–200 000 configuration averages were taken to obtain each data point.

### A.. When ρ0<ρc

Analytically we obtained that the order parameter behaves as ${m}_{0}^{-}-{m}_{d}\left(\kappa \right)\sim {\left(\mathrm{\Delta }\kappa \right)}^{\beta }$ with $\beta =1/2$ in the thermodynamic limit, where $\mathrm{\Delta }\kappa \equiv {\kappa }_{c}^{-}-\kappa$.

The dashed curve in the inset of Fig. 6 are obtained from the analytical solution Eq. (8) by taking the limit $n\to \infty$. This means that the solution is relevant in the thermodynamic limit. This dashed curve follows a solid red line in the region $\mathrm{\Delta }\kappa <\mathrm{\Delta }{\kappa }^{*}\approx {10}^{-4.2}$ (this point was estimated and indicated by an arrow), whereas it is deviated from the red line for $\mathrm{\Delta }\kappa >\mathrm{\Delta }{\kappa }^{*}$. This means that the critical behavior ${m}_{0}^{-}-{m}_{d}\sim {\left(\mathrm{\Delta }\kappa \right)}^{1/2}$ appears in the region $\mathrm{\Delta }\kappa <\mathrm{\Delta }{\kappa }^{*}$. On the other hand, the data points in the main panel were obtained by simulations for different system sizes. The dashed curve was the same drawn in the inset. A particular point on the dashed curve at $\mathrm{\Delta }{\kappa }^{*}$ is indicated by an arrow, which locates at the same $\mathrm{\Delta }{\kappa }^{*}$ in the inset. As we may see, the data points obtained from $N=1.6384×{10}^{8}$ begin to deviate from the dashed curve at this $\mathrm{\Delta }{\kappa }^{*}$. This means that the data points obtained from the system of size $N=1.6384×{10}^{8}$ for $\mathrm{\Delta }\kappa <\mathrm{\Delta }{\kappa }^{*}$ belong to the critical region, whereas other data points do to the noncritical region. This means that finite-size scaling behavior of the order parameter that appears in the form of plateaus in Fig. 6 has to be checked with the data points for the systems of sizes $N\ge \left({10}^{8}\right)$. However, it would be impractical to perform simulations with such huge system sizes.

FIG. 6.
Plot of ${m}_{0}^{-}-{m}_{d}$ versus $\mathrm{\Delta }\kappa ={\kappa }_{c}^{-}-\kappa$ in a double logarithmic scale. Data are obtained for ER networks with mean degree $z=8$ with ${\rho }_{0}=2×{10}^{-3}$. The black dashed curve both in the main panel and in the inset represents the analytical solution. For $N=1.6384×{10}^{8}$ ($▵$), crossover behavior is likely to occur at $\mathrm{\Delta }\kappa \approx {10}^{-4.2}$, which is roughly close to the point from which the analytical solution (black dashed curve in the inset) of ${m}_{0}^{-}-{m}_{d}$ deviates from the straight red line with a slope of 0.5.

Following the conventional finite-size scaling theory,

${m}_{0}^{-}\left(\infty \right)-〈{m}_{0}^{-}\left(N\right)〉\sim {N}^{-\beta /\overline{\nu }}$
at ${\kappa }_{c}^{-}$. We check this relation in Fig. 7. For small system sizes seems to be about 0.2, whereas it is estimated to be $\approx 0.24$ for large $N$. Again the crossover occurs between the system sizes $N={10}^{7}$ and ${10}^{8}$. We could obtain a more reliable value for the exponent ratio $\beta /\overline{\nu }$ for somewhat larger system sizes, but that is impractical.
FIG. 7.
Plot of ${m}_{0}^{-}-〈{m}_{0}^{-}\left(N\right)〉$ versus $N$ at $\kappa ={\kappa }_{c}^{-}$ in a double logarithmic scale. Data are obtained for ER networks with mean degree $z=8$. ${\rho }_{0}=2×{10}^{-3}$ is used. The slope of the data point corresponds to $\beta /\overline{\nu }$. As the system size is increased, crossover behavior appears in the slope from $-0.2$ to $-0.24$, which indicates that $\beta /\overline{\nu }\approx 0.24$ in the limit $N\to \infty$.

The fluctuation of the order parameter $\chi \left(\kappa \right)\equiv N\left(〈{m}_{d}^{2}〉-{〈{m}_{d}〉}^{2}\right)$ diverges as $\sim {\left({\kappa }_{c}^{-}-\kappa \right)}^{-\gamma }$. For finite systems of size $N$, it is expected that $\chi \sim {N}^{\gamma /\overline{\nu }}$ at $\kappa ={\kappa }_{c}^{-}$. From the simulation data, we obtain $\gamma /\overline{\nu }\approx 0.5$, as shown in Fig. 8.

FIG. 8.
Plot of the susceptibility $\chi$ versus $N$ at $\kappa ={\kappa }_{c}^{-}$ in a double logarithmic scale. Data are obtained for ER networks with mean degree $z=8$. ${\rho }_{0}=2×{10}^{-3}$. Here the slope of the data points corresponds to $\gamma /\overline{\nu }$, which is estimated to be $\approx 0.5±0.01$.

With the measured values $\beta /\overline{\nu }\approx 0.24$ and $\gamma /\overline{\nu }\approx 0.5$ and the analytic result $\beta =1/2$, we guess $\overline{\nu }=2$ and then $\gamma =1$. If we use those values, then the hyperscaling relation $2\beta +\gamma =\overline{\nu }$ would hold.

### B.. When ρ0=ρc

At ${\rho }_{0}={\rho }_{c}$, the jump in the order parameter does not appear, and ${m}_{0}^{+}={m}_{0}^{-}$ at $\kappa ={\kappa }_{c}$ in the thermodynamic limit. In finite systems, however, the order parameter can still exhibit a jump in some samples. Thus, the order parameter distribution $p\left(m\right)$ accumulated over different samples exhibits two separate peaks, as shown in Fig. 9. We regard the data points of $p\left(m\right)$ in the region $m<{m}_{0}$ ($m>{m}_{0}$), where ${m}_{0}$ has the theoretical value $0.171405\cdots$, as those obtained from for different samples. At $\kappa ={\kappa }_{c}$, in finite systems, we obtain the power-law behaviors ${m}_{0}-〈{m}_{0}^{-}\left(N\right)〉\sim {N}^{-\beta /\overline{\nu }}$ with $\beta /\overline{\nu }\approx 0.153$ (Fig. 10) and $〈{m}_{0}^{+}\left(N\right)〉-{m}_{0}\sim {N}^{-{\beta }^{\prime }/{\overline{\nu }}^{\prime }}$ with ${\beta }^{\prime }/{\overline{\nu }}^{\prime }\approx 0.164$ (Fig. 11).

FIG. 9.
Plot of the histogram of the order parameter $p\left(m\right)$ at ${\kappa }_{c}\approx 0.108021$. Data are obtained for ER networks of $N=2.048×{10}^{7}$ with mean degree $z=8$ and ${\rho }_{0}={\rho }_{c}\approx 0.00747762$. Even though simulations were performed at ${\kappa }_{c}$ and ${\rho }_{c}$, the distribution of the order parameter exhibits two peaks, a prototypical pattern of a discontinuous transition due to the finite size effect. As $N$ is increased, we expect that the two peaks converge and become a single peak.
FIG. 10.
Plot of ${m}_{0}-〈{m}_{0}^{-}\left(N\right)〉$ versus $N$ at ${\kappa }_{c}\approx 0.108021$. Data are obtained for ER networks with mean degree $z=8$. ${\rho }_{0}={\rho }_{c}$ is taken as $\approx 0.00747762$. $\beta /\overline{\nu }$ is estimated to be $\approx 0.153±0.003$.
FIG. 11.
Plot of $〈{m}_{0}^{+}\left(N\right)〉-{m}_{0}$ versus $N$ at ${\kappa }_{c}\approx 0.108021$. Data are obtained for ER networks with mean degree $z=8$. ${\rho }_{0}={\rho }_{c}$ is taken as $\approx 0.00747762$. ${\beta }^{\prime }/{\overline{\nu }}^{\prime }$ is estimated to be $\approx 0.164±0.005$.

For $\kappa <{\kappa }_{c}$, the fluctuation of the order parameter $\chi \equiv N\left(〈{m}_{d}^{2}〉-{〈{m}_{d}〉}^{2}\right)$ behaves as $\sim {N}^{\gamma /\overline{\nu }}G\left[\left({\kappa }_{c}-\kappa \right){N}^{1/\overline{\nu }}\right]$ with a certain scaling function $G$. On the other hand, for $\kappa >{\kappa }_{c}$, we obtain that ${\chi }^{\prime }\equiv N\left(〈{m}_{u}^{2}〉-{〈{m}_{u}〉}^{2}\right)$ behaves as $\sim {N}^{{\gamma }^{\prime }/{\overline{\nu }}^{\prime }}{G}^{\prime }\left[\left(\kappa -{\kappa }_{c}\right){N}^{1/{\overline{\nu }}^{\prime }}\right]$ with a certain scaling function ${G}^{\prime }$. We numerically obtain $\gamma /\overline{\nu }\approx 0.69$ (Fig. 12) and ${\gamma }^{\prime }/{\overline{\nu }}^{\prime }\approx 0.6$ (Fig. 13).

FIG. 12.
Plot of the susceptibility $\chi$, the fluctuation of the order parameter ${m}_{d}$, as a function of the system size $N$ at ${\kappa }_{c}\approx 0.108021$. $\gamma /\overline{\nu }$ is estimated to be $\approx 0.69±0.005$. Data are obtained for ER networks with $z=8$. ${\rho }_{0}$ is taken as ${\rho }_{c}\approx 0.00747762$.
FIG. 13.
Plot of the susceptibility ${\chi }^{\prime }$, the fluctuation of the order parameter ${m}_{u}$, as a function of the system size $N$ at ${\kappa }_{c}\approx 0.108021$. ${\gamma }^{\prime }/{\overline{\nu }}^{\prime }$ is estimated to be $\approx 0.6±0.005$. Data are obtained for ER networks with mean degree $z=8$. ${\rho }_{0}={\rho }_{c}$ is taken as $\approx 0.00747762$.

On the basis of the numerically obtained values $\beta /\overline{\nu }\approx 0.53$ and $\gamma /\overline{\nu }\approx 0.69$, and the theoretical value $\beta =1/3$, we estimate $\overline{\nu }\approx 2.179$ and $\gamma \approx 1.5$. Those values are confirmed in Fig. 14 for $\chi {N}^{-\gamma /\overline{\nu }}$ versus $\left({\kappa }_{c}-\kappa \right){N}^{1/\overline{\nu }}$, in which the data collapse well with the choices of $\overline{\nu }\approx 2.13±0.1$ and $\gamma /\overline{\nu }\approx 0.69$. The measured values of the exponents satisfy the hyperscaling relation $\left(2\beta +\gamma \right)/\overline{\nu }\approx 0.996$ well. Similarly, for $\kappa >{\kappa }_{c}$, on the basis of the numerical values ${\beta }^{\prime }/{\overline{\nu }}^{\prime }\approx 0.164$ and ${\gamma }^{\prime }/{\overline{\nu }}^{\prime }\approx 0.6$, and the theoretical value ${\beta }^{\prime }=1/3$, we obtain ${\overline{\nu }}^{\prime }\approx 2.03$ and ${\gamma }^{\prime }\approx 1.218$. Data for ${\chi }^{\prime }$ for different system sizes collapse well into a single curve with the choices of ${\overline{\nu }}^{\prime }=2.13±0.1$ and ${\gamma }^{\prime }/{\overline{\nu }}^{\prime }=0.6$ (Fig. 15). These values yield $\left(2{\beta }^{\prime }+{\gamma }^{\prime }\right)/{\overline{\nu }}^{\prime }\approx 0.91-0.93$, which deviates slightly from the expected value of unity that would satisfy the hyperscaling relation. To obtain those results, we used the numerical values ${\rho }_{c}\approx 0.00747762$ and ${\kappa }_{c}\approx 0.108021$. We remark that $\beta ={\beta }^{\prime }=1/3$ is obtained analytically.

FIG. 14.
Scaling plot of the susceptibility $\chi$ for $\kappa <{\kappa }_{c}$ in the form $\chi {N}^{-\gamma /\overline{\nu }}$ versus $\left(\kappa -{\kappa }_{c}\right){N}^{1/\overline{\nu }}$, where $\overline{\nu }\approx 2.13$ and $\gamma \approx 1.47$ are used. Data are obtained for ER networks with mean degree $z=8$. ${\rho }_{0}={\rho }_{c}$ is taken as $\approx 0.00747762$.
FIG. 15.
Scaling plot of the susceptibility ${\chi }^{\prime }$ for $\kappa >{\kappa }_{c}$ in the form ${\chi }^{\prime }{N}^{-{\gamma }^{\prime }/{\overline{\nu }}^{\prime }}$ versus $\left(\kappa -{\kappa }_{c}\right){N}^{1/{\overline{\nu }}^{\prime }}$, where ${\overline{\nu }}^{\prime }\approx 2.13$ and ${\gamma }^{\prime }\approx 1.28$ are used. Data are obtained for ER networks with mean degree $z=8$. ${\rho }_{0}={\rho }_{c}$ is taken as $\approx 0.00747762$.

## V.. SUMMARY AND DISCUSSION

We investigated the properties of phase transitions in the SWIR model with a finite density ${\rho }_{0}$ of initially infected seeds [7]. A node in the state $S$ can change its state to weakened ($W$) or infected ($I$) when it comes in contact with an infected node from the same or a different root. A weakened node can also change its state to infected ($I$) when it contacts an infected node from the same or a different root. The reaction probabilities $\kappa$ and $\mu$ in Eqs. (1) and (2), respectively, serve as control parameters. For convenience, we take $\kappa =\mu$. We found that for a given network, there exists a critical density of seeds ${\rho }_{c}$ such that for ${\rho }_{0}<{\rho }_{c}$, the order parameter, the density of nodes in state $R$ in the absorbing state, increases continuously with the critical exponent $\beta =1/2$ as $\kappa$ is increased up to a transition point ${\kappa }_{c}^{-}$ and then jumps to a finite value, followed by a continuous increase. Accordingly, the order parameter behaves as $m\left(\kappa \right)={m}_{0}^{-}-b{\left({\kappa }_{c}^{-}-\kappa \right)}^{1/2}$ for $\kappa <{\kappa }_{c}^{-}$, where $b$ is a positive constant. At ${\kappa }_{c}^{-}$, the order parameter is discontinuous by $\mathrm{\Delta }m={m}_{u}\left({\kappa }_{c}^{-}\right)-{m}_{0}^{-}$. Thus, the order parameter itself exhibits a hybrid phase transition. This pattern is different from that for the single-seed case, in which the order parameter jumps from $m=0$ to a finite value, and thus $\beta =0$. The fluctuation of the order parameter diverges at the transition point ${\kappa }_{c}^{-}$ according to a power law with the exponent $\gamma$. For the correlation size exponent $\overline{\nu }$ measured in finite systems, we find that the hyperscaling relation $2\beta +\gamma =\overline{\nu }$ holds reasonably well.

As ${\rho }_{0}$ is increased, the jump shrinks and becomes zero at ${\rho }_{c}$. For ${\rho }_{0}={\rho }_{c}$, the transition becomes continuous. We determined a complete set of critical exponents describing the phase transition at ${\kappa }_{c}$. The critical exponents are listed in Table I.

TABLE I.
List of the critical exponents of the SWIR models with a single seed and with multiple seeds.
Type
Single seed03
Multiple seeds1/212
Multiple seeds1/31/3

## ACKNOWLEDGMENT

This work was supported by the National Research Foundation of Korea by Grant No. NRF-2014R1A3A2069005.

## Appendices

### DERIVATION OF THE CRITICAL EXPONENT $\beta$ AT ${\rho }_{c}$

Here we introduce an analytical method to determine the critical exponents $\beta =1/3$ at ${\rho }_{0}={\rho }_{c}$. It was already noted in Sec. III B that for ${\rho }_{0}={\rho }_{c}$,

$F\left({\kappa }_{c},{m}_{0}\right)=\frac{dF}{dm}{|}_{{\kappa }_{c},{m}_{0}}=\frac{{d}^{2}F}{d{m}^{2}}{|}_{{\kappa }_{c},{m}_{0}}=0.$
We consider a line of the solution $F\left(\kappa ,m\right)=0$ near $\left({\kappa }_{c},{m}_{0}\right)$ by expanding $F\left({\kappa }_{c}+\delta \kappa ,{m}_{0}+\delta m\right)$ as
$\begin{array}{ccc}& & F\left({\kappa }_{c}+\delta \kappa ,{m}_{0}+\delta m\right)\hfill \\ & & \phantom{\rule{1em}{0ex}}\simeq \frac{\partial F}{\partial \kappa }{|}_{{\kappa }_{c},{m}_{0}}\left(\delta \kappa \right)+\frac{1}{6}\frac{{\partial }^{3}F}{\partial {m}^{3}}{|}_{{\kappa }_{c},{m}_{0}}{\left(\delta m\right)}^{3}+\cdots =0\phantom{\rule{2.em}{0ex}}\hfill \end{array}$
where only nonzero terms are considered. Since $\delta \kappa$ and ${\left(\delta m\right)}^{3}$ are two lowest terms in Eq. (A2) and their coefficients have the opposite sign to each other, $\delta m\sim {\left(\delta \kappa \right)}^{1/3}$ when $\delta \kappa \ll 1$. Thus for both cases of $\kappa <\left(>\right){\kappa }_{c}$, the critical exponents $\beta ={\beta }^{\prime }=1/3$.

## REFERENCES

[1]

S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Rev. Mod. Phys.80, 1275 (2008).RMPHAT0034-6861

[2]

N. Araújo, P. Grassberger, B. Kahng, K. J. Schrenk, and R. M. Ziff, Eur. Phys. J. Special Topics223, 2307 (2014).1951-6355, doi: 10.1140/epjst/e2014-02266-y

[3]

D. Lee, Y. S. Cho, and B. Kahng, J. Stat. Mech. (2016) P1240021742-5468, doi: 10.1088/1742-5468/2016/12/124002

[4]

R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Rev. Mod. Phys.87, 925 (2015).RMPHAT0034-6861

[5]

D. J. Watts, Proc. Natl. Acac. Sci. USA99, 5766 (2002).PNASA60027-8424

[6]

P. S. Dodds and D. J. Watts, Phys. Rev. Lett.92, 218701 (2004).PRLTAO0031-9007

[7]

H.-K. Janssen, M. Müller, and O. Stenull, Phys. Rev. E70, 026114 (2004).PLEEE81539-3755

[8]

P. L. Krapivsky, S. Redner, and D. Volovik, J. Stat. Mech. (2011) P120031742-5468, doi: 10.1088/1742-5468/2011/12/P12003

[9]

S. Ni, W. Weng, and H. Zhang, Physica A390, 4528 (2011).PHYADX0378-4371

[10]

G. Bizhani, M. Paczuski, and P. Grassberger, Phys. Rev. E86, 011128 (2012).PLEEE81539-3755

[11]

N. Crokidakis and S. M. D. Queirós, J. Stat. Mech. (2012) P060031742-5468, doi: 10.1088/1742-5468/2012/06/P06003

[12]

L. Hébert-Dufresne, O. Patterson-Lomba, G. M. Goerg, and B. M. Althouse, Phys. Rev. Lett.110, 108103 (2013).PRLTAO0031-9007

[13]

L. Chen, F. Ghanbarnejad, W. Chai, and P. Grassberger, Europhys. Lett.104, 50001 (2013).EULEEJ0295-5075

[14]

S. Melnik, J. A. Ward, J. P. Gleeson, and M. A. Porter, Chaos23, 013124 (2013).CHAOEH1054-1500

[15]

T. Hasegawa and K. Nemoto, J. Stat. Mech. (2014) P110241742-5468, doi: 10.1088/1742-5468/2014/11/P11024

[16]

W. Cai, L. Chen, F. Ghanbarnejad, and P. Grassberger, Nat. Phys.11, 936 (2015).1745-2473, doi: 10.1038/nphys3457

[17]

L. Hébert-Dufresne and B. M. Althouse, Proc. Natl. Acad. Sci. USA112, 10551 (2015).PNASA60027-8424

[18]

Q. Zhang and D. Wang, Intl. J. Environ. Res. Public Health12, 9750 (2015).1660-4601, doi: 10.3390/ijerph120809750

[19]

A. S. Mata and S. C. Ferreira, Phys. Rev. E91, 012816 (2015).PLEEE81539-3755

[20]

H.-K. Janssen and O. Stenull, Europhys. Lett.113, 26005 (2016).EULEEJ0295-5075

[21]

K. Chung, Y. Baek, M. Ha, and H. Jeong, Phys. Rev. E93, 052304 (2016).1539-3755, doi: 10.1103/PhysRevE.93.052304

[22]

W. Choi, D. Lee, and B. Kahng, Phys. Rev. E95, 022304 (2017).1539-3755, doi: 10.1103/PhysRevE.95.022304

[23]

M. A. Pires and N. Crokidakis, Physica A467, 167 (2017).PHYADX0378-4371

[24]

R. Albert and A.-L. Barabási, Rev. Mod. Phys.74, 47 (2002).RMPHAT0034-6861

[25]

R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett.86, 3200 (2001).PRLTAO0031-9007

[26]

D. Mollison, J. R. Stat. Soc. B39, 283 (1977).

[27]

M. E. J. Newman, Phys. Rev. E66, 016128 (2002).1063-651X, doi: 10.1103/PhysRevE.66.016128

[28]

P. Erdős and A. Rényi, Publ. Math. Debrecen6, 290 (1959).

[29]

T. Hasegawa and K. Nemoto, Phys. Rev. E93, 032324 (2016).1539-3755, doi: 10.1103/PhysRevE.93.032324

[30]

D. Lee, W. Choi, J. Kertéz, and B. Kahng, arXiv:1608.00776 (2016).

[31]

J. Chalupa, P. L. Leath, and G. R. Reich, J. Phys. C12, L31 (1979).JPSOAW0022-3719

[32]

S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Phys. Rev. Lett.96, 040601 (2006).PRLTAO0031-9007

[33]

G. J. Baxter, S. N. Dorogovtsev, K. E. Lee, J. F. F. Mendes, and A. V. Goltsev, Phys. Rev. X5, 031017 (2015).2160-3308, doi: 10.1103/PhysRevX.5.031017

[34]

D. Lee, M. Jo, and B. Kahng, Phys. Rev. E94, 062307 (2016).1539-3755, doi: 10.1103/PhysRevE.94.062307

[35]

S. V. Buldyrev, R. Parshani, G. Paul, H. E. Stanley, and S. Havlin, Nature (London)464, 1025 (2010).NATUAS0028-0836

[36]

S.-W. Son, P. Grassberger, and M. Paczuski, Phys. Rev. Lett.107, 195702 (2011).PRLTAO0031-9007

[37]

D. Zhou, A. Bashan, R. Cohen, Y. Berezin, N. Shnerb, and S. Havlin, Phys. Rev. E90, 012803 (2014).PLEEE81539-3755

[38]

S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Genio, J. Gómez-Gardeñes, M. Romance, I. Sendina-Nadal, Z. Wang, and M. Zanin, Phys. Rep.544, 1 (2014).PRPLCM0370-1573

[39]

D. Lee, S. Choi, M. Stippinger, J. Kertesz, and B. Kahng, Phys. Rev. E93, 042109 (2016).1539-3755, doi: 10.1103/PhysRevE.93.042109

[40]

T. Hasegawa and K. Nemoto, arXiv:1611.02809 (2016).

Citing articles via
https://www.researchpad.co/tools/openurl?pubtype=article&doi=10.1103/PhysRevE.95.062115&title=Critical behavior of a two-step contagion model with multiple seeds&author=Wonjun Choi,Deokjae Lee,B. Kahng,&keyword=&subject=Articles,Statistical Physics,