Statistically consistent coarse-grained simulations for critical phenomena in complex networks

# Statistically consistent coarse-grained simulations for critical phenomena in complex networks

## Abstract

We propose a degree–based coarse graining approach that not just accelerates the evaluation of dynamics on complex networks, but also satisfies the consistency conditions for both equilibrium statistical distributions and nonequilibrium dynamical flows. For the Ising model and Susceptible-Infected-Susceptible epidemic model, we introduce these required conditions explicitly and further prove that they are satisfied by our coarse-grained network construction within the annealed network approximation. Finally, we numerically show that the phase transitions and fluctuations on the coarse-grained network are all in good agreements with those on the original one.

###### pacs:
05.50.+q, 89.75.Hc, 05.10.-a

## I Introduction

Complex network has been one of the most active research topics in statistical physics and many other disciplines Albert and Barabási (2002); Dorogovtsev and Mendes (2002); Newman (2003); Boccaletti et al. (2006); Arenas et al. (2008). It describes not only the pattern discovered ubiquitously in real world but also a unified theoretical framework to understand the inherent complexity in nature. However, the investigation of large networks, such as human brain that composes of about neurons and synapses Rabinovich et al. (2006), requires tremendous time-demanding efforts. The phenomenological description may capture certain properties of system, but always neglects microscopic information. A promising alternative is to develop coarse-grained (CG) methods, aiming at significant reducing the degree of freedom while proper preserving the microscopic information of interest.

Several CG schemes have been proposed. Renormalization transformation has been used to simplify self–similar networks, and the reduced networks often preserve some topological properties of the original ones Kim (2004); Song et al. (2005); Goh et al. (2006); Radicchi et al. (2008). Spectral coarse graining technique has been proposed, in which the eigenvalues of Laplace matrix of network are almost unchanged, such that the dynamics of random walk and synchronization are preserved Gfeller and Rios (2007, 2008). Equation–free multiscale computational framework has been developed to accelerate simulation using a coarse time stepper Kevrekidis et al. (2003). This approach has been applied to study the CG dynamics of oscillators network Moon et al. (2006), gene regulatory network Erbana et al. (2006), and adaptive epidemic network Gross and Kevrekidis (2008). However, to our best knowledge, no attempt has been made for developing CG simulation method to study critical phenomena (usually described by stochastic models) on complex networks. The size–dependent and scaling behaviors in these systems are studied so far mainly by such as Monte Carlo (MC) and kinetic MC (KMC) simulations Dorogovtsev et al. (2008); Hinrichsen (2000). Apparently, these microscopic approaches are often too expensive. It is noticed that the CG stochastic models have been proposed to study reaction–diffusion processes on regular lattices Katsoulakis et al. (2003). However, the existing methods are largely inapplicable to critical phenomena on complex networks with diversified heterogeneity. Moreover, the crucial issue concerning the methodology development as to what criteria should be met to make the CG model statistically consistent with the microscopic one is yet to be addressed.

In this paper, we propose the degree-based CG (-CG for short hereafter) be a statistically consistent scheme, within the annealed network approximation (ANA) Dorogovtsev et al. (2008); Boguá and Pastor-Satorras (2003); Caldarelli et al. (2002). It may therefore be an efficient and reliable CG method for evaluating the stationary properties of phase transitions and studying size effect on complex networks. We put forward the conditions of statistical consistency (CSC) on both equilibrium and nonequilibrium properties, exemplified with the Ising model and Susceptible–Infected–Susceptible (SIS) model, respectively. We show that the -CG approach that merges together the nodes of similar degrees does warrant the CSC within ANA. The stochastic Ising spin–flip and epidemic spreading dynamics can therefore be evaluated faithfully and efficiently with the -CG networks. The calculated phase diagrams, fluctuation dynamics, and system size–scaling behaviors are all shown in excellent agreements with the corresponding microscopic MC and KMC results.

The rest of this paper is organized as follows. In Sec.II, we present a general scheme for coarse graining network and -CG approach, and further prove that this approach satisfies statistical consistency. Extensive numerical demonstrations of the CG approach are performed on diverse networks in Sec.III. At last, main conclusions and discussion are addressed in Sec.IV.

## Ii Coarse Graining Procedure

### ii.1 Network Coarse Graining

Let us start with the basic ingredients of network coarse graining. Consider a network consisted of nodes whose connectivity is given by the adjacency matrix , in which if nodes and are connected and otherwise. Merging nodes together into a CG-node (denoted by ) leads to a CG–network with CG-nodes. We adopt the mean–field definition of the CG connectivity between and , as the average number of links connecting any two nodes inside and . The adjacency matrix of the CG network is then

 Acμν=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩2qμ(qμ−1)∑i,j∈Cμ;i

For illustration, Fig.1 depicts an example of coarse graining a network with six nodes into a weighted CG network with three CG-nodes.

### ii.2 Csc

We address the issue on the statistical consistency of a CG scheme with the microscopic network, in terms of both the equilibrium distribution and the nonequilibrium flow. We exploit the Ising model and the SIS model as the paradigmatic examples for equilibrium and nonequilibrium systems, respectively.

For the Ising model defined on a network, the Hamiltonian is given by , with the spin variable , and the ferromagnetic interaction parameter . The probability of a given microscopic spin configuration is given by canonical distribution , where is the inverse temperature and is the partition function. The corresponding CG Hamiltonian assumes the sum of pairwise interactions inside and between the CG-nodes, i.e. ,

 Hc1 =−J2∑μAcμμ[n+μ(n+μ−1)+n−μ(n−μ−1)−2n+μn−μ] =−J2∑μAcμμ(η2μ−qμ), (2) Hc2 =−J∑μ<νAcμν(n+μn+ν+n−μn−ν−n+μn−ν−n−μn+ν) =−J∑μ<νAcμνημην. (3)

Here, is the CG spin variable, and is the number of up/down spins inside . Note that above is a closure expression at the CG level that depends only on and . To make the CG model consistent with the microscopic one, it demands that the probability of any given CG configuration in equilibrium be the sum of the probabilities of all microscopic configurations that contribute to it. That is

 g({ημ})e−βHc=∑{si}∏μδ(ημ−∑i∈Cμsi)e−βH, (4)

where is the number of microscopic configurations corresponding to .

We now consider the nonequilibrium scenario, exemplified by the SIS network for epidemic spreading dynamics. At the microscopic level, the SIS network nodes represent individuals being either susceptible () or infected (). A susceptible individual can get infected at rate , where is the spreading flow from an infected individual to . On the other hand, an infected node can recover to susceptible state at rate . Without loss of generality, we set hereafter to scale the infection parameter . The SIS model exhibits a nonequilibrium dynamical phase transition at from an absorbing state (all recovered) to an active state (disease spreading persistently) Hinrichsen (2000); Keeling and Eames (2005). For the CG–SIS network, we define the CG variables as the number of infected individuals inside a CG-node by . Written in a closure form, the CG recovery rate is , and the CG spreading flow between two CG-nodes is . For this nonequilibrium system, the CG model is consistent with the microscopic one if the CG flow matches the microscopic flows:

 fcμ→ν=∑i∈Cμ,j∈Cνfi→j (5)

Note that the CG recovery rate of holds trivially.

### ii.3 Degree-based CG Scheme

In our approach, we merge the nodes with similar degrees together. We shall show the resulting -CG scheme does satisfy the CSC, as defined by Eq. (4) for the Ising model and Eq. (5) for the SIS model, within the ANA for the ensemble averaged dynamics Dorogovtsev et al. (2008); Boguá and Pastor-Satorras (2003); Caldarelli et al. (2002). In many previous studies Bianconi (2002); Pastor-Satorras and Vespignani (2002), ANA has been extensively confirmed to be a useful tool for describing quenched networks, as in case of the present paper. Although ANA is just an approximation, it still gives a reasonable description of the average behavior of nodes of the same degree. Moreover, one can consider that the ANA is a statistical characterization of large number of quenched networks with the same degree distributions. Nevertheless, in Ref.Lee et al. (2009) it has been pointed out that there exists some discrepancies between considering annealed approximation for quenched networks and considering annealed network models by themselves. According to the ANA, one can replace the dynamics on a given network of nodes by that on a weighted graph of the full connectivity , where and are the degrees of nodes and , respectively, and is the mean degree. In an ideal -CG scheme, the microscopic nodes in a single CG-node are of the same degree: . The CG connectivity is then , for the CG–network dynamics treated at the ANA level. As results, Eqs.(2) and (3) become, respectively (noting that )

 Hc1 =−J2N⟨k⟩∑μK2μ[∑i,j∈Cμsisj−qμ] =−JN⟨k⟩∑μK2μ∑i,j∈Cμ;i

Their sum can be written as

 Hc=−JN⟨k⟩∑μ,νKμKν∑i∈Cμ,j∈Cν;i

It is identical to the microscopic Hamiltonian at the ANA level. In other words, the Hamiltonian of any -CG configuration equals to the collection of its contributing microscopic configurations. The pre-exponential terms in two sides of Eq.(4) are both the degeneracy of CG configuration, as well as the identical energy factors (). We have thus proved that the ideal -CG approach to the Ising model obeys the CSC of Eq.(4) exactly. It is also easy to prove that the nonequilibrium CSC of Eq.(5) is true for the SIS model in consideration; both sides there equal to . Therefore, the -CG approach satisfies the CSC for both equilibrium probability distributions and nonequilibrium dynamical flows. Certainly it is anticipated that the CSC holds approximately if merged together are the nodes of similar degrees rather than the exactly same ones.

## Iii Numerical Demonstrations

### iii.1 CG-MC and CG-KMC simulations

The MC simulation with the Metropolis dynamics Landau and Binder (2000) and the KMC simulation Gillespie (1977) are applied to the Ising model and SIS model, respectively, at both the microscopic and the CG levels. In the CG-MC simulation, a CG-node is randomly chosen, followed by the Metropolis try for spin-flip process. The probabilities of flipping a up/down–spin is , where , and is the change of CG Hamiltonian resulting from the flip of a up/down spin. It is easy to confirm that the CG-MC simulation obeys the detailed balance condition. In the CG-KMC simulation, following the Gillespie algorithm Gillespie (1977), a CG process to be executed is randomly selected based on the transition rates of all processes. The configuration and transition rates are updated for executing the next CG process.

### iii.2 Scale-free networks

We first consider the Barabási–Albert (BA) scale-free network Barabási and Albert (1999), with the degree distribution follows a power-law with scaling exponent . Figure 2 plots typical time evolutions of the magnetization in Ising model at (in unit of ) and the density of infected nodes in SIS model at , where and are used. For both the microscopic and CG simulations, the systems attain the steady states associated with fluctuating noise after transient time. It is clear that there are in good agreement in the steady-state values of and , as well as their fluctuating amplitudes for both simulations cases.

For the Ising model, we construct initial configurations by preparing each node with a random spin value or with an equal probability. As the simulation proceeds, the system quickly relaxes to an equilibrium state. With the same temperature we run at least times of simulations corresponding to different initial configurations and network realizations. In each simulation, MC steps (MCS: each spin is attempted to flip once on average during each MCS) are performed and the last MCS are used to investigate the system’s behavior. As decreases the value of the magnetization undergoes a transition from zero to nonzero at the critical temperature . Below we notice that due to finite-size effects the system can switch between two stable states via a nucleation mechanism, resulting in the oscillations of Aleksiejuk et al. (2002). Above , fluctuates around zero, and the susceptibility per node has a maximum at phase transition, which can be used to determine as we shall show in Fig.3(b). The susceptibility is related to the magnetization fluctuation via the fluctuation–dissipation theorem. To avoid the offset of due to its oscillations in simulations, we will use instead the absolute value of in the following. and are shown in Fig.3(a) and (b), respectively, as a function of the temperature for different and . Apparently, the -CG results (reported by symbols) are in excellent agreements with the microscopic–level counterparts (by solid and dashed lines). As contrast, we also report the results of random-merging (RM) CG scheme (by dotted lines) in Fig.3(a) and (b), where each CG–node includes (, ) nodes selected randomly from the whole network. Evidently the random scheme fails badly in reproducing the microscopic behaviors. We also used the same in the random scheme as in the case of -CG scheme, and found that the two random schemes produce the same results.

It is important to merge nodes with similar degrees together, as already shown both analytically and numerically. Strikingly, even when the original network is reduced to one with only CG-nodes, the -CG scheme still faithfully reproduces the phase transition curves and fluctuations properties. Since is largely reduced compared to , a considerable speed–up of CPU time, about a factor of 40, is realized for . A significantly higher gain can be expected as the network gets larger, allowing the computational study of network size effect very affordable.

Using the -CG approach with a fixed size of CG–networks , we calculate the dependence of phase–transition critical on the network size , as reported in Fig.3(c). It had been well established that the Ising model with ferromagnetic interactions on BA scale-free network undergoes a phase transition from ferromagnetism to paramagnetism at a critical temperature that increases as the logarithm of network size Aleksiejuk et al. (2002); Bianconi (2002); Leone et al. (2002); Dorogovtsev et al. (2002). Our result is consistent with the theoretical expression Bianconi (2002) that , where for the present case of study. Note that the critical phenomenon disappears when .

For the SIS model there is an epidemic threshold on a finite size BA scale-free network network Pastor-Satorras and Vespignani (2001, 2002). Our numerical simulation starts from a random configuration with about half nodes being infected. After an initial transient regime, the system will evolve into a steady state with a constant average density of infected nodes. The steady density of infected nodes is computed by averaging over at least different initial configurations and at least different network realizations with the same parameter . The epidemic threshold occurs at (absorbing state) if and (active state) if Pastor-Satorras and Vespignani (2001). Due to finite size effects, the fluctuation can drive the system to the absorbing state, especially in the vicinity of . Once the absorbing state is arrived, the system will never leave it. Based on the consideration, we use, in practice, a nonzero tolerance in (with the order of ) as the boundary of the phase transition point. Reported in Fig.4(a) and (b) are the calculated results of and its fluctuation , respectively, as a function of , obtained by the CG model and the microscopic model. Again, while the RM scheme fails when the relative infection rate parameter , the agreement between the -CG and the microscopic results remains excellent. Fixing in the -CG approach, the resulting shown in Fig.4(c) is proportional to , also consistent with the theoretical prediction Pastor-Satorras and Vespignani (2002).

Figure 5 show that the simulation results of the microscopic and CG levels on scale-free networks with other two scaling exponents, and . For larger , becomes lower, while gets larger. It is clearly observed that our -CG approach is still applicable. In addition, many other types of networks such as random network and small-world network are also tested, and all results show that the validity of our -CG approach does not depend on network topology.

### iii.3 Degree correlated networks

It is worthy noting that the above numerical demonstrations are carried out on degree uncorrelated networks. We will show that our -CG approach is valid to reproduce critical behaviors on degree–degree correlated networks as well. It has been witnessed that many real networks display different degree–mixing patterns Newman (2002). To measure the degree of the correlation, in Ref.Newman (2002) Newman introduced a degree-mixing coefficient: , where and are the remaining degrees at the two ends of a link and means the average over all links. indicates that there is no degree correlation, while () indicates that a network is assortatively (disassortatively) mixed by degree. Previous studies have revealed that degree–mixing pattern plays an important role in dynamical behaviors on networks, such as percolation Newman (2002), epidemic spreading Boguá et al. (2003), synchronization M. Chavez and Boccaletti (2006). To generate different degree–mixing networks, we employ a algorithm proposed in Xulvi-Brunet and Sokolov (2004). At each elementary step, two links in a given network with four different nodes are randomly selected. To get an assortative network, the links are rewired in such a way that one link connects the two nodes with the smaller degrees and the other connects the two nodes with the larger degrees. Multiple connections are forbidden in this process. Repeat this operation until an assortative network is generated without changing the node degrees of the original network. Similarly, a disassortative network can be produced with the rewiring operation in the mirror method. We start from BA scale-free networks with a neutrally degree-mixing pattern, and produce some groups of degree-mixing networks by performing the above algorithm. Figure 6 displays the results of Ising model and SIS model for three different values of . For each , the simulation results of CG models agree well with those of microscopic ones. In Ising model, shifts to right and at becomes smaller as increases. In SIS model, both and the fluctuation of at decrease with .

## Iv Conclusions and Discussion

In summary, we propose an approach for coarse graining the phase transition dynamics on complex networks described by stochastic models for both equilibrium and nonequilibrium systems. The -CG approaches via degree-based merging scheme are feasible since the reliable microscopic information such as phase transition behaviors and fluctuations are preserved. We have verified that our -CG approach supports, in consistent with the microscopic models, the equilibrium distributions for Ising model and the nonequilibrium dynamical flows for SIS model. Stochastic description, as exemplified here, is ubiquitously important in the study of phase transition dynamics and complex networks for a wide range of realistic systems. Moreover, this work also suggests the development of other promising CG statistical models satisfying CSC.

It is interesting to compare our -CG approach with heterogeneous mean-field theory (HMFT) that successfully predicts and on heterogeneous scale-free networks Dorogovtsev et al. (2002); Pastor-Satorras and Vespignani (2001). Based on the ansatz that nodes with the same degree share the same dynamical properties, HMFT derives a series of coupled mean-field equations for degree-dependent quantities. Recently, Langevin approach together with the HMFT has been developed in Refs.Boguá et al. (2009); Caccioli and Dall’Asta (2009), which have confirmed that such an approach is responsible for the effect of fluctuations in a finite-size network that often play important roles in the vicinity of phase transitions. With regard to the present study, we develop a reliable CG simulation approach, that not only can correctly predict the critical phenomena and fluctuations information, but also can be applied to diverse networks. Especially, it is shown that the validity of our -CG approach on the application of correlated networks; however, in this case the HMFT becomes, in general, very difficult to deal with. This is because that specific formulation of degree–degree correlation are unknown for most of correlated networks, such as , that is the conditional probability of a node of degree being connected to a node of degree . On the other hand, the CSC discussed in this work may provide a solid understanding of the physical mechanism behind the basic assumption of HMFT. Our analysis here may lead to the advancement in efficient and consistent CG approaches for dynamics on surfaces and soft lattices.

Note that the present study is limited to the case of quenched networks, that is, the connectivity of networks is frozen in time. While for the case of annealed networks, i.e. the networks themselves are dynamical objects, our CG approach will encounter some difficulties in application. In this case, since the network connections are frequently reshuffled, the adjacency matrix of the resulting CG–network and CG variables should be accordingly updated. This will lead to the very inefficiency of our CG approach in simulations. However, an important advancement in studying critical phenomena on annealed networks has been made in Refs.Noh and Park (2009); Lee et al. (2009) by means of finite-size scaling theory. Developing coarse-grained simulation methods on annealed networks and adaptive networks Gross and Blasius (2008) deserve further investigations.

###### Acknowledgements.
This work was supported by the National Science Foundation of China under Grants No. 20933006 and No. 20873130.

### References

1. R. Albert and A.-L. Barabási, Rev. Mod. Phys. 74, 47 (2002).
2. S. N. Dorogovtsev and J. F. F. Mendes, Adv. Phys. 51, 1079 (2002).
3. M. E. J. Newman, SIAM Review 45, 167 (2003).
4. S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.-U. Hwang, Phys. Rep. 424, 175 (2006).
5. A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Phys. Rep. 469, 93 (2008).
6. M. I. Rabinovich, P. Varona, A. I. Selverston, and H. D. I. Abarbanel, Rev. Mod. Phys. 78, 1213 (2006).
7. B. J. Kim, Phys. Rev. Lett. 93, 168701 (2004).
8. C. Song, S. Havlin, and H. A. Makse, Nature 433, 392 (2005).
9. K.-I. Goh, G.Salvi, B. Kahng, and D. Kim, Phys. Rev. Lett. 96, 018701 (2006).
10. F. Radicchi, J. J. Ramasco, A. Barrat, and S.Fortunato, Phys. Rev. Lett. 101, 148701 (2008).
11. D. Gfeller and P. D. L. Rios, Phys. Rev. Lett. 99, 038701 (2007).
12. D. Gfeller and P. D. L. Rios, Phys. Rev. Lett. 100, 174104 (2008).
13. I. G. Kevrekidis, C. W. Gear, J. M. Hyman, P. G. Kevrekidis, O. Runborg, and C. Theodoropoulos, Comm. Math. Sci. 1, 715 (2003).
14. S. J. Moon, R. Ghanem, and I. G. Kevrekidis, Phys. Rev. Lett. 96, 144101 (2006).
15. R. Erbana, I. G. Kevrekidis, D. Adalsteinsson, and T. C. Elston, J. Chem. Phys. 124, 084106 (2006).
16. T. Gross and I. G. Kevrekidis, Eur. Phys. Lett. 82, 38004 (2008).
17. S. N. Dorogovtsev, A. V. Goltseve, and J. F. F. Mendes, Rev. Mod. Phys. 80, 1275 (2008).
18. H. Hinrichsen, Adv. Phys. 49, 815 (2000).
19. M. A. Katsoulakis, A. J. Majda, and D. G. Vlachos, J. Comp. Phys. 186, 250 (2003).
20. M. Boguá and R. Pastor-Satorras, Phys. Rev. E 68, 036112 (2003).
21. G. Caldarelli, A. Capocci, P. D. L. Rios, and M. A. M. noz, Phys. Rev. Lett. 89, 258702 (2002).
22. M. J. Keeling and K. T. Eames, J. R. Soc. Interface 2, 295 (2005).
23. G. Bianconi, Phys. Lett. A 303, 166 (2002).
24. R. Pastor-Satorras and A. Vespignani, Phys. Rev. E 65, 035108R (2002).
25. S. H. Lee, M. Ha, H. Jeong, J. D. Noh, and H. Park, Phys. Rev. E 80, 051127 (2009).
26. D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistcal Physics (Cambridge University Press, Cambridge, 2000).
27. D. T. Gillespie, J. Chem. Phys. 81, 2340 (1977).
28. A.-L. Barabási and R. Albert, Science 286, 509 (1999).
29. A. Aleksiejuk, J. A. Holysta, and D. Stauffer, Physica A 310, 260 (2002).
30. M. Leone, A.Vázquez, A. Vespignani, and R. Zecchina, Eur. Phys. J. B 28, 191 (2002).
31. S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Phys. Rev. E 66, 016104 (2002).
32. R. Pastor-Satorras and A. Vespignani, Phys. Rev. Lett. 86, 3200 (2001).
33. M. E. J. Newman, Phys. Rev. Lett. 89, 208701 (2002).
34. M. Boguá, R. Pastor-Satorras, and A. Vespignani, Lect. Notes Phys. 625, 127 (2003).
35. J. M. M. Chavez, D.-U. Hwang and S. Boccaletti, Phys. Rev. E 74, 066107 (2006).
36. R. Xulvi-Brunet and I. M. Sokolov, Phys. Rev. E 70, 066102 (2004).
37. M. Boguá, C. Castellano, and R. Pastor-Satorras, Phys. Rev. E 79, 036110 (2009).
38. F. Caccioli and L. Dall’Asta, J. Stat. Mech. 10, P10004 (2009).
39. J. D. Noh and H. Park, Phys. Rev. E 79, 056115 (2009).
40. T. Gross and B. Blasius, J. R. Soc. Interface 5, 259 (2008).
104166