Absorbing transition in a coevolution model with node and link states in an adaptive network: Network fragmentation transition at criticality

We consider a general model in which there is a coupled dynamics of node states and links states in a network. This coupled dynamics coevolves with dynamical changes of the topology of the network caused by a link rewiring mechanism. Such coevolution model features the interaction of the local dynamics of node and link states with the nonlocal dynamics of link-rewiring in a random network. The coupled dynamics of the states of the nodes and the links produces by itself an absorbing phase transition which is shown to be robust against the link rewiring mechanism. However, the dynamics of the network gives rise to significant physical changes, specially in the limit in which some links do not change state but are always rewired: First a network fragmentation occurs at the critical line of the absorbing transition, and only at this line, so that fragmentation is a manifestation of criticality. Secondly, in the active phase of the absorbing transition, finite-size fluctuations take the system to a single network component consensus phase, while other configurations are possible in the absence of rewiring. In addition, this phase is reached after a survival time that scales linearly with system size, while the survival time scales exponentially with system size when there is no rewiring. A social interpretation of our results contribute to the description of processes of emergence of social fragmentation and polarization.


INTRODUCTION
Modeling a complex system of interacting agents requires specifying the network of interactions and the state of the agents, represented as nodes of the network. The links of the network can also have a state, representing for instance attractive or repulsive interactions. In addition, the network might not be fixed, but adaptive with a time dependent topology. A general dynamical model includes the coupled dynamics of the states of the nodes, the states of the links and the topology of the network. In terms of social collective emergent properties a key question is if the asymptotic state of the dynamics describes some sort of consensus or alternatively, group formation with social fragmentation and/or polarization.
A first class of models are those that do not include link states, other than existence or non-existence of the link between two nodes. Among these, there are models that describe dynamics on the network: evolution of the states of the nodes in a fixed network, while other describe dynamics of the network, as for example dynamics of network formation with no reference to the state of the nodes. We use the term coevolution models [1,2] for those describing a coupled evolution of these two dynamics: node state dynamics and network topology dynamics, so that the structure of the network is no longer a given, but a variable. In these coevolution models the network plasticity is a parameter measuring the ratio of time scales of network evolution and evolution of the states of the nodes. In many cases coevolution models feature a network fragmentation transition [3,4], associated with an absorbing phase transition, in which a single component network fragments into several disconnected components. Coevolution models have been analyzed in the context of social differentiation [5], voter model and opinion formation [3,[6][7][8] or cultural polarization [9]. They have also been used to describe empirical data on community structure of online games [10] and more recently twitter data on echo chambers and polarization dynamics [11] as well as echo chambers on political parliamentary voting [12].
A different class of models are those that focus on the state of the link [13]. While the existence of a link in the context of opinion formation has been often interpreted as a positive interaction like friendship, trust or collaboration, it has also been argued that negative interactions are a major driving force in collective social behavior [14,15]. Therefore, considering positive and negative links should be an essential ingredient of interacting agents models. Links states were pioneered in Heider's social balance theory [16][17][18][19][20][21]. Coevolution of link states and network topology, with no node state dynamics, has also been considered [22]. More attention has been recently paid to the coupled dynamics of nodes and link states in a time independent network (fixed topology) in different contexts [23][24][25][26][27]. In particular, in [27], inspired by processes of opinion formation, we considered an interacting agents model in which each node can be in either of two states and links can also be in either of two states associated with a friendly or unfriendly relation. Two nodes in a different state linked by a friendly relation or two nodes in the same state linked by an unfriendly relation are considered as unsatisfying pairs and they evolve to satisfying pairs by changing the state of one of the nodes or by changing the state of the link. This model exhibits an absorbing phase transition from a dynamically active state with persistent dynamics to a frozen absorbing configuration. The transition occurs for a critical value of the relative time scale for node and link state updates, that depends on the average degree of the network. The system also shows a finite-size topological transition associated with group splitting. Groups are here defined as a set of nodes in the same state connected by friendly links among themselves and by unfriendly links to the nodes in other groups. The existence of these groups, when they are poorly connected among them, reflects social fragmentation, with each group acting as an echo chamber. While [27] considers dyadic interactions, a related model in [26] implements an optimization dynamics of social balance triangular relations.
In this paper we go beyond the model of [27] by also considering adaptive dynamics of the network, so that we have a full dynamical model including dynamical evolution of the states of the nodes and the states of the links as well a time dependent network topology. This is a coevolution model, in the sense that the coupled dynamics of node-link states coevolves with dynamical changes of the topology of the network caused by a link rewiring mechanism. Note also that this model includes a coupling between the local mechanism of changes of node and link states and the nonlocal mechanism of random rewiring in the network. We find in this model an absorbing phase transition which is rather independent of the network dynamics. However, network dynamics gives rise to significant new phenomena. In particular, link rewiring can produce a network fragmentation transition on the critical line of the absorbing transition and it modifies essentially the nature of the ac phase: a consensus phase can be reached by finite-size fluctuations.
The paper is organized as follows. First, we discuss our general coevolution model and its rate equation formulation. Next, we describe the absorbing transition predicted by the rate equations, as well as results on this transition obtained by Monte Carlo simulations. In the following sections we focus on the special case of largest rewiring probability, the approach to the absorbing states and calculation of survival times, and the topological transitions among the final frozen states. The final section summarizes our main results. An Appendix discusses the derivation of the rate equations. We consider an initial uncorrelated random network with N nodes and an average degree µ [28]. Each node is endowed by a binary state value 1 or −1. Each link also holds a binary state variable + or − representing, respectively, attractive or repulsive interactions. In an opinion formation model, the state of the nodes represent two different opinions and attractive or repulsive interactions stand for friendly or unfriendly relations. The states of the nodes are represented in the figures by filled-in blue and white circles, while the states of the links are represented by solid (+) or dashed (−) lines. We consider repulsive links connecting nodes with different (resp. the same) state and attractive links connecting nodes with the same (resp. different) state as satisfying (resp. unsatisfying ) relations. In Fig. 1, we show all 6 possible types of pairs and identify a, e and c as unsatisfying pairs and b, d and f as satisfying ones [27].
We implement a dynamical model that tends to minimize the number of unsatisfying pairs. The change from unsatisfying pairs to satisfying ones is done through two different mechanisms: (i)Flipping node or link states (rightward arrows in Fig. 1): This is a local mechanism in which an unsatisfying pair is turned into a satisfying one by changing the state of one of the nodes or the state of the link. With probability p a link state is changed, and with the complimentary probability (1-p) a node state is changed. This gives rise to a coupled dynamics of node and link states in which p measures the ratio of time scales of node and link evolution. The pair a is turned into a satisfying pair b via a link update with probability p or it is turned into a satisfying pair f through a node update with complimentary probability (1-p), as shown in Fig. 1. Note that in the node update a to f, it doesn't matter which of the nodes in the pair a changes its state, both result in f. The pair c is similarly transformed into a satisfying pair: a link update transforms it into the pair d and a node update (from any side) transforms it into the pair f. For the case of the unsatisfying pair e, a link update with probability p transforms it into the satisfying pair f. With complimentary probability (1-p), a node update takes place on pair e; either its white node with probability 1−p 2 flips its state and e becomes d, or alternatively the blue node with probability 1−p 2 changes its state and e becomes b.
(ii)Link rewiring (leftward arrows in Fig. 1). In this nonlocal mechanism [3], one of the two nodes, chosen at random, of an unsatisfying pair of type e breaks the existing link and reconnects to another node of its own type, again chosen randomly amongst all possible nodes. Hence the unsatisfying pair of type e can either become a satisfying pair of type b or of type d according to which one of the two nodes (white or blue, respectively) is selected to find a matching node in the network. Summing up, when an e-type pair has been chosen, the following actions can occur: (i) change the state of the link from attractive to repulsive, with probability (1 − r)p; (ii) the node holding the blue opinion turns into white, with probability 1−p 2 ; (iii) the node holding the white opinion turns into blue, with probability 1−p 2 ; (iv) the node holding the white opinion breaks the link and binds to another white node in the network, with probability rp 2 , (v) the node holding the blue opinion breaks the link and binds to another blue node in the network, with probability rp 2 .
The parameter r sets the relative time scale of network evolution. In principle one could also think of a rewiring mechanism that does not change link state, but instead turns pairs a and c into a satisfying pair f. We do not include this possibility here: In terms of an opinion formation model it corresponds to a search in the network to establish an unfriendly relation with an agent with an opposite opinion. It can be socially argued that this process would occur with a much smaller probability than the one considered here of searching to establish a friendly relation with an agent with the same opinion.
The three parameters of the model, in addition of system size N , are the probabilities p and r and the average degree of the network µ. The limiting case r = 0 was considered in [27]. In that case there is coupled evolution of node and links states by local dynamics, but there is no network dynamics. In our coevolving case, the network evolves by a nonlocal mechanism coupled to the evolution of node and link states. In the limiting case p = 0, there is only dynamics of the states of the nodes and no network dynamics. Still this limit does not correspond to the well known voter model in a fixed network [29] where updates from pairs a and c to pair f do not occur. In the case p = 1 there is only link dynamics (including rewiring).
In a Monte Carlo implementation of this coevolution model, the system evolves asynchronously. A single pair is selected at random at each time step: If it is satisfying, nothing happens. If it is unsatisfying, it adapts according to the dynamics described in Fig. 1. A Monte Carlo step is counted after a sequence of time steps equal to the total number of links in the network. In an Erdős-Rényi network [30] with N nodes and average degree µ the total number of links is L = 1 2 µN , a quantity that is kept constant during the evolution. Defining L i , i ∈{a, b, c, d, e, f} as the number of pairs of type i, the associated densities are ρ i = L i /L, satisfying the obvious normalization condition The dynamical rules are such that the unsatisfying pairs (ρ a , ρ c and ρ e ) are turned into satisfying. An absorbing configuration of the system [31] is reached when there are no unsatisfying pairs, ρ a = ρ c = ρ e = 0. The system is then in a frozen state where there can be no further evolution. Note that, although at each step of the dynamics an unsatisfying pair becomes a satisfying pair, this does not mean that the total number of unsatisfying pairs decreases, because an update that converts a pair from unsatisfying to satisfying by changing an individual node state might also change the status of another pair, to which the node involved in the update also belongs to, from satisfying to unsatisfying. The main question is when and how an absorbing or frozen configuration is reached. To answer this question we will analyze the dynamical evolution of the densities ρ i (t), i ∈{a, b, c, d, e, f}. This dynamics is analyzed by Monte Carlo simulations and by a set of rate equations.

RATE EQUATIONS
The densities {ρ a , ρ b , ρ c , ρ d , ρ e , ρ f } obey six coupled differential rate equations (RE) describing the model sketched in Fig. 1. These equations are derived in the Appendix. Each equation indicates the time evolution of one of the densities. The rate equations are a combination of linear and nonlinear terms. The linear terms describe the change of densities due to the direct update of pairs shown in Fig. 1, while the nonlinear terms are the outcome of the node update causing a change in the densities of the pairs linked to the updated node. There are two main assumptions in the derivation of these equations to be taken into account when comparing their predictions with Monte Carlo simulations. First, they are derived in the thermodynamic limit (N → ∞), so that they cannot describe finite-size fluctuations. Second, they assume a random regular network in which all nodes have exactly µ random neighbors.
The rate equations have two sets of fixed points or stationary solutions. The first solution is given by , i.e. there are no unsatisfying pairs, and the relative fractions of satisfying pairs are arbitrary, but still need to fulfill the normalization condition The first solution, Eq. (1), described by generally non-vanishing values of the densities that depend on on r, p and µ, is approached asymptotically on time and it is independent of initial conditions. Given that in this solution there exist unsatisfying pairs, it corresponds to a dynamically active state. The second set of solutions, Eq. (2), includes a continuous of solutions which depend on initial conditions. However, in all of them there are no unsatisfying pairs, and therefore they correspond to absorbing states with frozen dynamics.
The condition ρ st e ≥ 0 ensures that all other densities are non-negative as well. Therefore p = p c is the critical value for an absorbing transition with order parameter ρ st e : For p ≤ p c (µ, r) the system reaches a stationary dynamically active phase in which ρ st e is nonzero for all values of r, µ and p. In this regime a linear stability analysis indicates that this solution is linearly stable. On the other hand, for p > p c (µ, r) the system reaches an absorbing or frozen state which turns out to be marginally stable. This analysis shows that the absorbing phase transition found in [27] for r = 0 is robust under rewiring network dynamics.
In order to check our theoretical predictions of the active-to-absorbing transition we have conducted extensive Monte Carlo simulation of our coevolution model. The simulations have been carried out taking as initial condition an Erdős-Rényi network with the desired average connectivity µ. Initial conditions for the states of the nodes and states of the links are chosen randomly distributed with initial densities x 0 and l 0 , respectively, of blue nodes and attractive links (see the Appendix ). Results of the simulations and the numerical integration of the rate equations are shown in Fig. 2 for three value of the rewiring probability (r = 0, 0.5, 1) as a function of µ for p = 0.8 fixed (left panel), and as a function of p for µ = 6 (right panel). The Monte Carlo simulations confirm the predicted absorbing transition, but quantitative agreement for the order parameter ρ st e in the active phase is better far from the transition point, that is, for large values of µ or small values of p. In addition, the transition point is shown to be rather independent of r. This fact is not captured by the approximate rate equations that, as mentioned earlier, assume a random regular network at all times. A comparison of the MC and RE results for the critical line in the (p, µ) plane for different values of r is also shown in in nodes in different states. This is a relevant quantity to monitor the effect of the network dynamics as measured by r: its value for all points in plane (p, µ) decreases when increasing the rewiring parameter r in such a way that in the active phase it becomes zero at r = 1. Moreover, the analytical solution of the rate equations, Eq. (1), predicts that the value of ρ st f at the critical line p = p c is independent of µ and it is only a function of r. At the critical value of p, we find independent of µ. Replacing this value of q c in the expression for ρ st f in Eq. (1) we find the specially simple relation This result, well confirmed by the MC simulations, is shown in Fig. 4, and it indicates a special behavior of the system for r = 1.
To be more precise, as we will discuss latter the mechanism of network fragmentation transition is rooted in the absence of pairs of type f at r = 1 where ρ c f r=1 = 0.

THE CASE R=1
When r = 1 the unsatisfying pair e only evolves by node state flipping or by link rewiring, but no link state flipping is allowed. This turns out to be a singular situation in several aspects of our study. Firstly, we note that Eq. (1) implies that in the active phase the densities ρ st a , ρ st c and ρ st f are zero. Therefore only one density of unsatisfying pairs (ρ st e ) and two equivalent densities of satisfying pairs (ρ st b and ρ st d ) survive in this active phase. These non-vanishing dynamical variables exhibit very large fluctuations as seen in 10 individual realizations of the process obtained by MC simulations (Fig. 5). Starting from a fixed initial condition, fluctuations of ρ b and ρ e increase as time goes on. A quantitative measure of such fluctuations is given by the standard deviation of ρ b at stationary state, σ[ρ st b ]. We show in Fig. 6 the dependence of this quantity on the rewiring probability, r, for two different system sizes. As illustrated, there is an abrupt increase of these fluctuations at r = 1, rather independent of system size.
Secondly, we note that in the absorbing phase it is always ρ only attractive links, or (ii) ρ st b = 0 and ρ st d = 0, which represents a fragmented network with disconnected groups, each of them with nodes in a given state connected by attractive links. This second solution, predicted by the rate equations in the thermodynamic limit, is a network fragmentation as a manifestation of criticality which occurs for r = 1 as discussed in the last section on Topological transitions.

DYNAMICS: APPROACH TO THE ABSORBING STATES AND SURVIVAL TIMES
In the absorbing phase the system approaches an absorbing state from given initial conditions by an ordering process. This process is described by the time evolution of ρ e as shown in Fig. 7 for a set of parameters in the absorbing phase and three different values of the rewiring parameter r. We observe an exponential approach to the absorbing state, ρ e ∼ e −t/τ0 . We also observe that τ 0 is independent of r, so that the rewiring process does not modify the approach to the absorbing state. In a finite system of size N the exponential approach lasts until ρ e ∼ 1/N . This occurs at a characteristic time τ ∼ log N . This characteristic logarithmic dependence is shown in panels (a) and (c) of Fig. 8 for three values of the rewiring probability and two sets of system parameters in the absorbing phase. Again, the logarithmic dependence is not modified by the rewiring process, although the value of τ decreases with increasing r, indicating that the rewiring process accelerates the ordering process.
In the active phase of a finite system the dynamical state has a survival time because a finite-size fluctuation will eventually take the system to an absorbing state. This survival time τ is a stochastic variable whose average value is shown in panels (b) and (d) of Fig. 8 as a function of the system size for a set of parameters in the active phase and different values of r. The curves for r = 0 can be well fitted by an exponential dependence τ ∼ e αN , while for r = 1 we find a linear dependence τ ∼ βN . For intermediate values of r the scaling of τ with N interpolates between these two extreme dependencies. For given r, we conjuncture that the system size scaling of the survival time can be fitted by a phenomenological function as τ = α 1 (α 2 N ) r e (1−r)α3N . The chosen points (system parameters) of panels (b) and (d) in the plane (p, µ) have the same distance from the critical line, and for both of them we find common fitting coefficients, α 1 , α 2 , α 3 . Therefore, in the active phase, as compared with the absorbing phase, the rewiring process has a much more significant effect in the dynamical processes changing the system size scaling. In addition a special behavior, the linear system size scaling of the survival time is found for r = 1.

TOPOLOGICAL TRANSITIONS
In a finite system, either in the absorbing phase or after the survival time in the active phase, the system reaches a fully satisfying configuration in which ρ a = ρ c = ρ e = 0. These configurations display a range of different non-trivial topological structures depending on the system parameters (p, µ, r). In this section we describe these configurations, the transitions among them and the relation with the active to absorbing transition discussed before.
The final, dynamically frozen, configurations can be classified according to the vanishing or non vanishing values of ρ b , ρ d , ρ f . There are two main categories, fragmented and connected networks. The connected networks can be organized into three subclasses: consensus (CP), two-group (TP) and split configurations (SP). In the fragmented configurations it is ρ f = 0, ρ b = 0, ρ d = 0: therefore the nodes organize in two disjoint networks, one network containing all blue nodes and the ther all white nodes. All links within the nodes of each independent network are attractive and there are no links between nodes belonging to different networks. An example of a fragmented configuration is shown in Fig. 9(a). In the connected consensus configurations, still ρ f = 0, but now either ρ b = 0 or ρ d = 0: all nodes are in the same state connected by attractive links as shown in Fiq. 10(b). Both two-group and split configurations are characterized by ρ b = 0, ρ d = 0, ρ f = 0. In a two-group configuration, two groups of nodes, each group in a different state, are connected by some repulsive links, while the internal links of each group are attractive. An example is shown in Fiq. 11(a). The split configurations are similar to the two-group configurations except by the fact that some very small groups are attached to two large groups by some repulsive links as shown in Fiq. 11(b). We first consider which configurations are obtained for the special case r = 1 for which we had interesting predictions from the rate equation analysis. The solution of the rate equations for r = 1 in the active phase, and also at the critical line, predicts that ρ f = 0, which corresponds to either a consensus or a fragmented configuration. The Monte Carlo simulations conclude that the fragmented network solution is only found at the critical line of the absorbing transition, see an example in Fig. 9(a). Further evidence of this result is shown in Fig. 12, where we plot, for different system sizes, the probability of occurrence of the fragmented structure in the (p, µ) parameter plane, showing that this probability is only non-zero at the critical line and it increases with system size. An example of how this network fragmentation solution disappears as we move away from r = 1 is shown in Fig. 9(b). On the other hand, the consensus configuration is the one found in the active phase. An example of this configuration is shown in Fiq. 10(b), together with a snapshot of the dynamically active network before it falls into the consensus absorbing state by a finite-size fluctuation. Evidence that, for r = 1, the consensus solution is the only one found in the active phase, and that it also exists at the critical line, is given in Fig. 13(a) that shows results of Monte Carlo simulations for the probability of occurrence of the consensus configuration in the (p,µ) parameter plane. As we move slightly away from r = 1, e.g. Fig. 13(d) for r = 0.98, the consensus configuration disappears for values of (p,µ) in the active phase close to the critical line.
In the absorbing phase for r = 1, as well as in the absorbing and active phases when r < 1, two-group and split configurations are found. The probability of occurrence of these solutions in the (p,µ) parameter plane are shown for r = 1 and r = 0.8 in Fig. 13. We observe a transition between two-group and split configurations, named Finite-size topological transition for r = 0 in [27]. Our results indicate that this transition does not exist in the active phase for r = 1, while for r < 1, and for the chosen symmetric initial conditions, it occurs for a value of the network average degree µ split which is essentially independent of p. In Fig. 14 we show the distribution of the number of groups for different parameter values. We use the precise definition of a group as a set of nodes in the same state, connected by attractive links among themselves and by repulsive links to the nodes of other groups. From this figure we estimate µ split ∼ 8 which, as compared with values for r = 0 [27], indicates that µ split decreases as a consequence of network dynamics (r = 0).

SUMMARY AND DISCUSSION
We have analyzed, by an analytical rate equation approach and by Monte Carlo simulations, a general dynamical model of interacting agents in a network in which nodes can be in two different states (blue or white) and links are either in an attractive or repulsive state. Unsatisfying pairs with repulsive links joining nodes in the same state, or with attractive links joining nodes in different states evolve to satisfying pairs by two different coevolving processes: i) Coupled dynamics of the state of the nodes and the state of the links; ii) network dynamics by random rewiring of an attractive link between two nodes in different states (pair e) to connect to a node in the same state (pairs b or d ). The rewiring parameter 0 ≤ r ≤ 1 is such that there is no dynamics of the topology of the network for r = 0, while for r = 1 pair e only evolves by updating the state of the nodes or by link rewiring, but there is no change of the state of the link from attractive to repulsive. We find a phase transition, predicted by the rate equations, between an active phase with persistent dynamics and an absorbing phase in which the system has reached an absorbing configuration as steady state. This phase transition is already found for r = 0 [27], but its existence is robust against the introduction of the network dynamics: The origin of the transition and the critical line do not significantly depend on the rewiring parameter r. However, the network dynamics in the limiting case r = 1 produces very significant physical changes at criticality and in the active phase. As predicted by the rate equations, we find network fragmentation on the critical line, and only on this line, so that fragmentation is a manifestation of criticality. In the active phase only links in the attractive state survive and the fluctuations of the density of attractive links connecting nodes in the same state diverges as r −→ 1 indicating the singularity of the case r = 1. Finite-size fluctuations in the active phase take the system, after a survival time, to a consensus frozen configuration with all nodes in the same state and connected by attractive links, while for r = 1 there are other possible configurations with a finite-size topological transition between two-group and split configurations. In addition, the survival time scales linearly with system size N , while in the absence of link rewiring it scales exponentially with N .
We have found four possible final frozen configurations of the system: fragmented, consensus, two-group and split. When our model is understood as an opinion formation model, these configurations, except the consensus one, can be interpreted as a manifestation of social polarization. Our results give support to the claim that negative interactions are a major cause of social polarization. Indeed, our solution for r = 1 in the active phase is such that repulsive links do not survive and the consequence is that the final state is that of consensus.

APPENDIX
Rate equations for the coupled evolution of node and link states and rewiring process The dynamical behavior of the different densities of pairs as a function of time can be obtained as a set of so-called rate equations whose derivation we now sketch. The basic ingredient is the evolution rules sketched in Fig. 1 as considered on a network in which each node has exactly µ links at all times. This assumption is a clear approximation because in the rewiring process the degrees of the nodes change over time. With this provision in mind, the rate equations are written as The only pairs chosen by the update rule are a, c and e. However, the update of one of these pairs will also change the number of pairs b, d and f , as in the process of node update, the states of links connected to the updated node will also changed. The nonlinear terms in the rate equation account for this mechanism. In the following we explain the derivation of the rate equations considering the linear and nonlinear terms separately.
The linear terms are obtained by the variation of densities due to the direct update of the state of nodes and links in a time step as where dρi dt i→j is the direct effect of the update from pair i ∈ {a, c, e} to the j ∈ {b, d, f }.
As an example, we now we derive in detail the linear terms of the equation for ρ e : In a given update step, with probability ρ e a pair e is randomly chosen and updated producing a density changes of ∆ρ e = − 1 N µ/2 . According to the local adaptation rule in Fig. 1, with probability (1 − r)p, pair e turns into the pair f. With probability 1−p 2 , the pair e becomes b and with probability 1−p 2 , it becomes d. In addition, according to the rewiring update rule in Fig. 1, with probability rp 2 , the pair e becomes b and with probability rp 2 , it becomes d. The time interval (measured in Monte Carlo steps) in any update step is given by ∆t = 1 N µ/2 . Therefore, the variation of ρ e due to the direct effect of update of the pair e is Non-linear terms are an indirect effect of the node update. Blue nodes can be an end to any of the links c, d, e and f while white nodes can be an end to any of the links a, b, e and f. Thus, node update will change the state of the connected links to the updated node so that As an example of a specific case of how to obtain the nonlinear terms, we derive the non linear term −(1 − p)(µ − 1) 1 2 ρeρ f 2χ in the equation for ρ f . This term is the consequence of a node update from the pair connection e to d. In a Monte Carlo step, with probability ρ e a pair e is randomly picked and with probability 1−p 2 a node update to became d takes place. Now, let us examine the change of the rate of ρ f under the update of e to d as presented in Fig. 15. The normalized number of whole pair connections attached to the one side of a pair e is µ−1 N µ/2 and from this portion, the fraction of the pair f that is attached to the link e is ρ f χ µ−1 N µ/2 . In addition, due to the asymmetry of pairs e and f, the update from any side of these pairs would result in a different pair. For instance, if the white opinion in the pair f is flipped it converts to c, while if the blue opinion is flipped, it turns to a. Thus, when we deal with pairs e and f, the contribution of nonlinear terms appear with probability 1 2 . Thus, the global change in the density of pairs f due to the node update from e to d is given by ∆ρ

Initial condition
In the numerical integration of Eqs. (6) it is important to ensure that the initial conditions satisfy the relations where µ b and µ w stand for the mean degree of the blue and white nodes, respectively. Note that, µ b and µ w are functions of time but the total number of links in the network, N µ 2 , is a conserved quantity. In the numerical integration of the rate equations and the MC simulation we choose regular random and Erdős-Rényi networks as initial configurations, respectively. In both cases, requirements on initial conditions are satisfied taking where 0 is the initial fraction of attractive links and x 0 the initial fraction of blue nodes. Finally, in the numerical integration of the rate equations, we employed the predictor-corrector method.