Generation of scale-free assortative networks via Newman rewiring for simulation of diffusion phenomena

By collecting and expanding several numerical recipes developed in previous work, we implement an object-oriented Python code, based on the networkX library, for the realization of the configuration model and Newman rewiring. The software can be applied to any kind of network and"target"correlations, but it is tested with focus on scale-free networks and assortative correlations. For generating the degree sequence we use the method of"random hubs", which gives networks with minimal fluctuations. For the assortative rewiring we use the simple Vazquez-Weigt matrix as a test in the case of random networks; since it appears to be not effective in the case of scale-free networks we then turn to another recipe which generates matrices with decreasing off-diagonal elements. The rewiring procedure is important also at the theoretical level, in order to test which types of statistically acceptable correlations can actually be realized in concrete networks. From the applicative point of view, its main use is in the construction of correlated networks for the solution of dynamical or diffusion processes by analyzing the evolution of single nodes, i.e., beyond the Heterogeneous Mean Field approximation. As an example, we report on an application to the Bass diffusion model, with calculations of the time $t_{max}$ of the diffusion peak. The same networks can additionally be exported in environments for agent-based simulations like NetLogo.


I. INTRODUCTION
The configuration model is a well-known method for constructing uncorrelated networks having an assigned degree distribution.Such networks can be used, for example, as a connectivity environment for the simulation of dynamical processes or diffusion phenomena.
In particular, one is often interested into scale-free networks, i.e., with degree distribution of the form P (k) = c γ k −γ , where P (k) is the probability that a randomly chosen node has degree k.There are standard commands in software packages like Python's networkX for generating samples from a scale-free distribution, and for performing on the corresponding "stubs" of nodes a wiring procedure uncorrelated from the degree.In Sect.IV A of this work we briefly recall this method and we compare its accuracy end efficiency with another method originally proposed in [1] ("random hubs" method) and here implemented using the object-based functions of networkX.(All codes are available at https://github.com/Ladilu/python-bass-accessible.) Much less is known about the possibility of generating scale-free networks with assigned degree correlations, especially assortative networks, which are usually employed to represent social networks, financial and economic networks, and also some traffic networks and technological networks [2].For this purpose it is possible in principle to apply to an uncorrelated scale-free network the Newman rewiring procedure with assortative target correlations [3].
It is not easy to tell in advance, however, if and when this procedure will be successful.The main purpose of this paper is to explore this issue and find some useful common principles and numerical recipes.
At which level can degree correlations be assigned in a network?Given a statistical degree correlation matrix P (h|k) having the correct properties of positivity, normalization and pseudo-symmetry (and there are many ways of constructing such matrices), when can we say that scale-free networks with those correlations exist, or at least that a statistical ensemble of networks exists, which exhibit those correlations at an average level?Answering such questions is important for the sake of network theory itself, but also because in the so-called Heterogeneous Mean Field (HMF) approximation [4][5][6] several general results concerning diffusion processes have been obtained postulating certain statistical correlations, or just requiring that the function knn (k) derived from the full P (h|k) has certain properties.We recall that knn (k) represents the average degree of the neighbours of a node of degree k: In general, networks with different correlations P (h|k) may have the same knn (k), see [7].
Rigorous results about the convergence of the knn (k) function in random graphs with given joint degree distribution of neighbor nodes and in the configuration model have been given by [8].
A typical example is the assortative matrix by Vazquez-Weigt [9].It has the form and its knn (k) function is kV az nn (k) = (1 − r) The Newman assortativity coefficient r enters explicitly in the definition (in general r spans the range [−1, 1], but here 0 ≤ r ≤ 1).Notice that this knn (k) is linear but does not start from the origin.Using the correlation matrix above one can quickly write diffusion equations in HMF approximation, namely a system of n coupled differential equations, one for each degree class, which e.g. for the Bass diffusion model have the form With p = 0 (no "publicity" or "broadcast" term) these equations reduce to those of the SI ("Susceptible-Infected") epidemic model.However, it turns out that it is not always possible to build real scale-free networks having this kind of correlations.
The outline of the work is the following.In Sect.II we introduce Newman's assortative rewiring by applying it to the case of Erdös-Renyi random networks.In this case the technique works quite well already with the assortative matrix by Vazquez-Weigt.In Sect.
III we show that the rewiring technique fails for Barabasi-Albert networks, partly due to their intrinsic correlations, but also to their scale-free degree distribution, as confirmed in the following, when other scale-free networks are considered.In Sect.IV we look at the general scale-free case with exponents between 2 and 3.For this we first need to implement in an efficient way the configuration model in order to build the uncorrelated networks from which the rewiring process begins.We note that the Chung-Lu method has some limitations under this respect, at least when the size of the networks is not very large.For this reason we introduce the method of the random hubs (Sect.IV B).In Sect.V we apply the rewiring procedure in the general scale-free case and we find that it works quite well with some target assortative matrices which were defined in previous work; they do not have an explicit expression like the Vazquez-Weigt, but the construction algorithm is included in the present code.In Sect.VI we implement the network Bass model in first closure on the single nodes, i.e., with N differential equations coupled via the adjacency matrix (N is the number of nodes).We briefly discuss the outcome concerning the total diffusion curves and the peak time of the total diffusion rate, also comparing the uncorrelated and assortative case.In The general working principle of the Newman rewiring method is the following.First one chooses at random in the list of the links of the network two links (a, b) and (c, d), between nodes a, b, c, d, with excess degrees A, B, C, D. Then one computes the probability E 1 of these links according to a "target" correlation matrix e 0 , namely E 1 = e 0 AB e 0 CD and the analogous probability E 2 for the exchanged links, i.e. for the couple (a, c) and (b, d).After this, a sort of Metropolis-Monte Carlo criterion is applied: if E 2 > E 1 , then the rewiring is performed with probability 1, otherwise it is performed with a probability proportional to the ratio E 1 /E 2 .
In our code the target correlations are defined in terms of the nodes degrees by assigning P (h|k) before the rewiring cycles, and translated into the correlations e 0 ij of the excess degrees with the formula e 0 ij = P (i + 1|j + 1)(j + 1)P (j + 1)/⟨k⟩, written as e0[i,j]=(Phk[i+1,j+1]*(j+1)*probability seq[j+1])/aver degree.The array probability seq gives the normalized degree distribution of the random graph G and is obtained via the function degree dist(G), which in turn uses the command nx.degree histogram(G).Note that the indices i and j are in the range 0, ..., n − 1, while the indices h and k are in the range 1, ..., n.
The knn (k) function of the target correlation matrix is also computed, according to the definition 1, for comparison with the average knn (k) of the rewired ensemble.A counter variable returns for each sub-cycle the number of rewiring steps which have been accepted.
(See typical values in Figs. 1, 2) The choice of the edges for the rewiring is not performed by creating the edge list, but instead the module discrete sequence is called.The latter returns a sample sequence of length n from the discrete cumulative distribution of the degrees.Then two numbers are picked from the sequence; if they are the same, the code skips.The choice among the neighbours for the rewiring is made by means of the built-in python module random, which allows to generate and choose random numbers in a list.The vertices are selected so that they are different from each other.Before finally performing the rewiring, after the condition is proven, another control is made, in order to check whether the edges that are going to be added already exist; in fact, they must not be already present in the graph, or the networkX package will not perform the rearrangement as prescribed.
At the end of each rewiring sub-cycle the r coefficient of the network is computed and stored for a final evaluation of the average and standard deviation of r over all sub-cycles.
As we shall see, this gives only a rough "integral" check of the convergence of the degree correlations to the target correlations.

Moreover, at the end of each sub-cycle the knn (k) function of the current rewired network
G2 is computed with the command knn=nx.averagedegree connectivity(G2,...) and accumulated into a dictionary knn1= dict(sorted(knn.items()))for subsequent evaluation of its final average.All the graphs of knn (k) for each sub-cycle are plotted together in a "cloud" graph which gives a visual impression of their fluctuations.Finally, the ensemble average of knn (k) and the target knn (k) are plotted for comparison.The final size of the giant connected component is also computed.

FIG. 1: The knn functions of an ensemble of 200 networks obtained by Newman rewiring
with target correlations of the Vazquez-Weigt type (r = 0.3), starting from a random network of 4000 nodes, probability connection 0.001.The Newman coefficient of the ensemble is r = 0.20 ± 0.01.The giant component after the last rewiring is 97.9 %.In each of the 200 rewiring sub-cycles, the number of accepted rewiring steps was about 218'000.

III. NEWMAN REWIRING FOR BARABASI-ALBERT NETWORKS
A rewiring procedure similar to the one starting from random networks can be applied to Barabasi-Albert (BA) networks generated by preferential attachment to α existing nodes with the method G=nx.barabasi albert graph(N,alpha).In this case, however, the rewired network turns out to be not assortative, and so the procedure fails.See Figs. 3, 4.
This is due to a limitation of the Newman rewiring when applied to scale-free networks (see Sect.V) and, in addition, to the fact that BA networks have intrinsic degree correlations which cannot apparently adapt well to the Vazquez-Weigt correlation matrix.In fact, although the r coefficient of BA networks is close to zero, their knn (k) function is decreasing at small degrees and increasing at large degrees.
We notice in this regard that the Newman rewiring works quite well, on the contrary, when one starts from an uncorrelated scale-free network with γ = 3 and takes the BA correlation matrix as target [1].
FIG. 2: Average of the knn functions of the ensemble of random networks shown in Fig. 1 (dotted plot).The continuous line shows the target knn (k) function (eq.( 3), with r = 0.3).
The matching between target correlations and correlations of the rewired network is usually assessed by comparing the respective knn functions, at least for degrees which are not too large (the structural disassortativity prevents any real knn to increase indefinitely like a theoretical target assortative knn , because there are simply not enough hubs in a scale-free network for connecting hubs preferably with hubs).A common phenomenon occurring in scale-free networks, however, is the strong fragmentation of the rewired network, especially if the smallest degree is k min = 1.In that case, the nodes with k = 1 represent a vast majority of all nodes, and if P (1|1) ̸ = 0 in the target correlations, then a large number of isolated couples are formed; if P (2|1) ̸ = 0 a large number of isolated triples, and so on, leaving a small principal connected component.This phenomenon can be observed using the method G2.subgraph(max(nx.connectedcomponents(G2), key=len)).The situation improves markedly if the minimum degree is equal to 2. In this case, the principal component is really a "giant" component, close to 100% of all nodes.

IV. CONFIGURATION MODEL FOR SCALE-FREE NETWORKS
We are now going to illustrate a method for constructing uncorrelated scale-free networks with given exponent 2 < γ ≤ 3 and for applying Newman rewiring to them in order to approach some target assortative correlations as well as possible.It is known that many relevant real networks are assortative, in particular social and economic networks; it is therefore very useful, in order to simulate dynamical processes on this kind of networks, to generate them through some algorithm.This issue is also interesting in itself, for "foundational" purposes, in order to understand if the scale-free property is fully compatible with assortative (or possibly disassortative) correlations.
The established paradigm of scale-free complex networks has been recently re-discussed [10,11] and it turns out that in many cases it is impossible, due to statistical errors and fit uncertainties, to state whether certain networks are scale-free or not.Since the preferential attachment method leads to an exponent γ = 3, for the generation of scale-free networks with a different exponent it is usually necessary to employ the configuration model.In Sect.
V B we shall present a novel alternative method [12] in which the excess degree correlations e jk are assigned at the beginning and a good approximation of scale-free degree distributions is obtained as a consequence.3), with r = 0.3).

A. Test of the implementation of the Chung and Lu model
Before describing our implementation of the configuration model, we will discuss the limitations of a method already available in networkX and based on the model by Chung and Lu [13,14].This is a weaker version of the configuration model, treated also by Newman in his book [15] in comparison to the standard configuration model with "fixed degree distribution", in which it is possible to prove analytically several general properties.
Suppose we want to build a network with N nodes, starting from a list of N values for the degrees of the nodes, randomly extracted from a power law distribution with given exponent γ.It is possible to generate such a list with various stochastic methods, in a similar way as for the generation of samples of the normal distribution or of other known distributions.A general procedure makes use of a probability transformation method.One defines first a vector F k = k j=1 P (j) where k = 1, ..., n and P (j) denotes the normalized degree distribution (the power law in our case, with n the maximum degree).The values of F k define breakpoints of the unit interval (0, 1).After generating a random number ξ in this interval, a new node is introduced into the list with degree k if F k−1 < ξ < F k , and the procedure is repeated N times.
In networkX there is an auxiliary function which does just this, namely S=nx.utils.powerlawsequence(N, gamma).The list of values obtained is real, in the interval (1, +∞).For transforming this into a list of degrees, a proper rounding to integers is necessary.It can be checked that for large N the degree succession obtained respects well the scale-free criterion of the ratio of probabilities, namely P (k 1 )/P (k 2 ) = (k 2 /k 1 ) γ .
In the Chung and Lu method (networkX command: G = nx.expecteddegree graph(S, selfloops=False)) one actually starts from the not-rounded scale-free degree succession and connects nodes i and j having degrees k i and k j (i, j = 1, ..., N ) with probability proportional to k i k j .It is possible to prove that the network obtained must be scale-free in the limit N → ∞.The algorithm by Miller and Hagberg employed by networkX [14] runs in an expected time of order N .It is clear that the exact degrees of the nodes cannot be assigned a priori like in the standard configuration model, but become random variables themselves.
Our numerical tests of this method (code at https://github.com/Ladilu/python-bass-accessiblgive, at least for values N up to 10 5 , distributions of the final degrees which deviate significantly from the scale-free criterion (while the initial degree list, for the same values of N , is accurately scale-free).

B. Random hubs method
Therefore we chose to implement the standard configuration model in a custom function named configuration model, in which a multigraph is constructed starting from the assigned degree sequence.In this work we only made use of the code for undirected network.Once the degree sequence is given, a list containing all nodes is built, each node being repeated as many times as its degree determines.Then the list is reshuffled by means of the random module functions and split in half.The vertices for the edges are picked in an ordered way from the two halves of the list, so that the method cannot go back and accidentally select one node more times than allowed.This is possible with the built-in methods of the Python itertools package, which are programmed for operating on iterable objects such as arrays and lists.The edges are finally added to the graph until the list is finished, in such a way that the original degrees of the nodes are preserved.This function works with any assigned degree sequence.For example, it is also possible to get from networkX a random graph or a BA network like in the rewiring examples previously discussed, then extract their degree sequence and apply the configuration model function to it, obtaining an uncorrelated network with the same degree sequence (of course, the random graph should already be uncorrelated).In order to build a scale-free network one could apply the function configuration model to a degree sequence generated with the command S=nx.utils.powerlawsequence(N, gamma) as described above, but we chose instead to define a scale-free degree sequence with "minimal fluctuations" through a technique that we call the "random hubs method".This allows to build networks of relatively small size which are as close as possible to an ideal scale-free distribution, and works as follows.
If the average number of nodes with degree k is smaller than 1, i.e., N P (k) = X < 1, then a node with this degree will be created with probability X. Extending the procedure to all degrees, a random variable ξ ∈ (0, 1) is generated for each value of k, and then shows how many nodes of each degree, starting from k min , must be present in the network in order to satisfy the power law -except for the minimal fluctuations mentioned above.
For the highest degrees, most values will be zero, because only a few hubs will actually be present.This is what one also observes in the degree sequence of real scale-free networks, e.g.BA networks: some hubs, with random degrees, are present, and the remaining high degrees are missing.
The maximum degree n considered for the list s1 is connected to the number of nodes N through the Dorogotsev-Mendes relation.After the list has been built, it is necessary to check that the total number is even, otherwise the random wiring procedure cannot connect all nodes respecting their degrees.From the list s1 another list is built with the command list of nodes = to stublist kmin(s1, kmin), giving the degree sequence of the nodes.
For example, if we are generating a network with exponent γ = 2 (this extreme value of γ is normally excluded, and taken here only for simplicity) and there are 100 nodes with degree 1, the nodes of degree 2 must be 25, or let's say 24, 25 or 26, depending on the integer rounding and on the fluctuations due to the term ξ.The degree sequence list of nodes will start with 100 elements equal to 1, then 25 elements equal to 2, and so on.This degree sequence is fed to the wiring function configuration model, which creates the "stubs" of the network and adds links at random between them, until each node reaches the degree defined by the degree sequence.

V. RESULTS OF THE REWIRING PROCEDURE FOR SCALE-FREE NETWORKS
If we use as target the Vazquez-Weigt correlations, the average r coefficient of the resulting network ensemble is very close to zero, independently from the r parameter of the target.
The ensemble average of knn is linearly increasing only for small values of k, similarly to what happens for the rewiring of BA networks.We can thus conclude that the Vazquez-Weigt assortative correlations are incompatible with a scale-free degree distribution, not only due to structural disassortativity at large degrees, but for all degrees except in a small range close to k min .
One might wonder whether it is possible at all to obtain positive ensemble values of r by rewiring a scale-free network.The answer is affirmative, as shown by Xulvi and Brunet [16] with their "empirical" rewiring method (links pairs are exchanged if the differences between their degrees increase) and also in [1,17] with a maximally assortative rewiring based on a formula for the variation ∆r.
It is also possible to obtain assortative scale-free networks by assortative rewiring with correlations different from the Vazquez-Weigt recipe.This has the advantage, compared to empirical rewiring or maximally assortative rewiring, of maintaining some control about the functional form of the correlations.We shall illustrate two methods, denoted for brevity as BM1 and BM2.

A. BM1 assortative matrices
This set of matrices is built through a procedure [18] in which the matrix elements of P (h|k) are chosen to be largest on the main diagonal and decreasing elsewhere, then properly adjusted in order to satisfy the Network Closure Condition and normalized columnby-column.By performing the Newman rewiring with this kind of matrices as target, one obtains ensembles with r ≃ 0.2 − 0.3 and average knn as shown in Figs. 5, 6.

B. BM2 assortative matrices
An alternative procedure for constructing assortative matrices [12] is to define them with an explicit formula in terms of the excess degrees, like e.g.
where c δ is a normalization constant.From this one can obtain the degree distribution P (k) in the form of a series and check that it behaves with good approximation as a scale-free distribution.The P (h|k) matrix is obtained from e jk and P (k) as usual.The target knn function is, unlike in the previous cases, increasing but definitely non-linear (see Fig. 7).
The disadvantage of this method is that the scale-free exponent γ cannot be fixed a priory but depends on the δ exponent in the definition of e jk .

VI. NUMERICAL SOLUTIONS OF THE BASS DIFFUSION EQUATION ON ASSORTATIVE REWIRED NETWORKS
As mentioned in the Introduction, when we have a concrete realization (actually, a statistical ensemble of realizations) of a network having assigned degree distribution and correlations as close as possible to a certain theoretical target, we can solve systems of differential equations which describe a certain dynamical process unfolding on the network.The equations in the system are as many as the nodes.This means that the behavior of single agents is described much more accurately than in the HMF approximation, at the cost of course of a greater computational complexity which limits the network size to a few thousands of nodes.
The results of the numerical solutions still need to be aggregated in order to examine them.Nevertheless, one observes some interesting differences in comparison to the HMF results.To some extent these differences are due, in the scale-free case, to the real random occurrence of large hubs in the network, while in the HMF approximation all hubs contribute "virtually" to dynamics, each one with a small probability.
We have focused our attention on the network Bass model, which written on the single nodes in first closure [15,19] takes the form where A is the adjacency matrix of the network and the variable X l has to be understood as the expectation ⟨ξ l ⟩ of the non-adoption (ξ l = 0) or adoption (ξ l = 1) state of node l over many stochastic evolutions of the system.The code for the solution is available at https://github.com/Ladilu/python-bass-accessibleand makes use of the Python method odeint which can be imported from the package scipy.integrate.The adjacency matrix is handled by networkX as a sparse matrix, with the commands m = nx.adjacencymatrix(G), A=m.todense().Like in the HMF equations, the imitation coefficient q is normalized to the average connectivity ⟨k⟩, in such a way to allow comparisons between networks of different kind while re-scaling the dominant effect of the average degree on diffusion times.
The Bass model is reduced to a standard SI epidemic model when the publicity coefficient p is set to zero; however, the p-term allows to produce a meaningful diffusion dynamics also starting from null initial conditions, i.e., without initial adopting nodes (the equivalent of infected nodes in the SI model).One can thus define a characteristic peak adoption time t max as the time at which the adoption rate f (t) is maximum.f (t) is the derivative of F (t), which represents the fraction of total adopters as a function of time and has a typical S-shaped plot (F (t) = N l=1 X l (t)).It is also possible to define other characteristic times, like the takeoff time, see [12].
In general one observes that for scale-free networks which have the same critical exponent γ and have been rewired using the same target correlations, t max depends strongly on the degree k max of the largest hub present in the network (at least for the networks with 1000 nodes examined).The presence of a super-hub makes diffusion visibly faster, as could be intuitively expected, and it also appears to make the assortative rewiring process less efficient,    For comparable k max , the peak times are larger than with assortative networks.In the rewiring code, a Vazquez-Weigt matrix with r = 0 is chosen as target.The resulting effective r is slightly negative due to structural disassortativity.
code is also available at https://github.com/Ladilu/python-bass-accessibleand it consists of a few procedures.NetLogo allows the uploading of networks created by means of applying the Python rewiring algorithm with the criterium of setting a target BM1 matrix to an initial graph G.The graph G is built costructing the degree sequence from a power law and then applying the configuration model procedure.The Python code is able to produce a graphml output for the nodes and links of the network .There exists a corresponding method in the NetLogo code , wire graphml, which loads the graphml file and then the procedure create-network uses the method to effectively reconstruct the corresponding network.Another method in the procedure create-network is setup-nodes, which places the network in the screen within the area of a circle so that the coordinates of the points corresponding to the agents are well readable.There exists the possibility of setting an arbitrary number of seed adopters.In the procedure adopt the Bass model is implemented.
The p parameter for the broadcast influence is compared to a random variable, so that according to its value the agents which are subject to publicity adoptions are given the red color.The q parameter for a given node is set using the social influence parameter normalized to the number of neighbors divided by the mean degree of the network, which is computed separately.The adopt procedure is then used in the go procedure that can be set in the interface to make the diffusion process start in the simulation.The time-based approach is chosen, by means of resetting the ticks count before any new start of the Bass model simulation.The process can be monitored on the screen as represented in Fig. 8.A detailed description of the simulations and their results will be presented in a forthcoming paper.
FIG. 8: Some NetLogo views of the agent dynamics for Bass adoption in an assortative scale-free network.The agents that have not adopted are initially colored in blue, then "broadcast" (publicity) adoptions are represented by red color and "network" (imitation) adoptions by yellow color.The network has 500 nodes, γ = 2.5 and is obtained via Newman rewiring with a BM1 target matrix.

FIG. 3 :
FIG. 3: The knn functions of an ensemble of 200 networks obtained by Newman rewiring with target correlations of the Vazquez-Weigt type (r = 0.3), starting from a BA-2 network of 5000 nodes.The Newman coefficient of the ensemble is r = −0.07± 0.06.The giant component after the last rewiring is 100 %.In each of the 200 rewiring sub-cycles, the number of accepted rewiring steps was about 192'000.

FIG. 7 :
FIG.5:The knn functions of an ensemble of 400 networks obtained by Newman rewiring with target correlations of the BM1 type, starting from an uncorrelated scale-free network of 10000 nodes, γ = 2.5.The Newman coefficient of the ensemble is r = 0.337 ± 0.002.The giant component after the last rewiring is 100 %.In each of the 400 rewiring sub-cycles, the number of accepted rewiring steps was about 145'000.

TABLE I :
Total adoption rates as a function of time for a network Bass model running on uncorrelated and assortative scale-free networks with 1000 nodes (system of coupled differential equations for the single nodes).

TABLE II :
[20]ples of peak times t max observed in assortative networks with correlations of the BM1 type, γ = 2.5, N = 1000, k min = 2.Note the strong dependence on the degree k max of the largest hub, as discussed in the text.leading to smaller values of the final r in correspondence of the same target correlations.If we want to assess the effect of assortative correlations we thus need to compare network realizations which have similar k max .Then it turns out that in assortative networks t max is systematically smaller than in uncorrelated networks; in addition, the diffusion rate reaches its peak earlier, and decreases faster after the peak.See Tabs.I, II, III.This outcome differs from the predictions of HMF theory, according to which uncorrelated networks give in general a smaller t max than assortative networks[20].In order to reconstruct dinamically the evolution of the Bass diffusion on an assortative network, preliminary agent-based simulations have been carried out with NetLogo.The

TABLE III :
Same as in Tab.II, but with uncorrelated networks.