Network Rewiring Models

Recently we showed that a simple model of network rewiring could be solved exactly for any time and any parameter value. We also showed that this model can be recast in terms of several well known models of statistical physics such as Urn model and the Voter model. We also noted that it has been applied to a wide range of problems. Here we consider various generalisations of this model and include some new exact results.


Introduction
Graphs with a constant number of edges and vertices but which evolve by rewiring those edges are a classic network model as exemplified by Watts and Stogatz [1] ( [2,3] provide further examples). Such network evolution may also be recast as other types of statistical physics models (for example see [4,5,6,7]). As many real systems are effectively of constant size, non growing networks can also be used to model a wide range of data: the transmission of cultural artifacts such as pottery designs, dog breed and baby name popularity (as in [8,9,10,11]), the distribution of family names in constant populations, and the diversity of genes. In this paper we look at various extensions to a model of network rewiring for which an exact solution [12] was presented at ECCS06 [13] and in more detail in [14].

The Model
We will study the rewiring of a bipartite graph consisting of E 'individual' vertices connected by one edge only to any one of N 'artifact' vertices, as shown in Fig. 1. At each time step two choices are made. With probability Π R an individual is chosen. It is the artifact end of its single edge, connected to the 'departure' artifact, which is to be rewired. An 'arrival' artifact is also selected with probability Π A . Only after the choices are made is the network altered by rewiring the chosen edge so that its artifact end is moved from the departure to the arrival artifact. Note we do not explicitly exclude the possibility that the departure and arrival artifacts are the same. The individual vertices always retain one edge while the degree k of the artifact vertices is changing in time, only its average degree k = E/N is constant. It is the distribution of the artifact vertices at time t, n(k, t), and its probability distribution p(k, t) = n(k, t)/N , that we study. This process can be viewed in many other ways [14].
The evolution of the degree distribution in the mean field approximation is described by the master equation [12,13,14] Fig. 1. The bipartite graph has E 'individual' vertices, each with one edge. The other end of the edge is connected to one of N 'artifact' vertices. If the degree of an artifact vertex is k then this artifact has been 'chosen' by k distinct individuals. At each time step a single rewiring of the artifact end of one edge occurs. An individual is chosen (number (2i − 1) here) with probability ΠR which gives us the departure artifact (here D). At the same time the arrival artifact is chosen with probability ΠA (here labelled B). After both choices have been made the rewiring is performed (here individual (2i − 1) switches its edge from artifact D to B).
For our physical problem the removal probability must always satisfy Π R (k = 0) = 0 and Π R (k = E) = 1.
In addition for physical solutions we must have n(k, t) = 0 if k < 0 or k > E. The presence of the factors of (1 − Π) ensure that if the degree distribution initially satisfies its physical boundary condition, n(k, t = 0) = 0 if k < 0 or k > E, then this boundary condition is automatically satisfied at all times 4 . The factors of (1 − Π) are not seen in the master equations of the literature [2,3,4,15] and correspond to events where the arrival and departure artifacts are chosen to be the same 5 .
In general the master equation (1) gives the evolution only in the mean-field approximation because we are taking an ensemble average over many instances of the stochastic evolution and using the product of averages where we should have the average of products. However when the normalisations of probabilities Π A and Π R are constant then the master equation (1) may be exact. The most general Π R and Π A for which this is true is We will restrict ourselves to these forms and therefore our analytic results are exact [14]. Thus we are choosing our arrival edge with a mixture of preferential attachment (probability (1 − p r )) and random attachment (probability p r ). The removal artifact is found by choosing the artifact end of a randomly selected edge -'preferential removal'. The use of probabilities proportional to k can emerge naturally through short range searches of many networks, since the probability of arriving at a vertex on a random graph is proportional to its degree [16,17,18]. Not only is the master equation exact for our chosen probabilities (2) but the exact solution for the degree distribution n(k, t) may be found for any finite parameter value. This may be done in terms of (E + 1) eigenfunctions ω (m) (k) and their corresponding generating functions G (m) (x) 4 For this to be true it is absolutely vital that we have the factors of (1 − ΠR(k)) to ensure that with the condition ΠR(k = E) = 1 we do not include processes where an artifact with E edges is lost because we are adding another edge (the third term in (1) for k = E). 5 These events occur with probability (ΠRΠA). Since the network is unchanged by such events, we must exclude such events from the evolution of n(k, t) and the factors of (1 − Π) implement this. It is an approximation to drop these terms. In our model this is not justified for certain parameter values.
The solution is found to consist of simple combinations of the Hypergeometric function F (a, b; c; x) [14] with corresponding eigenvalues, The eigenvalues satisfy λ m > λ m+1 except for p r = 0 when λ 0 = λ 1 = 1.
There is a unique long time equilibrium distribution which can be of one of two phases, as Fig. 2 shows. For p r E −1 we get a condensate, most individuals attach to a single artifact. For p r ≫ E −1 we get a power law degree distribution of unit slope, with an exponential cutoff n(k) ≈ k −1 e −ζk (ζ = − ln (1 − p r )). There is a smooth transition between the two except in the E → ∞ limit.  Using the fact that the number of artifacts N and the number of edges E are constant gives c 0 = N and c 1 = 0 so that the eigenmode numbered one (m = 1) never contributes. Thus the approach to equilibrium of most quantities occurs on a timescale as illustrated in Fig. 3. This means that if we have rewired most of the edges once and almost never used random attachment, i.e. p r E −1 , then the approach to equilibrium is slow, τ 2 = O(E 2 ). However for other cases, p r ≫ E −1 , the small amount of randomness gives a rapid approach to equilibrium after every edge has been rewired just a few times. The initial conditions determine the remaining c m (m > 1).
Of particular interest are the homogeneity measures F n (t) which are the probability that n randomly chosen but distinct edges all share the same artifact. These are given by  Fig. 3. Plots of p(k) from simulations (data points) and the exact analytic results (lines) for E = N = 100 with pr = 10/E on the left and pr = 0.1/E on the right. The results are shown at four different times: t ≈ τ2 (red, crosses), t ≈ 2τ2 (green, circles), t ≈ 3τ2 (blue, stars) and to equilibrium (magenta, squares). The initial configuration has one edge per artifact. The data points are averages over 10 5 runs while the lines are the exact analytic results.
The properties of the Hypergeometric function mean we can express the solutions in terms of fixed fractions of a large number of fixed Gamma functions with all the dependence on time and the initial conditions is carried by factors of c m (λ m ) t . Also the n-th homogeneity measure F n (t) has contributions only from the eigenfunctions m ≤ n.
For instance the most useful homogeneity measure is F 2 (t): while the initial conditions set F 2 (0).

Phase Transitions of Unipartite Graphs in Real Time
The construction of Molloy and Reed [19] gives a unipartite graph of a given degree distribution but is otherwise random. In terms of our bipartite graph this is equivalent to taking pairs of individual vertices and merging the edge ends ('stubs') coming out of these individual vertices. The individual vertices are then thrown away. This is illustrated in Fig. 4. The rewiring of our bipartite model is then equivalent to a rewiring of the projected unipartite graph with the same linear attachment and removal probabilities, also illustrated in Fig. 4. Since the degree distribution of our artifact vertices is also the degree distribution of the unipartite graph, all our results can be applied directly to such graphs. For instance for p r = 1 we capture the degree distribution of the original Watts and Stogatz model 6 [1].
Analytic expressions for the global properties of such random graphs in the infinite N limit depend on a ratio, z, of the second and first moments of the degree distribution [16,17,18,19] There is a phase transition in the properties of such infinite random graphs at z = 1. This occurs when there is one tadpole (an edge connected at both ends to the same vertex) in the unipartite graph. In particular for z > 1 the average distance between two vertices in the giant component, l , may be estimated to be 7 [18] For simplicity we consider graphs where N = E which start with each artifact connected to only one individual so n(k, t = 0) = Eδ k,1 . The projected unipartite graph has k = 1 and initially steps. Only when we start to get a high degree node, so a condensate is forming and p r O(E −1 ), do we get a slower approach to equilibrium on a time scale τ 2 = O(E 2 ). The phase transition in infinite random graphs occurs at z = 1. In our case, our projected graphs start from z(0) = 0 but they reach z = 1 very quickly at . That is even if the evolution to the equilibrium distribution is slow, provided a reasonably large degree node exists, i.e. there is significant amount of copying, a large component emerges in the projected unipartite graph quickly, typically at t 1 ≈ E/2, since The numerical results for the evolution of the properties of the projected unipartite graph are shown in Fig. 5. The parameter z reaches the value 1 at (t 1 /E) ≈ 0.5 ± 0.0002 as expected. This is close to, but not exactly equal to, (t p /E) = 0.535 ± 0.005, where t p is the time at which the average distance and diameter of the largest component peak. The second derivative in time of the number of vertices in the largest component also suddenly switches sign at exactly the same time t p . The value of z at this time is z(t p ) = 1.06 ± 0.01.
Motivated by the approximate expression for the distance in the largest component of a large random graph (13), we find that the inverse distance for the parameters used in Fig. 5 is well fitted by the form a ln(z − 0.06) + b + cz + dz 2 but with different values either side of the peak time 8 . The fit is shown in Fig. 6. The deviations from the predicted z = 1 transition point seem to be finite size effects. The peaks 6 Strictly speaking we choose a random edge to rewire while Watts and Stogatz [1] rewired in a systematic manner. 7 This formula must be adapted from [18] to take account of the existence of vertices of zero degree. Analytic derivations of such global properties use an ensemble of graphs over which there is always a finite probability of getting from any one vertex of degree ki > 0 to a vertex of degree kj > 0 in a finite number of steps. In any one graph this need not be true. On the other hand numerically we measure the average distance in the largest component of one graph considering only vertices in the largest component. We then average this result over the ensemble of graphs. Numerical evidence suggests that this numerical measurement has the same qualitative behaviour as the analytic formula.   Fig. 5. Properties of the undirected random graph formed using a Molloy-Reed type projection [19]. The underlying bipartite graph has N = E = 10 5 starting from F2(0) = 0 and rewired with pure copying (pr = 0.0). Results are calculated for each instance and then averaged over a total of 1000 runs. . are sharper and closer to z = 1 as the network get larger 9 but with the same k , p r and F 2 (0) = 0. The network shown is evolved with pure copying p r = 0 so in this case the equilibrium distribution, a complete condensate F 2 = 1, will emerge only on a long time scale of One way to look at this transition is to use the interpretation of the model in terms of cultural transmission [8,9,10,11,12,13,14]. In this case the bipartite graph represents individuals who are choosing artifacts by either copying the choices made by another individual (preferential attachment) or by making their own innovation (random attachment). Suppose we now imagine that each person has two copies of an artifact. The unipartite graph is then one expression of the relationship between objects as defined by the choices made by individuals. For instance one could imagine asking people to categorise their two favourite pairs of shoes and each artifact could represent a different category, e.g. one artifact might represent black leather lace up shoes. The unipartite projection gives a metric in artifact space as defined by the choices made by the individuals. The phase transition in the unipartite network then marks the point where the individuals have reached some sort of consensus as the artifacts now form a Giant Connected Component given the metric provided by the individuals' choices.

Voter Models and Individual Networks
One possible generalisation of our rewiring model is to add a second graph connecting the individual vertices which we will call the Individual graph. When an individual rewires using preferential attachment they copy the artifact chosen by one of their neighbours in the Individual network. With p r = 0 and N = 2 we obtain the basic Voter model [6,7]. Our model corresponds to having a complete graph for the individual network but with the addition of a random rewiring process, p r > 0, and an arbitrarily large number of choices, N ≥ 2. Neither of these cases is studied in the Voter model literature where the focus is on different types of individual networks and any analytic results are only available for the E → ∞ limit [7]  Results for the equilibrium distribution show that it is qualitatively unchanged by the type of individual graph 11 except for the case of a one dimensional ring, as shown in Fig.7.
We will use two quantities to study the behaviour of the model. Our quantity F 2 of (10) is a measure of the global homogeneity. An equivalent measure which takes account of the local properties of the Individual network is the average interface density, ρ , the probability that any two individual vertices connected by the Individual graph have a different artifact. If the graph is complete or if the Individual graph is ignored (p r = 1) then from (11) we have ρ t = 1 − F 2 (t). Otherwise for two reasons we expect that ρ t = 1−F 2 (t) and that both would both differ from value obtained for a complete Individual graph as derived from (11). First because of the explicit reference to the Individual graph in the definition of ρ t but not in F 2 . Second, the structure imposed by the Individual graph will, in general, effect both the evolution timescale and, for p r > 0, equilibrium degree distributions as compared to the complete Individual graph case.
We can see these differences if we compare the equilibrium values reached on lattices of different dimensions but with some randomness present (otherwise ρ t = (1 − F 2 (t)) because both are zero). As Fig.8 and table 1 show, the local and global homogeneity measures ρ and (1 − F 2 ) are close to the analytic result for large dimension lattices with short network distances. As we take lattices of smaller dimension, F 2 gets much larger than the analytic result, and ρ much smaller. Table 2 shows a similar effect as we increase p r .
The time scale of the approach to consensus is often studied in Voter models. If p r > 0, so there is no absolute consensus, the approach to equilibrium, as measured by F 2 and ρ, is controlled solely by the time scale of the second eigenvalue τ 2 of (9) if the Individuals are connected by a complete network. For general Individual networks the evolution of the global F 2 or local ρ takes the same form as (11), a exp(−t/t 0 ) + c. However, as one might expect, the formation of small patches of consensus between nearest neighbours, measured by ρ , happens faster than the emergence of a global consensus, as measured by F 2 . This is accentuated if there is a large distance between individuals as the comparison between lattices of different dimensions in Fig.8 and in table 1 show. Varying p r also shows that local equilibration is faster than global but there does seem to be a marked difference between p r ≪ 1/E and p r 1/E as table 2 shows. 12 For p r E ≪ 1, local equilibration is a little slower than occurs on complete graph. However for p r E 1, this randomness brings local equilibrium an order of magnitude faster than was the case with a complete graph. It shows that a little bit of randomness can speed up local equilibration but not if an overwhelming consensus is going to emerge. In table 3 we see that the time scales for the exponential decay obtained from fitting our data are roughly in line for the predictions made for the completion time in this model [6,7,20] on a lattice, t 0 ∼ O(E) for a one-dimensional lattice, t 0 ∼ O(ln(E)) in two dimensions. However some discrepancies suggest more work is needed.
The main result to draw from table 3 is that the evolution towards equilibrium, its time scale and final value, are independent of the number of artifacts. This is to be expected at small p r given our linear attachment probabilities as this gives our model certain scaling properties [14]. Suppose we have the consensus emerging picking out one of our N artifacts and we merge the remaining N − 1 artifacts into one artifact. The probability of an individual copying the consensus artifact or one of the remaining artifacts is exactly the same as if we had a model with N = 2 and the same (1 − p r ). The only difference  Table 3. Table of time scales in units of τ2 for E individuals connected by a one-or two-dimensional torus found by fitting the data to a exp(−t/t0) + c. For two and ten types of artifact, pr = 0. Data was averaged over 1000 runs for the largest lattices down to 50 runs for the smallest lattices. Where the error is known reliably, the numbers in brackets specify the error in the last digit.
is that when a random innovation event occurs, with probability p r , the non-consensus artifacts are preferred to the single consensus artifact by a factor of (N − 1). Thus for N ≫ 2 the random events are more likely to destroy the emerging consensus than in the Voter model but only if p r ≫ 0. For the extreme case of p r = 0 we see the expected lack of dependence on N in table 3. The only effect of increasing the number of artifacts in our results comes from starting from a homogeneous initial condition so F 2 = 1/N which is further away from F 2 = 1 and consensus, see Fig. 9.

Two Types of Individual
Another variation of our original model is to introduce two types of individual, labelled X and Y . At each time step we first pick which type of individual to update; with probability q x we select at random one the X-type individuals. We rewire its artifact end, choosing its arrival artifact in one of three ways: at random, by copying the existing choice of one its own type of individual, or finally copying the existing choice made by a random individual of the opposite type. These arrival probabilities may be different for the two types so we have four independent arrival probabilities and one departure probability. Add in the freedom to choose different numbers of X and Y individuals, E x and E y , and N the number of artifacts, we find we have eight free parameters. The degree distribution is now n(k x , k y ; t), the number of artifact vertices which have k x (k y ) edges to X (Y ) type vertices at time t.
The question is can we solve this system analytically? By keeping our probabilities linear in degree and because our normalisations are constants of the evolution, our mean field equation is again exact, for the same reasons as in the original model [14]. Writing in terms of the generating function G(x, y, t) :=

Ex kx=0
Ey ky=0 x kx y ky n(k x , k y , t) we find that we can again split this into (E x + 1)(E y + 1) eigenfunctions which we label with a pair of indices (M, A): where f (MA) ij are constants. The eigenfunctions satisfy a two-dimensional second order PDE. We have not found a full solution but we can reduce this to a one-dimensional problem. Since we express our eigenfunctions in powers of (x − 1) and (y − 1), the eigenfunctions satisfy 13 f

Conclusions
In this paper we have looked at a variety of extensions to the basic network rewiring model of [12,13,14]. Studying the projection onto a unipartite graph gives us exact expressions for the time evolution of a finite sized system through a transition.
We have also shown that adding an Individual network leaves the qualitative behaviour of the model is unchanged in terms of F 2 . However quantitative differences are highlighted by comparison against the case of a complete graph for which our previous analytic work [12,13,14] provides exact analytic results. What we learn from this model is that the consensus (the condensate) may not be perfect, 1 p r E > 0, and it may emerge very slowly τ 2 ∼ O(E 2 ), but an effective consensus is always reached very quickly t 1 ∼ O(E).
Finally we have shown how some progress can be made on solving models with more than one type of edge. In particular we show how the various homogeneity measures F mn may be found exactly.

Additional Figures
These were not included in the proceedings. Data in some of the tables was derived from these curves.