Heterogeneous pair-approximation for the contact process on complex networks

Recent works have shown that the contact process running on the top of highly heterogeneous random networks is described by the heterogeneous mean-field theory. However, some important aspects such as the transition point and strong corrections to the finite-size scaling observed in simulations are not quantitatively reproduced in this theory. We develop a heterogeneous pair-approximation, the simplest mean-field approach that takes into account dynamical correlations, for the contact process. The transition points obtained in this theory are in very good agreement with simulations. The proximity with a simple homogeneous pair-approximation is elicited showing that the transition point in successive homogeneous cluster approximations moves away from the simulation results. We show that the critical exponents of the heterogeneous pair-approximation in the infinite-size limit are the same as those of the one-vertex theory. However, excellent matches with simulations, for a wide range of network sizes, are obtained when the sub-leading finite-size corrections given by the new theory are explicitly taken into account. The present approach can be suited to dynamical processes on networks in general providing a profitable strategy to analytically assess and fine-tune theoretical corrections.

The heterogeneous mean-field (HMF) approach for dynamical processes on complex networks has become widespread in the last years.This theory, formerly conceived to investigate the SIS dynamics on complex networks [16], assumes that the number of connections of a vertex (the vertex degree) is the quantity relevant to determine its state, neglects all dynamical correlations as well as the actual structure of the network.On the other hand, the quenched mean-field (QMF) theory [3,17] still neglects dynamical correlations but the actual quenched structure of the network is explicitly taken into account by means of the adjacency matrix A ij that contains the complete information of the connection among vertices [18].More recently, semi-analytic methods including dynamical fluctuations [5,6] and heterogeneous pair-approximations [8,19,20,7] have appeared as more accurate alternatives to HMF theory.
The CP [21] is the simplest reaction-diffusion process exhibiting a transition between an active and a frozen (absorbing) phase [22].The CP dynamics investigated in the present work is defined as follows [22]: A vertex i of an arbitrary unweighted graph can be occupied (σ i = 1) or empty (σ i = 0).At a rate λ, an occupied vertex tries to create an offspring in a randomly chosen nearest-neighbor, what happens only if it is empty.An occupied vertex spontaneously disappears at rate 1 (this rate fixes the time unit).Notice that in the SIS dynamics an occupied vertex creates ('infects' in the epidemiological jargon) an offspring in each empty nearest neighbor at rate λ.Even being equivalent for strictly homogeneous graphs (k i ≡ k ∀ i), these models are very different for heterogeneous substrates (see discussion in Ref. [23]).However, in both models the creation of particles is a catalytic process occurring exclusively in pairs of empty-occupied vertices, implying that the state devoid from particles is a fixed point of the dynamics and is called absorbing state.
After an intense discussion [13,12,24,25,26], the HMF theory showed up as the best available approach to describe scaling exponents associated to the phase transition of the CP on networks [27].However, some important questions remained unanswered.The transition point λ c = 1 predicted by the HMF theory [24] does not reflect the dependence on the degree distribution observed in simulations [24,27,11].Most intriguingly, it was observed a good accordance between simulations and a heuristic modification of the strictly homogeneous pair-approximation (HPA) (see Ref. [22] for a review) where the fixed vertex degree is replaced by the average degree of the network [11,27]: Finally, sub-leading corrections to the finite-size scaling, undetected by the one-vertex HMF theory, are quantitatively relevant for the analysis of highly heterogeneous networks, for which deviations from the theoretical finite-size scaling exponents were reported [27].Dynamical correlations is an important factor to ascertain the accuracy of the analytical results.The simplest way to explicitly consider dynamical correlations is by means of a pair-approximation [22].In this paper, we present a pair HMF approximation for the CP on heterogeneous networks.We have found that this theory yields great improvements in relation to the one-vertex counterpart but, however, may be farther from the simulation thresholds than Eq. ( 1).This apparent contradiction is solved showing that the higher-order homogeneous cluster approximations overestimate the actual transition point implying that the proximity is only a coincidence.We also show that pair HMF theory yields the same critical exponents as the one-vertex HMF theory, but with different corrections to the scaling.These corrections allow an almost perfect match with simulations constituting a great improvement in relation to the one-vertex mean-field theories [15,14].
The paper is organized as follows: Pair HMF theory is proposed and the transcendental equation that gives the transition points is derived in section 2. The numerical analyses of the thresholds and the comparisons with quasi-stationary simulations are presented in section 3. The critical exponents are analytically determined and compared with simulations in section 4. Our concluding remarks are drawn in section 5.

Pair HMF theory
In this section, we develop the pair HMF theory where the evolution of the system is given by the average behavior of vertices with the same degree.So, let us introduce the notation based on Ref. [7]: [A k ] is the probability that a vertex of degree k is in the state A; [A k B k ] is the probability that a vertex of degree k in state A is connected to a vertex of degree k in state B; [A k B k C k ] is the generalization to three vertices such that the pairs [A k B k ] and [B k C k ] are connected through a node of degree k and so forth.An occupied state is represented by 1 and an empty one by 0. The pair-approximation carried out hereafter uses the following notation: Obviously, we have that ψ kk = ψ k k , ω kk = ω k k , and φ kk = φk k .Independently of the dynamical rules, the following closure relations can be derived from simple probabilistic reasonings: The master equation for the probability that a vertex with degree k is occupied takes the form where the conditional probability P (k |k), which gives the probability that a vertex of degree k is connected to a vertex of degree k , weighs the connectivity between compartments of degrees k and k .The first term of Eq. ( 3) represents the spontaneous annihilation and the second term reckons the creation in a vertex of degree k due to its nearest neighbors.The dynamical equation for φ kk is The first term represents the annihilation in the vertex of degree k , the second one includes the creation in the vertex of degree k due to the connection with the neighbor of degree k and the third one is due to the annihilation of the vertex with degree k.These terms represent the reactions inside pairs with degrees k and k , that create or destroy a configuration [0 k , 1 k ].The fourth and fifth terms represent changes due to creation in vertices with degree k and k, respectively, due to all their neighbors except the link between the vertices of the pair itself, which is explicitly included in the second term.
The one-vertex mean-field equation proposed in Ref. [24] is obtained factoring the joint probability φ kk ≈ (1 − ρ k )ρ k in Eq. ( 3).Details of one-vertex solution can be found elsewhere [14].Finally, the factor k − 1 preceding the first summation in Eq. ( 4) is due to the k neighbors of middle vertex except the link of the pair [0 k 0 k ] (similarly for k − 1 preceding the second summation).
We now approximate the triplets in Eq. ( 4) with a standard pair-approximation [28,29] to find Substituting Eqs. ( 2) in ( 6) and performing a linear stability analysis around the fixed point ρ k ≈ 0 and φ kk ≈ 0, one finds The next step is to perform a quasi-static approximation for t → ∞, in which dρ k /dt ≈ 0 and dφ kk /dt ≈ 0, to find Finally, we plug Eq. ( 8) in Eq. ( 3) to produce where the Jacobian L kk is given by with δ kk being the Kronecker delta symbol.The absorbing state is unstable when the largest eigenvalue of L kk is positive.Therefore, the critical point is obtained when the largest eigenvalue of the Jacobian matrix is null.Let us focus only on uncorrelated networks where Since C kk > 0 is irreducible (all compartments have non-null chance of being connected) and u k > 0, the Perron-Frobenius theorem [18] warranties that Λ is the largest eigenvalue of C kk .The transition point is, therefore, given by −1 + Λ = 0 that results the transcendent equation that can be numerically solved for any kind of network (section 3).
To check the consistency of the theory, we consider the random regular networks (RRNs) that are strictly homogeneous networks with vertex degree distribution P (k) = δ k,m and connections done at random avoiding self and multiple edges [31].Upon substitution of P (k) in Eq. ( 12), one easily shows that the transition point is that is the same of the simple homogeneous pair-approximation.Simulations of CP on RRNs with m = 6 yield a critical point λ c = 1.2155(1) [32], slightly above the pair-approximation prediction λ c = 1.2.

Threshold for arbitrary random networks
In this section, we compare the thresholds given by Eq. ( 12) with simulations of the CP dynamics on random networks generated by the uncorrelated configuration model (UCM) [33].Power law degree distributions P (k) ∼ k −γ , where γ is the degree exponent, with minimum degree k 0 and structural cutoff k c = N 1/2 , the latter rendering networks without degree correlations [30], were used.This choice is suitable for comparison with the pair HMF theory where such a simplification was adopted.We investigated networks with either k 0 = 3 or 6.The latter is to compare with the results of Ref. [27] and to remark the improvement of the present theory.Networks of sizes up to N = 10 7 and degree exponents γ = 2.3, 2.5, 2.7, 3.0 and 3.5 were analyzed.
The thresholds for heterogeneous pair-approximations were determined for each network realization and averages done over 10 independent networks.Sample-to-sample fluctuations of the threshold positions become very small for large networks.The thresholds against network size for two degree exponents are shown in Fig. 1.The results are compared with the heuristic formula inspired in the HPA theory given by Eq. ( 1).We performed simulations of CP dynamics on the same network samples used   17), used to extrapolate the infinitesize limit of the thresholds.Acronyms: PHMF (pair heterogeneous mean-field , HPA (homogeneous pair approximation), HTP (homogeneous triplet approximation).
to evaluate the mean-field theories.The standard simulation scheme was used [22]: An occupied vertex j is randomly chosen.With probability p = 1/(1 + λ) the selected vertex becomes vacant.With complementary probability 1 − p one of the k j nearestneighbors of j is randomly chosen and, if empty, is occupied.The time is incremented by ∆t = 1/[(1 + λ)n(t)], where n(t) is the number of particles at time t.To overcome the difficulties intrinsic to the simulations of systems with absorbing states [22], we used the quasi-stationary (QS) simulation method [34], in which every time the system tries to visit an absorbing state it jumps to an active configuration previously visited during the simulation (a new initial condition).Details of the method with applications to dynamical processes on networks can be found elsewhere [15,32].
The QS probability P (n), defined as the probability that the system has n occupied vertices in the QS regime, was computed after a relaxation t r = 10 6 during an averaging time t a = 10 7 .The transition point for finite networks was determined using the modified susceptibility [31] which is expected to have a diverging peak that converges to the transition point when the network size increases.
The choice of the alternative definition, Eq. ( 14), instead of the standard susceptibility χ = N ( ρ 2 − ρ 2 ) [29] is due to the peculiarities of dynamical processes on complex networks.For example, the CP on annealed networks §, for which the QS probability distribution at the transition point has the analytically known form [15] where Ω = N/g, g = k 2 / k 2 and f (x) is a scaling function independent of the degree distribution.It is easy to show [23] that Using the scaling properties of g [30], for cutoff scaling as k c ∼ N 1/ω , one concludes that, at λ = λ c , χ ∼ N ϑ and χ ∼ N ϑ where ϑ = min[(γ − 3 + ω)/2ω, 1/2] > 0 and ϑ = min[(γ − 3)/ω, 0] ≤ 0. So, the susceptibility χ always diverges at the transition point while χ does not.Typical susceptibility versus λ curves are shown in Fig. 2. The peak positions shift leftwards converging to a finite threshold as network size increases.Notice that the larger the degree exponent the narrower the susceptibility curves and the faster the convergence to the asymptotic threshold.The infinite-size threshold λ * c is estimated in QS simulations as well as in the mean-field theories using an extrapolation As one can see in Fig. 1, the curves λ c vs. N for different mean-field theories are only shifted indicating that the exponents b i are the same.They can be obtained using a continuous approximation in Eq. ( 1) to obtain b 1 = b 2 = (γ − 2)/ω for k c ∼ N 1/ω , where ω = max(2, γ − 1) for UCM networks [33].These b i exponents can also be derived directly from equation (12) in a more complex calculation that is omitted for sake of brevity.§ In annealed networks, the vertex degrees are fixed while the edges are completely rewired between successive dynamics steps implying that dynamical correlations are absent and HMF theory becomes an exact prescription [14].1.
Transition points of the CP on UCM networks with different degree exponents, minimum vertex degree k 0 = 3 or k 0 = 6, and structural cutoff k c = N 1/2 for pair heterogeneous mean-field (PHMF) theory and QS simulations (λ p ). Number in parenthesis are uncertainties in the last digit.
The exponents b i in QS simulations differ from those of the mean-field theories.They can be analytically estimated using the scaling theory presented in Refs.[14,15].The QS density at the transition point scales as where β = max[1, 1/(γ − 2)] [14].These scaling laws are confirmed in the pair HMF theory developed in section 4. Assuming that both scaling laws hold at λ p one obtains Using again the continuous approximation to compute g and neglecting higher order terms one finds where, We performed non-linear regressions using Eq. ( 17) with λ * c , a 1 and a 2 free and fixing b i according to the theoretical corrections.Excellent fits were obtained, as can be seen in Fig. 1 and the numerical estimates of λ * c are shown in table 1.As expected, pair HMF theory is a very good improvement when compared with the one-vertex approximation λ c = 1.However, for some values of γ the heuristic HPA theory is closer to simulations than the pair HMF theory, as can seen in Fig. 1.It is a surprising result since heterogeneity is expected to play an important role in dynamical correlations even for degree distributions without a heavy tail as in the case γ > 3.
The puzzle behind this apparent paradox is that cluster approximations underestimate the real threshold and the convergence is expected only in the limit of large cluster approximations.A homogeneous triplet approximation (HTA) for the CP on unclustered networks yields the threshold [32]: Comparing this approximation with simulations, figure 1, one sees that HTA thresholds are, as expected, higher than the HPA ones but overestimate the simulation thresholds for all investigated networks with k 0 , more evidently for k 0 = 3.This result shows that the homogeneous cluster approximations will converge to a threshold above the correct one and they are, in principle, not applicable to the CP dynamics on heterogeneous networks as previously done [11,27].The proximity between HPA theory and simulations is therefore a coincidence.

Critical exponents
In this section, the critical exponents of the CP in the pair HMF theory are derived and compared with results of QS simulations.

Critical exponents in the pair HMF theory for infinite networks
It is well known that cluster approximations of higher orders improve the critical point estimates but does change the critical exponents in lattice systems [29].As expected, the pair HMF theory for the CP yields the same scaling exponents as the one-vertex approximation [12,15,14], changing only the amplitudes and the finite-size corrections to the scaling as we will show in this section.In a pair level, the scaling exponents associated to the absorbing state phase transition can be derived from Eqs. ( 3) and ( 6) keeping terms up to second order.

Assuming again uncorrelated networks, the dynamical equations become
and The quasi-static approximation with dρ k /dt ≈ 0 and dφ kk /dt ≈ 0 leads to which is inserted in Eq. ( 24) to result and, consequently, the stationary density where Θ i are given by and where ρ = k P (k)ρ k while A i are constants of order 1 given by and The rightmost sides of Eqs. ( 29)-(31) were obtained using Θ i < ρ (the proofs of these bounds are simple) and Eq. ( 28) in a self-consistent approach [35].The constants a i are of order 1/ k 2 and their explicit forms are omitted.
Multiplying Eq. ( 27) by P (k) and summing over k (k c → ∞) one finds where ] and F (a, b, c, x) is the Gauss hypergeometric function [36].Near to the critical point, ψ i 1, we can use the asymptotic form of F (a, b, c, x) to finally find where αi (λ), i = 1, 2, are positive parameters whose details are omitted for sake of conciseness.
The stationary density close to the transition point is given by An expansion around λ = λ c yields where the identity λ c = A 1 (λ c ) comes from Eq. (12).Considering only the leading term in ρ one finds At the transition point λ = λ c , equation (37) becomes which yields ρ ∼ t −δ where δ leading to a relaxation time scaling as with a γ-independent exponent ν = 1.The exponents (β, δ, ν ) obtained in this section are exactly the same of the one-vertex HMF theory [24].

Finite-size scaling critical exponents
Finite-size scaling (FSS) exponents associated to the QS state can obtained using a mapping of the CP dynamics in a one-step process [37] as proposed in Ref. [12].For finite-size systems the condition ρk c 1 is applicable for long times and very close to the transition point.So, we approximate Eq. ( 28) by ρ k λkΘ 1 / k which is inserted in Eq. ( 35) to find where the factor g is given by The first term proportional to ρ in Eq. ( 44) represents an annihilation n → n − 1 whereas the second one a creation event n → n + 1.Following the interpretation of Ref. [12], in a mean-field level Eq. ( 44) represents a one-step process defined by a transition rate W (n, m) from a state with m to another with n particles given by At the critical point, we have the additional simplification λ c = A 1 and the transition rate becomes equal to that of the one-step process associated to the CP dynamics in a one-vertex HMF theory [12], with the factor g = k 2 / k 2 replaced by g, given by Eq. (45).The QS analysis of this critical one-step process with the original g factor was done in Ref. [15], whose results are presented below.
The QS probability distribution P (n) is given by Eq. (15) with Ω = N/g.The QS density ρ and the characteristic time τ , defined as ρ = 1 N n nP (n) and τ = 1/P (1) [34], respectively, scale as ρ ∼ (gN ) −1/2 and τ ∼ (N/g) 1/2 .(47) Nevertheless, the factor g has exactly the same asymptotic scaling properties as the factor g, which are given by Eq. ( 16), and therefore the same FSS exponents of the one-vertex HMF are obtained in pair HMF approximation.The scaling laws ρ ∼ N −ν and τ ∼ N α with ν = max[(5 − γ)/2, 1/2] and α = max[(γ − 1)/4, 1/2] are obtained for UCM networks with a structural cutoff k c ∼ N 1/2 [33].Despite of the same asymptotic scaling, the sub-leading corrections in the new factor g are not negligible as one can see in Fig. 3.Moreover, the finite-size corrections in the critical point position observed for pair HMF theory as well as in QS simulations (Fig. 1) suggest that we must compute the critical quantities at λ p (N ) and not λ * c as done previously done [27].Figure 4 shows double-logarithmic plots for the FSS of the critical QS density and characteristic time following this strategy.For the wide range of degree exponents analyzed, the exponents obtained from power law regressions ρ ∼ (gN ) −Sν and τ ∼ (N/g) Sα are in remarkable agreement with the theoretical prediction S ν = S α = 1/2, as one can verify in table 2. Most importantly, the scaling   laws hold for the entire range of investigated sizes in contrast with the analysis for a fixed λ = λ * p and using the old factor g, for which large deviations of the theoretical scaling laws are observed at small sizes, the more evident for more heterogeneous networks (γ ≤ 2.5) [27].Noticeably, the exponent of the characteristic time for γ = 2.3 is in great agreement with the theory if factor g is used in contrast with a poor accordance observed for a similar degree exponent reported in Ref. [27].It is worth stressing that the almost perfect match is found only if both factor g and corrections in λ p (N ) are used concomitantly.In particular, for the k 0 = 3 case the scaling laws obtained in simulations are not consistent with HMF if this strategy is not used.Thus, we filled a missing gap showing that the critical exponents as well the sub-leading corrections to the FSS are very accurately predicted by the pair HMF theory.

Conclusions
The dynamics of the contact process on the top of complex networks was investigated using a pair heterogeneous mean-field theory in which the vertices are grouped accordingly their degrees.We compared the theoretical results with QS simulations and showed that they represent great improvements in relation to the simple HMF approach.However, for a wide range of the degree distributions, a heuristic homogeneous pairapproximation [11,27] is still more accurate than our heterogeneous approach.To unveil this contradiction we compared simulations with a homogeneous triplet approximation that must be more accurate than homogeneous pair-approximations.We observed, however, that the HTA theory overestimates the simulation thresholds showing that successive homogeneous cluster approximations [28] converge to the wrong critical point and, therefore, that the agreement between HPA and simulations is only a coincidence.We also determined the critical exponents in the pair HMF approach.For the infinite size limit the exponents are the same as the one-vertex theory.However, the finite-size corrections to the scaling obtained in the pair HMF theory allowed a remarkable agreement with QS simulations for all degree exponents (2.3 ≤ γ ≤ 3.5) and network sizes (10 3 ≤ N ≤ 10 7 ) investigated, suppressing a deviation observed for low degree exponents in the one-vertex HMF theory [27].Our results strongly corroborate that HMF theories predict the correct scaling exponents of the CP on SF random networks.
Finally, the present theoretical approaches can be applied to other important dynamical processes on complex networks as generalized voter models [38], sandpiles [39] as well as more sophisticated structures as multiscale and multiplex network [40,41].
Our approach permits to explicitly derive analytical expressions whereas previous pair approximations for dynamical processes in complex networks [42,8] usually needs a numerical integration of the corresponding master equations that limits the analysis to relatively small systems.For example, the threshold of the SIS model in a pair HMF approximation can easily obtained: This threshold coincides with that of the susceptible-infected-recovered (SIR) model in a one-vertex HMF theory [35].This results was recently proposed in Ref. [6] using heuristic arguments.

Figure 3 .
Figure 3. Ratios between the factor g obtained in pair HMF theory, Eq. (45) with λ = λ c (N ), and the factor g = k 2 / k 2 of the one-vertex HMF theory, for k 0 = 3.

Figure 4 .
Figure 4. FSS of the characteristic time and critical QS density for k 0 = 3 (left) and k 0 = 6 (right).The dashed lines have slope ±1/2 as guides to the eyes.

Table 2 .
Critical exponents obtained in QS simulations of the CP on UCM networks with minimum degree k 0 = 3 or k 0 = 6 and cutoff k c = N 1/2 .The exponents were obtained in power law regressions ρ ∼ (gN ) −Sν and τ ∼ (N/g) Sα .