General protocol for predicting outbreaks of infectious diseases in social networks

Epidemic spreading on social networks with quenched connections is strongly influenced by dynamic correlations between connected nodes, posing theoretical challenges in predicting outbreaks of infectious diseases. The quenched connections introduce dynamic correlations, indicating that the infection of one node increases the likelihood of infection among its neighboring nodes. These dynamic correlations pose significant difficulties in developing comprehensive theories for threshold determination. Determining the precise epidemic threshold is pivotal for diseases control. In this study, we propose a general protocol for accurately determining epidemic thresholds by introducing a new set of fundamental conditions, where the number of connections between individuals of each type remains constant in the stationary state, and by devising a rescaling method for infection rates. Our general protocol is applicable to diverse epidemic models, regardless of the number of stages and transmission modes. To validate our protocol’s effectiveness, we apply it to two widely recognized standard models, the susceptible–infected–recovered-susceptible model and the contact process model, both of which have eluded precise threshold determination using existing sophisticated theories. Our results offer essential tools to enhance disease control strategies and preparedness in an ever-evolving landscape of infectious diseases.

We can express S kjn , I kjn , and R kjn using binomial distributions involving p k , q k , and v k for j infected neighbors, as well as w k , x k , and y k for n recovered neighbors where C k j represents a combination of k links taken j at a time without repetition and ŵk = w k /(1 − p k ), xk = x k /(1 − q k ), and ŷk = y k /(1 − v k ) denote the conditional probabilities of a link being connected to a recovered neighbor given that the link is not connected to an infected neighbor.Subsequently, S k , I k , and R k are calculated by summing S kjn , I kjn , and R kjn across the range of j and n.

B. Total infection, recovery, and loss-of-immunity rates
The total infection rate (σ), the total recovery rate (π), and the total loss-of-immunity rate (φ) in the population are defined as follows: In the steady state, all of these total rates become equal, σ = π = φ, as indicated by Eq. ( 4) in the main text.

C. Bond-detailed-balance equations
From Table I in the main text, we derive six equations of the detailed balance conditions by summing the product of each contribution to ∆Π •• and the corresponding probability.− kjn rI kjn • (k − j − n) + kjn hR kjn • j = 0, (10) ∆Π IR = kjn λS kjn j • n ∆Π RS = − kjn λS kjn j • n + kjn rI kjn • (k − j − n) It is important to note that Cai et al. employed ∆σ = ∆π = 0 to derive Eq. (8), for the SIS model [12].However, in the context of the SIRS model, Eq. ( 8) cannot be derived from ∆σ = ∆π = 0; it naturally emerges as one of the BDB conditions.Consequently, the BDB conditions play a crucial role in accounting for dynamic correlations.

D. Identities for the number of bonds with different types of nodes
The total number of infected neighbors of all susceptible nodes equals the total number of susceptible neighbors of all infected nodes, Π SI = Π IS , and the same applies for other bonds involving different types of nodes.
where we used pk = 1 − p k and wk = 1 − ŵk .Based on this, we can derive the following relation.
By utilizing this relation and the identity Π IR = Π RI , Eq. (15) above, we arrive at Eq. ( 9) as presented in the main text.We transform Eq. ( 9) in the main text into an equation involving only I k , which allows us to determine the threshold.We express the second and third terms Eq. ( 9) in the main text with I k by using the similar procedure as described for Eq. ( 16).
Similarly, we proceed with the calculation as follows: where we used ŵk = w k /(1 − p k ).
For expressing the first term in terms of I k , we need to use the identities Π SI = Π IS , Eq. (13).By utilizing Eqs. ( 8) and (13), the first term of Eq. ( 9) in the main text can be expressed as: Substituting Eqs. ( 18), (19), and (20) into Eq.( 9) in the main text, we can rewrite it as Eq. ( 10) in the main text:

F. Derivation of epidemic threshold
We define the left-hand side of Eq. (10) in the main text as f (p) and represent f (p) compactly as f (p) = L k I k + H k kI k by introducing the definitions of H and L: By differentiating f (p) and using the relations we derive the equation for λ c : leading to the epidemic threshold λc = λ c /r of the SIRS model as .

A. Scaled rates for simulations
In the main text, we derived the expression for λSIRS For simulation purposes, it is more convenient to express λSIRS c in terms of scaled rates, which satisfy the conservation relation λ + h + r = 1.In the following subsections, we will use λ, r, and h as the probabilities for the infection, recovery, and loss-of-immunity processes, respectively.
First, we can simplify the equation by eliminating r from Eq. (24) using r = 1 − λ c − h, which yields the following equation: Now, we define the coefficients as follows: Given that k 2 k and h ≤ 1, all three coefficients are positive.Consequently, we obtain a second-order equation in terms of λ c : The solutions for λ c are then given by: By using the relation λHMF c = k / k 2 , we can rewrite 4xz/y 2 as 4(1 − h) λHMF (30) Solving Eq. (30), we obtain λ c (r) as follows: However, the positive solution violates the condition λ c = λ SI c = 0 for r = 0.By considering the negative solution, we can express λ c (r) as: For the square root to yield a real value, r should satisfy the following condition: Additionally, from h = 1 − λ c − r ≥ 0, we require: We can rephrase the inequality of Eq. (34) in terms of r as: Consequently, r should satisfy one of the two conditions, either Eq.(33) or Eq.(35), depending on the magnitude of λSIS c .However, for λSIS that the inequality of Eq. ( 35) must be satisfied by r, restricting the maximum r (r max ) as r max = 1 − λSIS c .It is important to note that λ c = λ SI c = 0 when r = 0, as expected.As r increases, λ c (r) monotonically increases from 0 to the maximum value r max , above which h becomes negative.

B. The uncorrelated configuration (UC) algorithm
for the generation of networks In the quenched scale-free network (SFN) characterized by a degree distribution P (k) ∼ k −γ , the network is designed with predetermined minimum and maximum degrees, represented as k min and k max , respectively.The procedure for constructing these networks adheres to the algorithm outlined in Refs.[1,2].
The number of nodes with a degree of k, denoted as N k , is determined through the following deterministic calculation: . This is done while ensuring that the total number of nodes, denoted as N , adheres to the equation N = kmax k=kmin N k .Here, the notation int[x] signifies the integer part of the real number x.To allocate a degree of k to N k nodes, a random selection is made from the pool of nodes that have not yet been assigned degrees.Subsequently, the selected N k nodes are assigned the degree k.By repeating these steps for all integer values of k within the range [k min , k max ], each node is allocated a specific degree.This process generates the degree sequence K N = k 1 , ..., k N for all N nodes.
Finally, to create a quenched network, a procedure is employed where two nodes are chosen at random.The crucial requirement for this selection is that the two chosen nodes are not already linked and that their current degrees do not surpass their designated degrees as per the degree sequence K N .If these conditions are met, the two nodes are linked together.However, if the conditions are not satisfied, an alternative pair of nodes that adheres to the conditions is randomly chosen.This process of linking nodes continues until all nodes are connected to one another according to the predefined degree sequence K N .
Since the algorithm used is deterministic, all the generated quenched networks follow the same degree sequence K N .The sole variation among these networks lies in the connections formed between nodes.The process of randomly rewiring the connections of N nodes only results in the rearrangement of linked nodes, and it doesn't impact the distribution P (k).As a result, networks characterized by the same degree sequence K N yield identical network-level averages for any observable X, expressed as X = k XP (k).Hence, during our simulations, we utilize a single network and compute ensemble averages across multiple samples, commencing from distinct initial distributions yet sharing the same initial density of I and S.

C.
The SIRS on a quenched SFN of γ > 4 To validate the expressions for λ c (h) and λ c (r) as given by Eqs. ( 29) and (32), we perform simulations of the SIRS model on quenched SFNs.
We construct a scale-free network (SFN) with a powerlaw degree distribution characterized by γ = 5, N = 10 7 , k min = 5, and k max = 100.The uncorrelated configuration (UC) algorithm detailed earlier is used to generate this network configuration.The resulting network has an average degree k = 6.0942 and an average squareddegree k 2 = 41.547.For this particular network configuration, we have the following values for λHMF With these parameters in place, we proceed with simulations to analyze and verify the expressions for λ c (h) and λ c (r).
For this quenched SFN, simulations are conducted using a sequential updating scheme.The process involves creating a list of infected and recovered nodes, followed by the random selection of a node from the list.If the chosen node is infected, it transitions to the recovered state with a probability r (I → R).On the other hand, if the chosen node is susceptible, all susceptible nodes linked to the infected node may become infected with a probability λ (S → I).If a recovered node is selected, it transitions to the susceptible state with a probability h (R → S).The list is updated in response to changes brought about by the S → I and R → S processes.Importantly, the I → R process does not alter the list.
After updating the list, the Monte Carlo time t is incremented by dt = 1/M (t), where M (t) represents the total number of infected and recovered nodes, i.e., the size of the list, before the update.The initial conditions are established by randomly infecting I(0) nodes in the entirely susceptible population of size N , which results in S(0) = N − I(0) and R(0) = 0.For a specific initial condition defined by a given I(0) value, we run simulations until a predefined time t max and calculate the averages of the quantities S(t), I(t), and R(t).
Initially, we outline the process used to estimate λ c for h = 0.5 and r = 0.5.Subsequently, we compare these estimates with the theoretical values of λ c (h) and λ c (r), which are provided by Eqs.( 29) and (32), corresponding to Eqs.( 14) and (15) in the main text.Finally, we create two phase diagrams, namely λ c (h) versus h and λ c (r) versus r, based on the obtained estimates for h in the range of [0.1, 0.9] and r in the range of [0.1, 0.7].
(1) λ c (h) for h = 0.5 To investigate the off-critical scaling behavior, we conduct simulations up to t max = 10 4 using random initial conditions of S(0) = I(0) = 0.5.The quantity I(∆, t) is averaged over approximately 50 samples, and simulations are carried out for various values of λ ranging from 0.081 to 0.089.To obtain the steady-state value I(∆), we perform a time average within an interval of [τ, t max ], where τ denotes the time after which I(∆, t) stabilizes to the steady state.
We investigate the power-law decay of I(∆) with respect to ∆ as ∆ β for ∆ → 0. To determine λ c , we examine the linearity of I(∆) by varying a tentative value λ * c in a double logarithmic plot of I against ∆ (where ∆ = λ−λ * c ).If we adjust λ * c to be either larger or smaller than the true λ c , the plot of I(∆) displays nonlinearity as ∆ approaches zero.As a result, we estimate λ c by identifying the λ * c value that results in a linear I(∆) plot and calculate the slope of the least-square fitting line to determine the exponent β.
In Fig. 1(a), the double logarithmic plot of I against ∆ at λ * c = 0.0796 is depicted, revealing the clear linear decay of I as ∆ approaches zero.This observation leads us to estimate λ c as 0.0796 for h = 0.5.The theoretical value of λ c (h), denoted as λ th c (h), for h = 0.5 is λ th c (h = 0.5) = 0.07861.While our estimate is slightly higher than λ th c (h), the error remains minimal at only 1.253%.By applying the least-square fitting to the data in Fig. 1(a), we calculate the slope β = 0.9970 (9).The value in the parentheses signifies the error associated with the last digit of the fitting result.To ensure the stability of the slope estimation, we compute the effective slope β eff between consecutive data points, and the results are plotted in Fig. 1(b).As evident from the figure, β eff (∆) saturates and remains consistently close to 1, thereby confirming the precision of the estimated λ c .
To further verify the accuracy of λ c = 0.0796, we approach the estimation through an alternative method.We examine the critical decay of I(t) ∼ t −α at λ c = 0.0796.In this procedure, simulations are conducted up to t max = 10 3 , commencing from random initial conditions with S(0) = 0.2 and I(0) = 0.8, and the average I(t) is computed over 1000 samples.(2) λ c for r = 0.5 Similarly to the h = 0.5 case, we adopt a similar approach to estimate λ c for r = 0.5.This involves examining the off-critical scaling behavior of I(∆) and subsequently confirming the accuracy of the estimated λ c through the critical decay of I(t).
The simulations are initiated with various λ values ranging from 0.0975 to 0.101.The simulations begin with random initial conditions of S(0) = I(0) = 0.5 and continue up to t max = 3 × 10 4 , the precise value depending on λ.The quantity I(∆) is averaged over a span of up to 200 samples.
By fine-tuning λ * c in the double logarithmic plot of I(∆), the optimal linearity of I(∆) is achieved at λ * c = 0.09646.Fig. 3 depicts the double logarithmic plot of I(∆) and β eff (∆) for this value of λ * c .By analyzing the slope of the least-square fitting line in Fig. 3(a), we estimate λ c = 0.09646 and β = 1.00 (1).Furthermore, Fig. 3(b) demonstrates that β eff (∆) remains consistently close to 1, providing additional confirmation of the precision of the estimated λ c .
To confirm the power-law decay at the estimated λ c , simulations were conducted at λ c = 0.09646 with a runtime of up to t max = 500.These simulations began with random initial conditions of S(0) = I(0) = 0.5, and I(t) was averaged over 1000 samples.As depicted in Fig. 4(a), I(t) exhibited a power-law decay pattern after a brief transition period.Utilizing the scaling plot showcased in Fig. 4(b), an estimate of α = 0.992 was derived.
Based on the findings from the analyses of off-critical and critical scaling behaviors of I, it can be concluded that the estimates for r = 0.5 are λ c = 0.09646 and α = β = 1 for a network with γ = 5.The theoretical value of λ c (r) (λ th c (r)) for r = 0.5 is λ th c (r = 0.5) = 0.09497.Thus, the presented estimate is only 1.583% higher than the value of λ th c (r).
(3) The phase diagrams In a manner analogous to the estimation processes for h = 0.5 and r = 0.5, we perform estimations for λ c (h) by incrementally increasing h from 0.1 to 0.9 in increments of 0.1.Additionally, we estimate λ c (r) by increasing r from 0.1 to 0.7 in increments of 0.1.It is important to note that r is constrained by the maximum value of r (r max ) given by r max = 1 − λSIS c = 0.8281, as determined by λSIS c of Eq. (36).Values of r exceeding this limit result in h ≤ 0.
We have constructed the phase diagrams using the estimates of λ c (h) and λ c (r), as shown in the plot of λ c (h) versus h (Fig. 5(a)) and the plot of λ c (r) versus r (Fig. 5(b)).In both diagrams, one region where λ > λ c corresponds to the active/endemic phase, while the other region where λ < λ c corresponds to the absorbing/disease-free phase.Notably, all the esti- mates are in close proximity to the theoretical lines obtained through Eqs. ( 29) and (32), which correspond to Eqs. ( 14) and (15) in the main text.The errors in the estimates from the theoretical values are minimal, measuring less than 1.3% for λ(h) and less than 1.584% for λ(r), respectively.Consequently, all the estimated λ c values demonstrate excellent agreement with the theoretical predictions.
For the SIS model, both the HMF theory and the mean-field Langevin equation predict β HMF = 1 when γ > 4, and β HMF = 1/(γ − 3) for 3 < γ < 4 and α = β HMF [7][8][9].Applying the HMF theory to the SIRS model yielded the same results as those observed in the SIS model, indicating that the state R does not alter the mean-field exponents.Remarkably, our simulation results of α = β = 1 for γ = 5 align well with the predictions of the HMF theory for γ > 4.

D.
The SIRS on a quenched SFN of 3 < γ < 4 We generate a scale-free network with γ = 3.5, comprising N = 10 8 nodes.The minimum and maximum degree values for this network are set as k min = 5 and k max = 5724, respectively.The network construction procedure results in average degree values k = 7.5751 and k 2 = 100.1724.Based on these values, we calculate the corresponding λHMF For h = 0.5 and 0.9, we follow a similar process to estimate λ c (h) as previously done for γ = 5.
For h = 0.5, we increase λ by 5 × 10 −4 within the range of 0.0405 to 0.0425.Subsequently, we compute I(∆, t) up to t max = 10 4 and then average I over up to 160 samples.The simulations are initiated with random initial conditions, where both I(0) and S(0) are set to 0.5.
For λ c = 0.0381, we perform simulations using 180 samples, with the initial conditions set to I(0) = 0.3 and S(0) = 0.7.The simulation time is limited to t max = 100 since samples with I(t) = 0 start to emerge beyond this point due to finite-size effects.We estimate α = 1.79 from the scaling plot I(t)t α , which aligns closely with the β estimate.The patterns observed in the plots of I(t) and I(t)t α vs. t mirror those in Figs. 2 and 4, hence we omit these plots for brevity.Utilizing the values of λHMF c and λSIS c from Eq. (37), we derive λ th c (h) = 0.03923 for h = 0.5 using Eq. ( 29).This theoretical value is 2.76% higher than the estimated λ c = 0.0381.
For h = 0.9, we increase λ by 10 −4 across the range from 0.0077 to 0.0081.Subsequently, we measure I(∆, t) up to t max = 10 4 and average I over up to 200 samples, with the initial conditions set to I(0) = S(0) = 0.5.
Upon analyzing the data, we discern that I(∆) demonstrates the highest linearity at λ * c = 0.007235, with β eff (∆) remaining stable at approximately 1.9.Fig. 6(c) and 6(d) display the double logarithmic plots of I(∆) and β eff (∆) corresponding to λ * c = 0.007235.Hence, we estimate λ c = 0.007235 and β = 1.9051(5) using the slope of the least-square fitting line in Fig. 6(c).We conduct simulations with 270 samples, beginning with random initial conditions of I(0) = 1.0, and I(t) is measured up to t max = 700.We estimate α = 1.93 from the scaling plot I(t)t α , which again aligns well with the β estimate.Similar to the h = 0.5 case, we derive λ th c (h) = 0.00762 for h = 0.9 using the values of λHMF c and λSIS c from Eq. (37).This theoretical value is approximately 5% higher than the simulation value of λ c = 0.007235.
The estimated values of λ c for h = 0.5 and 0.9 agree well with the predicted values λ th (h) from Eqs. ( 29) and (32), which are presented in Eqs.( 14) and (15) of the main text.Additionally, it is worth noting that the estimate of β for h = 0.9 is even closer to β HMF than the estimate for h = 0.5.Here, β HMF = 1/(γ − 3) = 2 for the case of γ = 3.5.This phenomenon is a result of the SIRS model becoming more akin to the SIS model as the value of h increases, owing to the relationship described in Eq. ( 24), which is presented as Eq.( 12) in the main text.Considering the fact that larger networks are anticipated to yield β values that are closer to β HMF , we can be more confident in asserting that the simulation results for γ = 3.5 validate the predictions for β and α that the HMF theory puts forward, i.e. β HMF = 1/(γ − 3) and α = β HMF for 3 < γ < 4.
Through these simulation results, we come to the conclusion that while dynamic correlations do impact the thresholds, they do not exert the same influence on the critical exponents in quenched SFNs.

E. Verification of bond-detailed-balance conditions
We verify the equations of bond-detailed-balance conditions, which are represented by Eqs. ( 7) through (12), through simulations.Let us define n s as the count of susceptible nodes, given by n s = k − j − n.This signifies the number of susceptible nodes among the k neighbors for a given combination of j and n.By rearranging Eqs. ( 7) through (12), we obtain six equalities to be verified, as presented below.
∆Π SS = 0 : All terms on both sides of the six equalities are calculated up to t max = 500 and then averaged over 10 samples.The simulations are conducted on the same quenched SFNs with parameters γ = 5, N = 10 7 , k min = 5, and k max = 100.Specifically, the simulations are carried out for a case with λ = 0.087 and h = 0.5.The estimated value of λ c for h = 0.5 is λ c = 0.0796 (as shown in Fig. 1), making the considered λ to be close to λ c , leading to ∆ = λ − λ c = 7.4 × 10 −3 .Fig. 7 presents six plots depicting the left-handside and the right-hand-side of each equality from Eq. (38).For instance, Fig. 7(a) displays the comparison kjn λS kjn jn s = kjn hR kjn n s , where ∆Π SS = 0.As observed, both sides of each equality demonstrate a decay in time and eventually converge to each other within the stationary state, which occurs after around t ≈ 100.
Furthermore, we have performed additional simulations for cases where λ is significantly different from λ c , specifically λ = 0.15 and 0.3.The results of these simulations also show the alignment of both sides of the six equalities in the stationary state, similar to the pattern exhibited in Fig. 7. Therefore, based on the simulation results that encompass scenarios near and far from the threshold, it is convincingly demonstrated that the bonddetailed-balance conditions persist within the stationary state of the endemic phase in the SIRS model.

III. THE CONTACT PROCESS
The contact process (CP) stands as the archetype among nonequilibrium models that belong to the directed percolation universality class, exhibiting the absorbing phase transition on lattice structures [3][4][5].
In the CP model, individuals are divided into two categories: susceptible (S) and infected (I).The processes of infection and recovery take place in the following manner: (i) A node is randomly selected.Should the chosen node be infected, a neighbor of the selected node is subsequently chosen at random.(ii) If the chosen neighbor is susceptible, it becomes infected with a rate denoted as λ.
(iii) The selected node then undergoes recovery at a rate of r.It is worth noting that a key difference between the CP model and the SIS model is that in the CP model, a single neighbor is randomly selected for the infection process, which takes place only if the chosen neighbor is susceptible.This random selection constitutes a crucial characteristic of the CP model, causing its threshold and critical exponents to deviate from those of the SIS model on complex networks [10].
On complex networks, the Contact Process exhibits an absorbing phase transition at a threshold λc = λ c /r.In the context of heterogeneous mean-field theory, λc is predicted to be equal to 1 for networks with γ > 2 [10].However, the accuracy of HMF theory in predicting the CP model's threshold λc on quenched networks remains uncertain.Similar to the SIR, SIS, and SIRS models, the dynamic correlations between connected nodes are also expected to play a crucial role in determining the CP model's threshold λc on quenched networks.To address this, the Heterogeneous Pair-Approximation (HPA) was introduced for the CP model [11].The HPA focuses on pairs of nodes connected to each other and formulates rate equations for the density of total pairs consisting of infected-susceptible nodes, as opposed to the rate equations for infected nodes directly.However, the HPA only provides a self-consistent equation for λc in the steady state, and numerical methods are used to determine λc from this equation [11].The self-consistent equation from the HPA is expressed as follows: where P (k) represents the degree distribution.The authors of Ref. [11] compared the results of Monte Carlo simulations with the numerically obtained values of λc from Eq. (39) for quenched scale-free networks (SFNs) with γ values ranging from 2.3 to 3.5, showing a high degree of agreement.The HPA succeeds in capturing the impact of dynamic correlations, although it falls short of providing an explicit expression for the CP model's threshold.
For the Contact Process, the introduction of random selection of one neighbor poses a theoretical challenge in determining λc in comparison to the SIS and SIRS models.Hence, it is crucial to devise an approach that incorporates the dynamic correlations arising from this random selection process, in a manner similar to the utilization of the probability p k in the SIS and SIRS models.
Additionally, there exists uncertainty regarding whether the CP demonstrates a transition on quenched SFNs when γ < 3.In cases where k 2 diverges and a minority of hubs control a significant portion of total links, the network structure is highly irregular.This situation leads to the majority of infections occurring at hubs and their immediate neighbors.This intricate network structure emphasizes the heightened importance of dynamic correlations for quenched SFNs with γ < 3, rendering the predictions of the HMF theory for λc less reliable.
As is evident in the SIS and SIRS models, the value of λc is influenced by both the first and second moments of node degree.Given this, it is reasonable to anticipate that this dependency is also valid for the CP model.Thus, it becomes imperative to derive an accurate expression for the threshold λc of the CP model on quenched networks to address the challenges outlined above.
To derive an explicit expression for λc in the context of the CP model on quenched networks, our approach involves initially devising a means to represent the dynamic correlations stemming from the random selection process using the previously defined probability p k in the SIRS model.Subsequently, we apply the methodology outlined in the main text to attain our desired expression.
Similar to the approach employed for the SIRS model in the main text, we introduce the variables S kj and I kj which represent the number of nodes connected to j infected neighbors within the groups of S k and I k nodes, respectively.In light of the principle of maximum entropy, we assume uniform values for p k across all degrees, denoted as p k = p ∀k. Taking into account the impact of randomly selecting one neighbor, we can formulate the rate equation for I k (t) in the CP model through the Heterogeneous Mean-Field approximation as follows: where 1/k in the first line arises due to the stochastic selection of one neighbor from a pool of k neighbors and the second line makes use of the property P (k |k) = k P (k )/ k for uncorrelated networks.This rate equation effectively encapsulates the influence of dynamic correlations, stemming from the random neighbor selection, on the evolution of infected nodes in the CP model.
To incorporate the term p k into Eq.( 40), similar to what was done for the SIRS model, we aim to determine the relationship between p k and the second term in the final line of Eq. (40).In the stationary state, we arrive at a self-consistent equation for I k by setting dI k dt = 0 and employing S k + I k = N k : which yields the expression for I k : with the last equality valid near λc due to I k = 0 at λc .By utilizing the expression for I k from Eq. ( 42), the expression for p = p k can be deduced from the HMF theory: The first equality represents the expression of p as given by the HMF theory for the SIS and SIRS models.On the other hand, the summation in the second line of Eq. ( 40) can be redefined using the p expression from Eq. (43): This allows us to rewrite the rate equation dI k dt from Equation (40) as: which achieves our objective of incorporating p.This equation provides the stationary-state relationship: The equation for I k involving p is given by: As λ approaches λ c , Eq. ( 45) and ( 46) become increasingly accurate.Eq. ( 45) aligns with the corresponding equations of the SIRS model (Eq.( 2) in the main text) and the SIS model of Ref. [12] by introducing a scaled infection rate denoted as Λ.This scaled infection rate, defined as Λ = λ k / k 2 , allows for a seamless application of our method to the CP model.
The total infection rate (σ) and the total recovery rate (π) are given by: In Table 1, we summarize all possible changes of the numbers of all types of bonds in time interval dt, along with the corresponding probabilities.From Table I, we deduce three equations based on the bond-detailed-balance conditions.These equations result from summing the product of each contribution to ∆Π ab (where a and b can be either I or S) and the corresponding probability.By setting ∆Π ab = 0, we obtain the following equations: where The three equations derived from the bond-detailedbalance conditions can be simplified into a single equation, which corresponds to Eq. (50).To evaluate each term in Eq. ( 50), we use the identity that the total number of infected neighbors of all susceptible nodes is equal to the total number of susceptible neighbors of all infected nodes, Π SI = Π IS .This leads to another equation: The second term in Eq. ( 50) can be obtained from I kj j = kI k − S kj j, utilizing Eq. ( 52).Using the expressions for S kj and p = 1 − p, we calculate as follows: The first term in Eq. ( 50) is obtained as: where the relation rI k = ΛkpS k from Eq. ( 46) is used in the last line.Combining Eqs. ( 53) and (54), we can express Eq.(50) as: We finally arrive at the equation for I k by utilizing Eq. (46): To determine the threshold λc (which equals λ c /r), we define the left-hand side of Eq. (56) as f (p) and rearrange it as follows: The function f (p) is equal to 0 at p = 0 and becomes negative at p = 1.As a result, the slope of f (p) at p = 0 (df (p)/dp| p=0 ) should be positive in order to have a nontrivial solution for p ∈ (0, 1), and it should approach zero as λ approaches λc .Thus, the condition for the transition to occur is df (p)/dp| p=0 = 0, similar to the cases of the SIS and SIRS models.By differentiating I k (p) of Eq. ( 47) at p = 0, we obtain: (58) With these expressions, we can derive the equation for λc from df (p)/dp| p=0 : Solving Eq. ( 59), we finally obtain the threshold λc for the CP model: Unlike the SIS and SIRS models, the threshold of the CP model remains finite for γ > 2, even in cases where k 2 diverges.For 2 < γ < 3, k 2 diverges as the network size N approaches infinity, leading to λc = 1 for N → ∞.While the HMF theory for the CP predicts λc = 1 for scale-free networks with γ > 2 [10], our result aligns with λc = 1 only for 2 < γ < 3. On quenched SFNs with γ > 3 and random networks following a Poissonian degree distribution P (k), λc exceeds 1 and varies with both k and k 2 .
Furthermore, the threshold λc of the CP model on a regular network with a degree of k 0 and a degree distribution of P (k) = δ k,k0 is known to be [12]: When substituting k 2 = k 2 0 and k = k 0 for the regular network into Eq.(60), we correctly reproduce the threshold of the CP model on regular networks.

IV. SIMULATION RESULTS OF THE CP ON QUENCHED SCALE-FREE NETWORKS
In the derivation of the threshold of the CP model, a crucial approximation is the representation of dynamic correlations arising from the random selection of one neighbor by λp k / k 2 in Eq. ( 45), which corresponds to scaling λ with k / k 2 as Λ = λ k / k 2 .Subsequently, we obtain an equation for dI k (t)/dt that is analogous to the SIS model of Ref. [12] and the SIRS models discussed in the main text.This simplifies the subsequent calculations.To validate the accuracy of this approximation, it's essential to verify one of the bonddetailed-balance conditions, specifically ∆Π II = 0 of Eq. (50), i.e., Λ k,j kj 2 S kj = k,j jI kj , through simulations.Here, Λ = λ k / k 2 .
To verify the threshold λc as given by Eq. (60) (referred to as Eq. ( 16) in the main text), we conducted simulations on quenched SFNs with varying values of γ ranging from 2.25 to 10, and N fixed at 10 7 , along with k min = 5.The upper limit of the degree, k max , increases from 27 for γ = 10 to 5 × 10 3 for γ = 2.25.
We estimated λc for each value of γ using methods similar to those used for the SIRS model.As examples, simulation results for γ = 5 and 2.5 are depicted in Figs .9 and 10.At the estimated λc , we confirmed the power-law decay of I(t) with α = β.For the CP model on quenched SFNs, the exponent β HMF follows β HMF = 1 for γ > 3 and β HMF = 1/(γ − 2) for 2 < γ < 3, yielding β HMF = 2 for γ = 2.5 and β HMF = 1 for γ = 5 [9][10][11].Comparing β HMF with our estimates, we find excellent agreement for γ = 5, but discrepancies for γ = 2.5 as evident in Figs. 9 and 10.Similarly, we also observe underestimation of β for γ = 2.25 and 2.8.The measurement of scaling exponents often encounters substantial finite-size effects for γ < 3, contributing to the underestimation of β in this range [9,10].By estimating λc while incrementing γ from 2.25 to 10, we map out the phase diagram illustrated in Fig. 11.In this phase diagram, each theoretical λc is computed using the k and k 2 specific to the corresponding network, and plugged into Eq.( 60).The errors in the estimates are presented as percentages, calculated using the formula 100 × (estimate − λc )/ λc .As evident in the inset of Figure 11, errors sharply increase near γ = 3, reflecting the significant structural change of the network at this point.For other values of γ, errors remain quite small, ranging from 2% to 4%.Thus, the simulation outcomes robustly endorse the theoretical threshold of Eq. ( 60) and the validity of our approach for the CP model.

Fig. 2 (
a) presents a double logarithmic plot of I(t) against t.As clearly depicted, I(t) exhibits a decay consistent with t −α behavior following an initial transient time.To ascertain the decay exponent α, we analyze the

FIG. 2 :
FIG. 2: The critical decay for h = 0.5 : (a) Plot of the critical decay of I(t) at λ = 0.0796.The red straight line corresponds to the fitting line with a slope of unity.(b) The scaling plot of f (t, α) = I(t)t α with α = 1.014.The horizontal line indicates the average value of f (t, α) in the scaling region.
scaling plot of f (t, α) = I(t)t α , where α is varied.An accurate choice of α results in a constant value of f (t, α) within a scaling region governed by I(t) ∼ t −α .As shown in Fig.2(b), we observe that f (t, α) achieves a flat plateau for α = 1.014, which reinforces the precision of the estimated λ c = 0.0796.Consequently, we confidently determine the estimates for h = 0.5 as λ c = 0.0796 and α = β = 1 for γ = 5.
FIG. 7: The verification of ∆Π•• = 0 for h = 0.5 and λ = 0.087 (∆ = 7.4 × 10 −3 ) on the quenched SFN with γ = 5, N = 10 7 , kmin = 5, and kmax = 100 : Each panel displays both the left-hand and right-hand sides of the equalities from Eq. (38).These values are presented as solid and dashed lines, respectively, as specified in the legend.The y-axis title indicates the corresponding bond change.Across all panels, it is evident that the values on both sides of each equality align perfectly, thus serving as confirmation of the bond-detailedbalance conditions outlined in Eq. (38).

TABLE I :
The bond changes by all possible events in a time interval dt of the CP on a quenched network.Λ represents the scaled infection rate defined as Λ = <k> <k 2 > λ, and Σ = σ + π.Event S kj → I kj I kj → S kj Probability ΛS kj j/Σ rI kj /Σ ∆ΠSS