Correlations and analytical approaches to co-evolving voter models

The difficulty in formulating analytical treatments in co-evolving networks is studied in light of the Vazquez–Eguíluz–San Miguel voter model (VM) and a modified VM (MVM) that introduces a random mutation of the opinion as a noise in the VM. The density of active links, which are links that connect the nodes of opposite opinions, is shown to be highly sensitive to both the degree k of a node and the active links n among the neighbors of a node. We test the validity in the formalism of analytical approaches and show explicitly that the assumptions behind the commonly used homogeneous pair approximation scheme in formulating a mean-field theory are the source of the theory's failure due to the strong correlations between k, n and n2. An improved approach that incorporates spatial correlation to the nearest-neighbors explicitly and a random approximation for the next-nearest neighbors is formulated for the VM and the MVM, and it gives better agreement with the simulation results. We introduce an empirical approach that quantifies the correlations more accurately and gives results in good agreement with the simulation results. The work clarifies why simply mean-field theory fails and sheds light on how to analyze the correlations in the dynamic equations that are often generated in co-evolving processes.


Introduction
Although much progress has been made in the field of complex networks by studying the structural properties of networks and dynamic processes in static networks, scientists have realized that the structure is often coupled and co-evolving with the dynamic processes via a feedback mechanism. Studies in such co-evolving systems have taken center stage in complex system science in recent years. Many models on how the interactions among individuals in a networking environment would affect the network structure have been constructed and studied [1][2][3][4][5][6][7][8][9][10][11]. Instead of adding yet another model to the list, this work focuses on one particular adaptive voter model (VM) [1], with the aim of analyzing it thoroughly and studying how successful analytical treatments could be formulated. VM of opinion formation typically considered a population with agent's relationship defined by a static network [12][13][14][15] or a dynamic network [1,[16][17][18][19]. Co-evolving networks are usually driven by some adaptive mechanisms. In VM, an agent may switch his opinion depending on the current network structure through the connected neighbors' opinions [1,2,20], and in some models depending even on a longer range effect [21]. This opinion switching process changes the fraction of agents taking on a particular opinion, but keeps the network structure intact. An adaptive mechanism that allows an agent to cut a link to a neighbor of the opposite opinion and rewire it to another agent of the same opinion will change the network structure. Although the cut-andrewire process does not change the agents' opinions, the change in the network structure will affect the neighborhoods of the agents and, in turn, the dynamics of future opinion switching. This constitutes a co-evolving system, in which the agents' opinions and network structure that reflects the agents' relationship are coupled and co-evolving in time. Co-evolving systems based on the dynamic relationship among agents can also be set up via dynamic processes other than the voting mechanism. These models include the adaptive Prisoner's Dilemma [4,5], the dissatisfied adaptive snowdrift game (DASG) [6,7] and the co-evolving epidemics [8][9][10]22].
Results of the numerical simulations showed that the co-evolving opinion switching and rewiring mechanisms strongly affect the agents' opinions and the network structure [1,2,19,23]. Despite the complexity of the co-evolving process, attempts on understanding the results analytically have been made with various degrees of success. A common approach is to set up dynamic equations [2,6,8] to capture the effects of two coupled dynamic processes. The number of equations depends on the extent of spatial correlations to be included. The equations are difficult to solve analytically in general. Thus the equations should be closed by some closure scheme, resulting in a mean-field approximation [20]. Such a scheme often amounts to imposing a cutoff on the spatial correlation. The resulting equations can be studied by fixedpoint analysis [7,8], iterating them in time [2,22] or by using generating functions [19]. Meanfield approximations usually involve the consideration of the neighborhood of an agent on how different opinions, strategies or health status are distributed among the connected neighbors, depending on the context of the process. When the agents evolve in a network with no or weak degree correlations, mean-field theories by invoking a binomial closure scheme give reasonably good agreements with the simulation results [15,24]. Depending on the adaptive mechanisms in co-evolving dynamics, however, spatial correlations are generated and they pose a challenge to the formulation of a reliable theory. In the VM [1], the assumption of a binomial distribution of opinions in an agent's neighborhood gives the correct prediction of a connected network phase with mixed opinions and a disconnected network phase with opinion segregation, but the critical value of the rewiring probability p c at which the transition between the two phases occurs is way off when compared with the simulation results. Recently, the spatial correlations were studied by considering the motifs near the transition [20,25] and a better value of p c was obtained by using a Poissonian degree distribution, and the treatment has also been applied to directed networks [26] and multistate systems [27]. In the co-evolving DASG [22], the adaptive mechanism was found to generate weak spatial correlations so that the binomial closure scheme remains valid. In adaptive epidemic models, a closure scheme that involves pair approximation and a Possion-type closure has been studied [8,28]. In such models, one often needs to consider three-node correlations that describe an agent being connected to two neighbors that could be an infective (I ) or a susceptible (S). A common approximation is to decouple the number of three-node configurations L x yz (where x, y, z ∈ [I, S]) for an agent of a given type y into a product of two-node correlations, i.e. L x yz = L x y L yz /n y where n y is the number of agents of type y, and L x y and L yz are the numbers of links of types x y and yz, respectively. Better theories can be constructed by retaining variables of longer spatial correlations. Marceau et al [29] introduced an improved compartmental formalism for the adaptive epidemic model, in which the variables correspond to the different neighborhoods that an agent may encounter. This amounts to a larger number of variables and thus a larger set of dynamic equations. The method gives very accurate results for some co-evolving systems of two states [30]. In a recent study on a simplified VM [31], however, Durrett et al [32] showed that the approach of [30] can only give results in qualitative agreement with the simulation data. It is the current stage of a great number of co-evolving network models and different analytical approaches of various complexities and degrees of success that motivated the present work. Our idea is to analyze one typical model in depth, study the reason behind the deficiency of a current mean-field approach systematically, and seek ways to make improvements.
In this work, the VM proposed by Vazquez et al [1] is studied in detail. In section 2, we define the VM and identify the key issues. The VM may be numerically unstable for the calculation of the density of active links in finite systems, where the absorption of static state (or consensus state) induced by fluctuations can always be reached. Thus, one needs to be careful in comparing the simulation results and theories, and reliable comparisons are difficult to make. This problem is overcome by introducing a modified VM (MVM) in which an agent could change his opinion with a tiny probability q without any external influence. The q → 0 limit recovers the VM. The MVM is stable and provides a platform for detailed comparisons between analytical and numerical results. In real situations, the mutation of one's opinion may be induced by the public information other than the influence of the opinions of one's neighbors. In section 3, the formalism toward a mean-field theory is examined. It is pointed out that capturing the correlations in two averaged quantities, each involving two types of averaging process, is key to the success of a theory. The assumptions behind the binomial closure scheme are critically examined. In section 4, the failure of the binomial closure scheme is explicitly shown by comparing the simulation results with what is assumed in the closure scheme. The failure asks for taking an agent's neighborhood into consideration. In section 5, we formulate an improved mean-field theory (IMFT) incorporating explicitly the spatial correlation up to the nearest neighbors of an agent of a particular opinion for the VM and the MVM based on the set of variables N A k,l and N B k,l , which correspond to the numbers of agents having opinion A and B with exactly k neighbors and l of them are B-opinion neighbors, respectively. In section 6, an empirical approach of obtaining the two averaged quantities accurately is introduced. The approach gives results in good agreement with the simulation results and a better estimate of p c than the binomial approximation and the IMFT. We summarize the results in section 7.

Models and key problems
We study the co-evolving VM proposed by Vazquez et al [1] and its variation. Consider a system of N agents. Each agent has an evolving opinion that can be A or B. The relationships among the agents are characterized by a total of L tot links, corresponding to a mean degree of k = 2L tot /N for each agent. Initially, each agent is assigned an opinion α (= A or B) and k neighbors randomly, i.e. the initial network is a random network of uniform degree k . The behavior of the system, e.g. whether a consensus can be reached, depends on how the agents change their opinion and their neighborhood. In the VM, the time evolution follows the following dynamics. In a time step, an agent i is picked randomly as a target agent and let α be his current opinion. The target agent selects randomly a neighbor j as the reference agent. If agent j has the same opinion as agent i, then agent i will keep his opinion and nothing will happen. If agent j is of the opposite opinionᾱ, then agent i will have a probability p to cut the link to agent j and rewire to a new neighbor of the same opinion chosen randomly among the non-neighbors in the system. With a probability 1 − p, agent i adopts the opinion of neighbor j and switches his opinion tō α. The active links, therefore, are those that connect two agents of opposite opinions. Both the cut-and-rewire action, which changes agent i's neighborhood, and the opinion switching action, which changes agent i's opinion, serve to make a pair of neighboring agents take on the same opinion. However, the former process tends to separate opposite opinions into disjoint groups and the latter process represents a negotiation between opposite opinions within a connected group. It was reported in [1] that the system would evolve into a connected and mixed phase with agents of opposite opinions coexisting in a connected group, or a disconnected and segregated phase with separated consensus groups, depending on the value of p. These phases and the transition can be characterized by the active link density ρ, which is the ratio of the number of active links and L tot . In figure 1, we reproduced the results of ρ( p) from the simulations for k = 4 together with the mean field There exists a value p c ≈ 0.38 that separates the ρ = 0 mixed phase and the ρ = 0 segregated phase. This work is motivated by several urging questions related to the VM and its analytical treatment. (i) Detailed simulations actually revealed that ρ → 0 in a system with finite N for any value of p, as long as the simulation time is sufficiently long [1]. This indicates that the VM The solid lines are the results based on a mean-field theory using the binomial closure. The inadequacy of the mean-field approach asks for better analytical treatments. The results are obtained for a network with mean degree k = 4.
is not stable. A question is how the model can be stabilized. (ii) The mean-field theory in [1] at best can only capture the trend of ρ( p). The key questions are why the mean-field theory fails and how a better analytical approach can be formulated.
The main reason for ρ → 0 at long time for a finite system in the range of small p is due to fluctuations, which drive the system into a consensus. To stabilize the model, we modify the VM by introducing a mutation probability q that a target agent i switches its opinion without referring to the opinion of the other agents. With probability 1 − q, the target agent i follows the actions in the original VM, i.e. the target agent carries out the cut-and-rewire process with the probability (1 − q) p and the opinion switching process with the probability (1 − q)(1 − p) when the selected referencing agent is of opposite opinion. This MVM gives stable results for a long time even for very small q. Figure 1(b) shows the simulation results for a system with N = 10 000 agents and q = 0.001. The active link drops monotonically with p to nearly zero. This resembles the behavior of ρ( p) in the VM, only that ρ remains non-vanishing for finite q. The small-p phase is highly active with a high active link density and the largep phase corresponds to an almost inactive phase in which agents of the same opinion form groups that are separated from each other. The stability of the MVM allows us to investigate the origin of the invalidity of the analytical approaches. The q → 0 limit of the MVM recovers the original VM.

Formulating the framework of the mean-field theory
Within the VM and the MVM, changes occur when the target and the reference agents have different opinions. Therefore, a starting point for analytical treatment is to establish a dynamic equation for the active link density ρ to describe the evolutionary processes and equilibrium states. Consider a target agent chosen at time t of opinion α with k neighbors, among them there are n k neighbors of the opposite opinion α. The subscript k in n k serves to label the condition that there are k neighbors. With the mutation probability q, the target agent switches his opinion to α, leading to a net change of (k − 2n k ) in the number of active links. With probability (1 − q)(1 − p)n k /k, the target agent chooses a reference agent of opposite opinion and switches opinion to α, leading again to a net change of (k − 2n k ) in the number of active links. The factor n k /k is the active link density of the target agent. With probability (1 − q) pn k /k, the target agent cuts the active link with a neighbor of opinion α and rewires the link to an agent of opinion α, thus decreasing the number of active links by 1. The changes in the active links apply to α = A and B, with α being B and A, respectively.
Although the initial network has all the nodes with degree k , the evolution processes lead to a spread in the degree among the nodes up to some value k max , characterized by a degree distribution P(k) that formally also evolves in time. Among the nodes of degree k, they may have different number of neighbors of an opposite opinion and thus different number of active links. Let Q(n k , k) be the probability of having n k active links given that a node has k links. We refer to this conditional probability Q(n k , k) as the active link distribution. The distributions P(k) and Q(n k , k) satisfy the normalization conditions k max k=0 P(k) = 1 and k n k =0 Q(n k , k) = 1, respectively. Considering that the target node could have any degree k and any number 0 n k k, a dynamic equation for ρ can formally be written as where L tot = k N /2 is the total number of links and it is a constant.
Equation (1) involves the averages of the form with (· · · ) being n k , n k /k and n 2 k /k. The last form emphasizes that two averages are invoked: one over different values of n k for a given k using Q(n k , k) and another over different values of k using the degree distribution P(k). The quantity n k is the average number of active links of a node and thus N n k /2 is the total number of active links, which is also given by ρ L tot . Therefore, n k = ρ k . The difficulty in closing the equation lies in the averages n k /k ≡ U and n 2 k /k ≡ V . In terms of these averages, equation (1) Formally, the average U involves where ρ k = n k /k n k = n k n k /k is the active link density of the agents with k neighbors. Similarly, the average V involves where ρ n k = n k /k. These averages can be evaluated provided that Q(n k , k) and P(k) are known accurately.
In practice, these averages are approximated in terms of ρ so as to close the equation. In [1], it was assumed that Q(n k , k) follows the binomial distribution. Taking ρ k as the probability that a neighbor of a node of degree k has an opposite opinion, the mean n k n k = k n k =0 Q(n k , k)n k = kρ k and the second moment n 2 k n k = k n k =0 Q(n k , k)n 2 k = k 2 ρ 2 k + kρ k (1 − ρ k ), where the last term is the variance within the binomial approximation. The binomial closure scheme [1] then assumes U bino = ρ k k ≈ ρ and V bino = (k 2 ρ 2 The approximation amounts to assuming ρ k to be weakly correlated with the degree k and to itself. We call it homogeneous pair approximation scheme or degree-free binomial closure. Substituting the binomial approximations of U bino and V bino into equation (1) gives an equation for ρ. Setting dρ/dt = 0 gives the equation which can be solved for the steady state value of ρ for given p and q. The solutions to equation (6) for the VM (q = 0) and the MVM (q = 0) are shown in figure 1. The deficiency of the binomial closure scheme by assuming a binomial distribution of active links around a node according to ρ is apparent. In both cases, the binomial closure overestimates ρ by much for a wide range of p. It also overestimates the value of p that ρ drops to zero for the VM and drops to a negligible value for the VMV. Even for small p, the binomial closure does not give the right values. These discrepancies ask for a careful examination of the approach and a better analytical treatment.

Validity of the mean-field approach and binomial closure
There are two possible sources of errors in the steps that lead to equation (6). One is to question the arguments that lead to the rate equation (equation (1) or equation (3)) and the other is the validity of the homogeneous pair approximation scheme. For the first question, we note that U and V are quantities that can be extracted from the simulations. Taking the values U sim and V sim from the simulations and evaluating ρ via give results that are in excellent agreement with the simulations, as shown by the solid line in figure 2. This agreement justifies the approach that has led to equation (3). The validity of the homogeneous pair approximation scheme to Q(n k , k) can also be tested numerically. In figure 3, we show Q(n k , k) for different values of k as extracted directly from the simulations (open squares). The results are compared with the binomial form of Q(n k , k) (circles in figure 3), which can be evaluated by using the simulation result of ρ and the binomial distribution. The binomial Q(n k , k) in general deviates from the exact Q(n k , k), with better agreement only observed for k ≈ k . We could further consider that the densities of the active links are different for nodes of different degrees. In this case, we can take the values of ρ k from the simulations and use ρ k to calculate a degree-dependent binomial Q(n k , k) as shown by the triangles in figure 3. It is clear that the k-dependent binomial approximation fits the simulation data at low k region and k ≈ k well, but deviates from them at high k. The invalidity of the homogeneous pair approximation scheme is further illustrated by the difference between U sim and V sim from the simulations and U bino and V bino as assumed by the degree-free binomial closure  (7) using U and V obtained from the simulations. For comparison, the dashed line gives the results from the homogeneous pair approximation scheme (equation (6)). The results are obtained for a network with mean degree k = 4. Figure 4 shows these differences over the whole range of p. Ideally, U and V should vanish. However, U bino and V bino are large for a wide range of p, with particularly large deviations in the range near the value of p that ρ drops to a negligible value. We conclude that equation (3) and thus equation (7) are valid for theoretical calculations. The problem is how to obtain accurate values of U and V theoretically.

An improved mean-field theory
Equation (1) focuses on the active link density ρ = n k / k and requires accurate forms for the two distributions Q(n k , k) and P(k). To go beyond the spatial correlation included in considering the active link density, we formulate an improved analytical approach by tracking the time evolution of every possible local environment. This treatment has been used to study the time evolution of the epidemics on network [33], and applied to the VM by different updating process recently [20]. It gives better agreement with the simulation data than the homogeneous first order approximation. Here, we use the method to study MVM. At an instant of time, the agents can be classified by their opinion α = A or B and their local environments, i.e. the number of neighbors k and the number of neighbors l of opinion B. Thus, every agent is described by α k,l with α = A, B. Let N A k,l (N B k,l ) be the number of A-nodes (B-nodes) with l neighbors of opinion B among his k neighbors in the system. These variables vary with time, but they obey the sum rule k,l (N A k,l + N B k,l ) = N . The active link density is given by ρ = k,l (l N A k,l + (k − l)N B k,l )/(N k ).   For every dynamic process in the MVM, there are corresponding changes in the set of variables N A k,l and N B k,l . The rate equations are given by where the terms on the right-hand side represent the contributions of different dynamic processes in the MVM, as we now explain.
(i) Consider the process that a randomly selected node of A k ,l switches opinion to become B k ,l as shown in figure 5(a). The probability of selecting a node of A k ,l is N A k ,l /N . The two ways for a node of A k ,l to switch his opinion are: mutation with a probability q and randomly selecting a neighbor of opinion B and switch to B with a probability (1 − q)(1 − p)l /k . The switching of a node of A k ,l to B k ,l causes a local variation of −1 in N A k ,l and +1 in N B k ,l . However, the event will also vary the local environment of the node of A k ,l , since there are k − l A-neighbors and l B-neighbors. Formally, we do not know the exact neighborhood of the neighbors of the switching node A k ,l . Here, we assume that the nodes are connected to each other randomly in the system. A node of A k,l with k − l A-neighbors has the probability (k − l)N A k,l /N A A to be one of the A-neighbors of the A k ,l node, where N A A = k,l (k − l)N A k,l is the total number of A-A links. The average number of A k,l neighbors of a A k ,l node is (k − l )(k − l)N A k,l /N A A . Thus, the switching of a A k ,l node to a B k ,l node leads N A k,l to decrease by (k − l )(k − l)N A k,l /N A A . Correspondingly, the number N A k,l+1 increases by (k − l )(k − l)N A k,l /N A A . Similar considerations can be made for the neighbors of the A k ,l node with a different neighborhood and a different opinion. The total variations in N A k,l and N B k,l as a result of an A k ,l → B k ,l switching are given by where N B A = k,l (k − l)N B k,l and δ i j is the Kronecker delta function. In general, the terms in dN A (B) k,l /dt should involve local information that extends to the next-nearest-neighbors of the node taking actions. In order to close the equations, such longer spatial correlations are approximated in terms of N A k ,l and N B k ,l . (ii) Consider the process where a randomly selected node of A k ,l cuts a link to a B-neighbor and rewires to another A-node as shown in figure 5(b). This happens with a probability (1 − q) pl N A k ,l /(N k ) and results in a change in N A k ,l by −1 and a change in N A k ,l −1 by +1, since the selected node A k ,l becomes A k ,l −1 after rewiring. If the B-neighbor whose link is cut by the A k ,l node is a B k,l node, then N B k,l will change by −1 and N B k−1,l will change by +1. The rewiring process gives a node of A k,l a probability N A k,l /N A to be connected to when an A-node cuts-and-rewires, where N A = k,l N A k,l is the total number of A-nodes. If this happens, there will be a change in N A k,l by −1 and a change in N A k+1,l by +1. The total variations in N A k,l and N B k,l as a result of an A k l → A k ,l −1 cut-and-rewire process are given by (iii) Consider the process where a randomly selected node of B k ,l switches opinion to become A k ,l . The consideration is similar to the process (i) discussed above. The total variations in N A k,l and N B k,l due to a B k ,l → A k ,l switching are given by where N B B = k,l l N B k,l is the total number of B-B links. (iv) Consider the process where a randomly selected node of B k ,l cuts a link to an A-neighbor and rewires to another B-node. The consideration is similar to process (ii) discussed above. The total variations in N A k,l and N B k,l as a result of a B k l → B k ,l +1 cut-and-rewire process are given by where N B = k,l N B k,l is the total number of B-nodes. Equations (9) and (10), together with equations (11)-(18), form a set of coupled equations that can be solved by numerical integrations. Practically, there is a spread in the degrees of the nodes and thus a truncation at an upper cutoff k max is invoked [29]. To initialize the iterations, we assume a network of 50% A-nodes and 50% B-nodes randomly connected with a uniform degree k = k . Figure 6(a) shows the results of ρ (solid line) obtained by the set of equations, for the MVM with k = 4 and q = 0.001. The corresponding results for the VM (q = 0) are shown in figure 6(b). In both cases, the results (solid lines) are in better agreement with the simulation data than the binomial closure (dashed lines). The improved theory captures the behavior in the region of small p very well and gives a p c closer to the simulation results. By using the improved theory, the quantities U and V can be evaluated by The differences between the U and V from simulations and U imft and V imft are shown in figure 6(c) for the MVM. The differences are smaller than those of the theory by using the homogeneous pair approximation scheme (see figure 4). The reason is that the IMFT considers explicitly the coupling to the nearest neighbors and makes the random approximation for the next-nearest neighboring coupling, while the homogeneous pair approximation scheme approximates the coupling to the nearest neighbors. Now, we can see that even if one considers exactly the coupling up to the nearest neighbors, the IMFT still deviates much from the simulation data (see figure 6). It indicates that at q = 0 or very small q, the correlation between the active links n k and the degree k is long-range near the critical region where p c ∼ 0.38 in the VM for k = 4. The IMFT considers explicitly the coupling to the nearest neighbors and approximates the effects of the next-nearest neighbors on the center-node opinion by random distribution. The extent of correlation in the IMFT turns out to be insufficient to give accurate results when the system is near the critical region. When the probability q increases, the mutation events will destroy the tendency of the clustering of the  same opinions and lead to a distribution of the active links that is closer to random. Thus, we expect that a larger q would weaken the coupling between the degree of a node and the active links among the neighbors of the node. Figure 7 shows the result of q = 0.1. In this case, both the binomial closure and the IMFT give results that are in good agreement with the simulation data (see figure 7(a)). The deviations U and V of the two approximating schemes become small. Of course, the IMFT works better than the homogeneous pair approximation scheme as it includes the short-range correlation more precisely.  (19); and (b) simulation results of V sim can be fit to the form in equation (20). The parameters are given in the text.

An empirical approach of modeling U and V
The IMFT invokes a set of variables N A k,l and N B k,l and their time evolution. An alternative approach is to approximate U and V in equation (3) better by using fewer parameters. Here, we introduce an empirical approach that combines a reasonable expectation on how U and V should depend on ρ with fitting to the simulation data. As U is formally given by U = ρ k k , it involves an average over one quantity and therefore we assume Similarly, V is formally k ρ 2 n k n k k and it involves an average over a product of three factors. Anticipating some correlation between ρ n k with k and with itself, we assume a polynomial form The parameters a, b, c and d are to be determined by fitting equations (19) and (20) to U sim and V sim , respectively. Figure 8 shows that the assumed forms indeed fit the simulation results very well. The parameters for k = 4 and q = 0.001 are found to be a = 0.933, b = 3.946, c = −0.578 and d = 1.770. A consistency check is to substitute U emp and V emp into equation (7) and compare the calculated ρ( p) with the simulation data. The results, as shown in figure 9(a), are in very good agreement, with small discrepancies right at p = 0 and p ∼ 0.4, where the transition to a small ρ regime occurs. The difference U emp = U sim − U emp and V emp = V sim − V emp are small over the whole range of p, as shown in figure 9(b). Again, there are some discrepancies at p = 0 and p ∼ 0.4. Recall that for the VM with q = 0, the system enters an inactive state with ρ = 0 for p > p c ≈ 0.38. Typically, fluctuations become larger as the system goes closer to p c . This feature carries over to the MVM (q = 0) and leads to the discrepancy near p = 0.4. On contrasting figure 9(b) with figure 4 for the binomial closure scheme, both U and V have been greatly suppressed by assuming the forms in equations (19) and (20).
A by-product of the empirical approach is that we could obtain an expression for p c in the VM by finding the value of p that gives ρ = 0 after setting q = 0 in equation (7). Besides a trivial result of p = 1, the non-trivial solution gives which is insensitive to the parameters b and c. The binomial closure scheme corresponds to the case of a = 1 and d = 1, giving p c,bino = ( k − 2) / ( k − 1), a result reported in [1]. This result gives p c,bino = 2/3 for k = 4, corresponding to a large overestimation of p c ∼ 0.38. In principle, the parameters a and d should depend on the mutation probability q. If we use the parameters a and d as obtained in the MVM for a small value of q = 0.001 as a reasonable representation of the VM, equation (21) gives p c ≈ 0.33, which is much closer to p c ≈ 0.38 for the VM than the binomial closure predicts.

Summary and conclusion
In summary, we studied a co-evolving VM and its variation so as to analyze the origin of the failure of an existing mean-field approach systematically and to provide an insight on how better analytical treatments can be formulated. The instability of the VM makes detailed analysis hard to perform. We introduced a MVM with a tiny probability q that an agent would change his opinion, without any external influence. The q → 0 limit of the MVM recovers the VM. The MVM suppresses the effects of fluctuations and overcomes the instability problem in the VM. The stability of the MVM provides a platform for obtaining accurate numerical results for various quantities, thus making detailed comparisons with different approximations possible and reliable. After establishing a rate equation for the active link density ρ, it was observed that two averaged quantities (U and V ), each involving two different types of averaging process, were important for the accurate determination of ρ. The assumptions behind the popular binomial closure scheme were stated and the deficiency of the approximation was shown by comparison with the simulation results, for both the VM and the MVM. We re-examined the rate equation and showed that it would work well if the quantities U and V could be found accurately. The failure of the binomial closure scheme is explicitly shown by comparing the simulation results with what was assumed within the scheme. The failure implies that using the active link density alone without considering more details of the agent's neighborhood is insufficient. We, therefore, formulated an improved mean-field theory based on the set of variables N A k,l and N B k,l , which give the numbers of A-opinion and B-opinion agents having exactly k neighbors and l of them are B-opinion neighbors, respectively. The rate equations of these variables were written down by considering the effects of the switching and cutting-and-rewiring processes. This amounts to incorporating spatial correlation to the nearestneighbors more precisely than the homogeneous pair approximation and making the random approximation for the next-nearest neighboring correlation. Results of ρ, U and V from the improved MFT show a better agreement with the simulation results for both the VM and the MVM. We also introduced an empirical approach by proposing the functional forms of how U and V would depend on ρ, based on the definition of U and V . The forms were shown to be valid and the four parameters are determined by fitting to the simulation results. The resulting U emp and V emp give an active link density ρ( p) in good agreement with the simulation results. In addition, the parameters also give a better estimate of p c for the transition in the VM than both the binomial approximation and the improved MFT.
We close by remarking that the analysis presented here sheds light on why simple closure schemes in mean-field theories do not work well for evolutionary games in networks and coevolving dynamic processes in networks. Such dynamic processes often induce inhomogeneities in the agents' neighborhoods. A reliable analytical approach requires a careful account of such inhomogeneities, which are often ignored in mean-field approaches that focused on single-agent (single-site) or two-agent (two-site) quantities. The methods used here can be applied to a wide range of co-evolving two-state processes and could be generalized readily to problems involving multiple states.