Physical Review. E
American Physical Society
image
Critical behavior of a two-step contagion model with multiple seeds
Volume: 95, Issue: 6
DOI 10.1103/PhysRevE.95.062115
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 κ; 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 κ 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 o(N) [O(N)]; i.e., the system falls into a subcritical (supercritical) state. In between, an epidemic transition occurs at κ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 κc1 and κc2 as κ is increased. The density of nodes in state R is finite for κ>κc1, whereas the density of nodes in state S disappears for κ>κc2. Thus, there exists a state of coexisting nodes in states R and S between κc1 and κ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+IW+I and W+I2I, which occur in addition to the reactions S+I2I and IR 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+I2I 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(N). On the contrary, in the single-seed case, such reactions rarely occur until the system reaches a characteristic dynamic step nc(N)N1/3: When dynamic step n is less than nc(N), the reactions S+I2I and IR are dominant, but the number of nodes in R still remains as o(N). The contagion spreads in the form of a branching tree. When the dynamics reaches nc(N), the branching tree forms long-range loops due to finite-size effect. Once such loops form, the reaction W+I2I 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(N) 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].

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.
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ρ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 κ or to W with probability μ. If a neighbor is in the state W, it changes to I with probability η, where κ, μ, and η 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:

S+IκI+I,
S+IμW+I,
W+IηI+I,
I1R.

The order parameter exhibits a discontinuous transition at a transition point κc when ρ0 is less than a critical value ρc, and it shows a continuous transition at κc when ρ0=ρc for given parameter values z, μ, and η, 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(κ), the density of nodes in state R in an absorbing state, is written using the local tree approximation as

m(κ)=ρ0+(1ρ0)k=1Pd(k)=1kkq(1q)kPR().
The first term in Eq. (5), ρ0, is the initial density of infected nodes. In the second term, the factor (1ρ0) represents the probability that a node is originally in state S. Pd(k) 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, Pd(k)kq(1q)k is the probability that a node has degree k and of them are in state R in the absorbing state. PR() is the conditional probability that a node is finally in state R, provided that it was originally in state S and its neighbors are in state R in the absorbing state.

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

PS()=(1κμ).
Next, the probability PW() is given as
PW()=j=01(1κμ)jμ(1η)j1,
where j denotes the number of attacks that a node sustains before it changes to state W. Using the relation PS()+PW()+PR()=1, one can determine PR() in terms of PS and PW.

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

qn+1=ρ0+(1ρ0)k=1kPd(k)z=1k1k1×qn(1qn)k1PR()ρ0+(1ρ0)f(qn),
where the factor kPd(k)/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 Pd(k)=zkez/k!, f(qn) becomes
f(qn)=11+μηκμ×e(κ+μ)qnz+μηκμeηqnz.

Equation (8) reduces to a self-consistency equation for q for given epidemic parameter values in the limit n. Once we obtain the solution of q, we can obtain the outbreak size m(κ) using Eq. (5). For ER networks, however, m(κ) 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 μ=κ for convenience and define a function:

F(m,ρ0)f(m)m1ρ0+ρ01ρ0.
Using formula (9), we approximate F(m,ρ0) in the limit m0 as
F(m,ρ0)=ρ01ρ0+am+bm2+cm3+O(m4),
where
a=κz[1/(1ρ0)],
b=12κ(η2κ)z2,
c=16κ(4κ22ηκη2)z3.
For convenience, we neglect the higher-order terms and redefine F(m,ρ0) as
F(m,ρ0)ρ01ρ0+am+bm2+cm3.

Depending on the relative magnitudes of a and b, various solutions of the self-consistency equation F(m,ρ0)=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(m0,ρ0)=0 is stable if and only if F(m,ρ0)/m|m=m0<0.

The equation of state in the steady state can be obtained using F(m,ρ0)=0. We notice that F(m=0,ρ0)=ρ0/(1ρ0)>0 and F(m=,ρ0)= because c<0, as shown in Fig. 2. We examine the solutions of dF(m,ρ0)/dm=0, which are obtained as

m±=b3c±b29c2a3c,
where a, b, and c are given in formulas (12)–(14). Note that a depends on ρ0. At these extreme points m±, F(m,ρ0) has either a local maximum or a local minimum. For a given ρ0, z, and η, both m± values exist, and they are positive in the range of κd<κ<κa, where b2/9c2a/3c=0 at κ=κd, and a=0 at κ=κa. For a given z and η, diverse types of phase transitions occur depending on ρ0. When ρ0 is less than a certain value ρc, the order parameter jumps at a transition point. On the other hand, when ρ0ρc, the order parameter increases continuously with κ. At ρ0=ρc, m+=m=m0 and F(m0,ρc)=0 at κ=κd=κc, as schematically shown in the blue (lower) curve in Fig. 3.
Schematic plot of F(m,ρ0) versus m for 0<ρ0<ρc. Curves represent F(m,ρ0) for different κ. md, mm, and mu are the solutions of F(m,ρ0)=0, where md<mm<mu. m0± are the solutions of F(m,ρ0)=dF(m,ρ0)/dm=0 with m0−<m0+.
FIG. 2.
Schematic plot of F(m,ρ0) versus m for 0<ρ0<ρc. Curves represent F(m,ρ0) for different κ. md, mm, and mu are the solutions of F(m,ρ0)=0, where md<mm<mu. m0± are the solutions of F(m,ρ0)=dF(m,ρ0)/dm=0 with m0<m0+.
Schematic plot of F(m,ρc) versus m for ρ0=ρc. There exists a solution m0 at which the self-consistency equations F(m0,ρc)=0 and d2F(m,ρc)/d2m|m=m0=0 hold.
FIG. 3.
Schematic plot of F(m,ρc) versus m for ρ0=ρc. There exists a solution m0 at which the self-consistency equations F(m0,ρc)=0 and d2F(m,ρc)/d2m|m=m0=0 hold.

A.. When ρ0<ρc

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

    • (1) For κ<κc+, there exists one stable solution m=md(κ), which increases slowly with κ. It is obtained that mdρ0/(1κz)+O(ρ02).
      (a) Schematic plot of the order parameter m(κ) versus κ for ρ0<ρc. Thick solid (dashed) curves represent stable (unstable) solutions of the self-consistency equation. Dashed-dotted lines represent the pandemic probability P∞(κ). (b) Plot of m(κ) versus κ for ER networks with mean degree z=8 and ρ0=2×10−3. Solid curve represents analytic solution of the self-consistency equation. Red dots (blue squares) represent averaged values of mu (md) obtained by numerical simulations on ER networks of system size N=4.096×107. (c) Plot of P∞,N versus κ for various system sizes N. Data are obtained for ER networks with mean degree z=8 and ρ0=2×10−3. At κ=κc−≈0.11495, dP∞,N/dκ∼N1/2, which implies that P∞(κ) behaves like a step function in the limit N→∞, as depicted schematically in (a) with dashed-dotted lines.
      FIG. 4.
      (a) Schematic plot of the order parameter m(κ) versus κ for ρ0<ρc. Thick solid (dashed) curves represent stable (unstable) solutions of the self-consistency equation. Dashed-dotted lines represent the pandemic probability P(κ). (b) Plot of m(κ) versus κ for ER networks with mean degree z=8 and ρ0=2×103. Solid curve represents analytic solution of the self-consistency equation. Red dots (blue squares) represent averaged values of mu (md) obtained by numerical simulations on ER networks of system size N=4.096×107. (c) Plot of P,N versus κ for various system sizes N. Data are obtained for ER networks with mean degree z=8 and ρ0=2×103. At κ=κc0.11495, dP,N/dκN1/2, which implies that P(κ) behaves like a step function in the limit N, as depicted schematically in (a) with dashed-dotted lines.
    • (2) At κ=κc+, there exist two solutions, md and m0+ (md<m0+). However, m0+ is not accessible because md is stable.
    • (3) When κc+<κ<κc, there exist three solutions, md, mm, and mu, with relative magnitudes md<mm<mu; however, the solution mm is unstable. Thus, only md is accessible from the initial density ρ0<md. The order parameter behaves as m0md(κ)(κcκ)1/2 for κ<κc. Thus, the critical exponent of the order parameter is obtained as β=1/2.
    • (4) At κ=κc, there exist two stable solutions, m0 and mu. Thus, the order parameter jumps between the two values, exhibiting discontinuous behavior. Hence, a hybrid phase transition occurs at the point (κc,m0).
    • (5) For κ>κc, there exists one solution, denoted as mu(κ), which increases with κ as mu(κ)mu(κc)(κκc).

B.. When ρ0=ρc

When ρ0=ρc, there exists a reaction probability κc that satisfies the relation F(m0,ρc)=dF(m,ρc)/dm|m=m0=0, and b23ac=0. Thus, the two solutions, m0 and m0+, reduce to the same one, which is denoted as m0. The function F(m,ρc) versus m is shown in Fig. 3, and the order parameter m versus κ is shown with the analytic solution and simulation data in Fig. 5. At κc, singular behavior occurs, and the order parameter m behaves as mm0|κκc|1/3 on both sides. The derivation of this exponent 1/3 is presented in the Appendix.

Plot of m(κ) versus κ for ER networks with mean degree z=8 and ρ0=ρc≈0.00747762. Solid curve represents analytic solution of the self-consistency equation. Red dots (blue squares) represent average values of mu (md) obtained by numerical simulations on ER networks of the system size N=1.024×107.
FIG. 5.
Plot of m(κ) versus κ for ER networks with mean degree z=8 and ρ0=ρc0.00747762. Solid curve represents analytic solution of the self-consistency equation. Red dots (blue squares) represent average values of mu (md) obtained by numerical simulations on ER networks of the system size N=1.024×107.

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 μ is set equal to κ, and η=1/2. With these parameter values, we numerically solve Eq. (8) and determine ρc0.00747762 and κc0.108021. Specifically, we increase the value ρ0 from zero until the values m0 and m0+ become the same. Then ρc is determined. Moreover, at that point, the values κc and κc+ also become the same as κc0.108021. The values ρc and κc will be used in numerical analysis later. For ρ0<ρc, we take ρ0=2×103 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 m0md(κ)(Δκ)β with β=1/2 in the thermodynamic limit, where Δκκcκ.

The dashed curve in the inset of Fig. 6 are obtained from the analytical solution Eq. (8) by taking the limit n. This means that the solution is relevant in the thermodynamic limit. This dashed curve follows a solid red line in the region Δκ<Δκ*104.2 (this point was estimated and indicated by an arrow), whereas it is deviated from the red line for Δκ>Δκ*. This means that the critical behavior m0md(Δκ)1/2 appears in the region Δκ<Δκ*. 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 Δκ* is indicated by an arrow, which locates at the same Δκ* in the inset. As we may see, the data points obtained from N=1.6384×108 begin to deviate from the dashed curve at this Δκ*. This means that the data points obtained from the system of size N=1.6384×108 for Δκ<Δκ* 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(108). However, it would be impractical to perform simulations with such huge system sizes.

Plot of m0−−md versus Δκ=κc−−κ in a double logarithmic scale. Data are obtained for ER networks with mean degree z=8 with ρ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×108 (▵), crossover behavior is likely to occur at Δκ≈10−4.2, which is roughly close to the point from which the analytical solution (black dashed curve in the inset) of m0−−md deviates from the straight red line with a slope of 0.5.
FIG. 6.
Plot of m0md versus Δκ=κcκ in a double logarithmic scale. Data are obtained for ER networks with mean degree z=8 with ρ0=2×103. The black dashed curve both in the main panel and in the inset represents the analytical solution. For N=1.6384×108 (), crossover behavior is likely to occur at Δκ104.2, which is roughly close to the point from which the analytical solution (black dashed curve in the inset) of m0md deviates from the straight red line with a slope of 0.5.

Following the conventional finite-size scaling theory,

m0()m0(N)Nβ/ν¯
at κc. We check this relation in Fig. 7. For small system sizes N, β/ν¯ seems to be about 0.2, whereas it is estimated to be 0.24 for large N. Again the crossover occurs between the system sizes N=107 and 108. We could obtain a more reliable value for the exponent ratio β/ν¯ for somewhat larger system sizes, but that is impractical.
Plot of m0−−〈m0−(N)〉 versus N at κ=κc− in a double logarithmic scale. Data are obtained for ER networks with mean degree z=8. ρ0=2×10−3 is used. The slope of the data point corresponds to β/ν¯. As the system size is increased, crossover behavior appears in the slope from −0.2 to −0.24, which indicates that β/ν¯≈0.24 in the limit N→∞.
FIG. 7.
Plot of m0m0(N) versus N at κ=κc in a double logarithmic scale. Data are obtained for ER networks with mean degree z=8. ρ0=2×103 is used. The slope of the data point corresponds to β/ν¯. As the system size is increased, crossover behavior appears in the slope from 0.2 to 0.24, which indicates that β/ν¯0.24 in the limit N.

The fluctuation of the order parameter χ(κ)N(md2md2) diverges as (κcκ)γ. For finite systems of size N, it is expected that χNγ/ν¯ at κ=κc. From the simulation data, we obtain γ/ν¯0.5, as shown in Fig. 8.

Plot of the susceptibility χ versus N at κ=κc− in a double logarithmic scale. Data are obtained for ER networks with mean degree z=8. ρ0=2×10−3. Here the slope of the data points corresponds to γ/ν¯, which is estimated to be ≈0.5±0.01.
FIG. 8.
Plot of the susceptibility χ versus N at κ=κc in a double logarithmic scale. Data are obtained for ER networks with mean degree z=8. ρ0=2×103. Here the slope of the data points corresponds to γ/ν¯, which is estimated to be 0.5±0.01.

With the measured values β/ν¯0.24 and γ/ν¯0.5 and the analytic result β=1/2, we guess ν¯=2 and then γ=1. If we use those values, then the hyperscaling relation 2β+γ=ν¯ would hold.

B.. When ρ0=ρc

At ρ0=ρc, the jump in the order parameter does not appear, and m0+=m0 at κ=κ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(m) accumulated over different samples exhibits two separate peaks, as shown in Fig. 9. We regard the data points of p(m) in the region m<m0 (m>m0), where m0 has the theoretical value 0.171405, as those obtained from m0(N) [m0+(N)] for different samples. At κ=κc, in finite systems, we obtain the power-law behaviors m0m0(N)Nβ/ν¯ with β/ν¯0.153 (Fig. 10) and m0+(N)m0Nβ/ν¯ with β/ν¯0.164 (Fig. 11).

Plot of the histogram of the order parameter p(m) at κc≈0.108021. Data are obtained for ER networks of N=2.048×107 with mean degree z=8 and ρ0=ρc≈0.00747762. Even though simulations were performed at κc and ρ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. 9.
Plot of the histogram of the order parameter p(m) at κc0.108021. Data are obtained for ER networks of N=2.048×107 with mean degree z=8 and ρ0=ρc0.00747762. Even though simulations were performed at κc and ρ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.
Plot of m0−〈m0−(N)〉 versus N at κc≈0.108021. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as ≈0.00747762. β/ν¯ is estimated to be ≈0.153±0.003.
FIG. 10.
Plot of m0m0(N) versus N at κc0.108021. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as 0.00747762. β/ν¯ is estimated to be 0.153±0.003.
Plot of 〈m0+(N)〉−m0 versus N at κc≈0.108021. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as ≈0.00747762. β′/ν¯′ is estimated to be ≈0.164±0.005.
FIG. 11.
Plot of m0+(N)m0 versus N at κc0.108021. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as 0.00747762. β/ν¯ is estimated to be 0.164±0.005.

For κ<κc, the fluctuation of the order parameter χN(md2md2) behaves as Nγ/ν¯G[(κcκ)N1/ν¯] with a certain scaling function G. On the other hand, for κ>κc, we obtain that χN(mu2mu2) behaves as Nγ/ν¯G[(κκc)N1/ν¯] with a certain scaling function G. We numerically obtain γ/ν¯0.69 (Fig. 12) and γ/ν¯0.6 (Fig. 13).

Plot of the susceptibility χ, the fluctuation of the order parameter md, as a function of the system size N at κc≈0.108021. γ/ν¯ is estimated to be ≈0.69±0.005. Data are obtained for ER networks with z=8. ρ0 is taken as ρc≈0.00747762.
FIG. 12.
Plot of the susceptibility χ, the fluctuation of the order parameter md, as a function of the system size N at κc0.108021. γ/ν¯ is estimated to be 0.69±0.005. Data are obtained for ER networks with z=8. ρ0 is taken as ρc0.00747762.
Plot of the susceptibility χ′, the fluctuation of the order parameter mu, as a function of the system size N at κc≈0.108021. γ′/ν¯′ is estimated to be ≈0.6±0.005. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as ≈0.00747762.
FIG. 13.
Plot of the susceptibility χ, the fluctuation of the order parameter mu, as a function of the system size N at κc0.108021. γ/ν¯ is estimated to be 0.6±0.005. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as 0.00747762.

On the basis of the numerically obtained values β/ν¯0.53 and γ/ν¯0.69, and the theoretical value β=1/3, we estimate ν¯2.179 and γ1.5. Those values are confirmed in Fig. 14 for χNγ/ν¯ versus (κcκ)N1/ν¯, in which the data collapse well with the choices of ν¯2.13±0.1 and γ/ν¯0.69. The measured values of the exponents satisfy the hyperscaling relation (2β+γ)/ν¯0.996 well. Similarly, for κ>κc, on the basis of the numerical values β/ν¯0.164 and γ/ν¯0.6, and the theoretical value β=1/3, we obtain ν¯2.03 and γ1.218. Data for χ for different system sizes collapse well into a single curve with the choices of ν¯=2.13±0.1 and γ/ν¯=0.6 (Fig. 15). These values yield (2β+γ)/ν¯0.910.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 ρc0.00747762 and κc0.108021. We remark that β=β=1/3 is obtained analytically.

Scaling plot of the susceptibility χ for κ<κc in the form χN−γ/ν¯ versus (κ−κc)N1/ν¯, where ν¯≈2.13 and γ≈1.47 are used. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as ≈0.00747762.
FIG. 14.
Scaling plot of the susceptibility χ for κ<κc in the form χNγ/ν¯ versus (κκc)N1/ν¯, where ν¯2.13 and γ1.47 are used. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as 0.00747762.
Scaling plot of the susceptibility χ′ for κ>κc in the form χ′N−γ′/ν¯′ versus (κ−κc)N1/ν¯′, where ν¯′≈2.13 and γ′≈1.28 are used. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as ≈0.00747762.
FIG. 15.
Scaling plot of the susceptibility χ for κ>κc in the form χNγ/ν¯ versus (κκc)N1/ν¯, where ν¯2.13 and γ1.28 are used. Data are obtained for ER networks with mean degree z=8. ρ0=ρc is taken as 0.00747762.

V.. SUMMARY AND DISCUSSION

We investigated the properties of phase transitions in the SWIR model with a finite density ρ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 κ and μ in Eqs. (1) and (2), respectively, serve as control parameters. For convenience, we take κ=μ. We found that for a given network, there exists a critical density of seeds ρc such that for ρ0<ρc, the order parameter, the density of nodes in state R in the absorbing state, increases continuously with the critical exponent β=1/2 as κ is increased up to a transition point κc and then jumps to a finite value, followed by a continuous increase. Accordingly, the order parameter behaves as m(κ)=m0b(κcκ)1/2 for κ<κc, where b is a positive constant. At κc, the order parameter is discontinuous by Δm=mu(κc)m0. 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 β=0. The fluctuation of the order parameter diverges at the transition point κc according to a power law with the exponent γ. For the correlation size exponent ν¯ measured in finite systems, we find that the hyperscaling relation 2β+γ=ν¯ holds reasonably well.

As ρ0 is increased, the jump shrinks and becomes zero at ρc. For ρ0=ρc, the transition becomes continuous. We determined a complete set of critical exponents describing the phase transition at κ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ββγγν¯ν¯
ρ0=1/NSingle seed03
0<ρ0<ρcMultiple seeds1/212
ρ0=ρcMultiple seeds1/31/31.47±0.051.28±0.052.13±0.12.13±0.1

ACKNOWLEDGMENT

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

Appendices

DERIVATION OF THE CRITICAL EXPONENT β AT ρc

Here we introduce an analytical method to determine the critical exponents β=1/3 at ρ0=ρc. It was already noted in Sec. III B that for ρ0=ρc,

F(κc,m0)=dFdm|κc,m0=d2Fdm2|κc,m0=0.
We consider a line of the solution F(κ,m)=0 near (κc,m0) by expanding F(κc+δκ,m0+δm) as
F(κc+δκ,m0+δm)Fκ|κc,m0(δκ)+163Fm3|κc,m0(δm)3+=0
where only nonzero terms are considered. Since δκ and (δm)3 are two lowest terms in Eq. (A2) and their coefficients have the opposite sign to each other, δm(δκ)1/3 when δκ1. Thus for both cases of κ<(>)κc, the critical exponents β=β=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).

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,