Coevolution of agent’s behavior and noise parameters in majority vote game on multilayer networks

In many real-world systems, agents’ behaviors are usually coupled with environmental changes. To model their coevolutionary process, we study a non-equilibrium model known as majority vote model coupled with reaction-diffusion processes on a two-layer multiplex network. The dynamics of the noise parameter in noise layer is related to the voting behavior and represents the increase and decrease of social tension in the context of social dynamics. We perform Monte Carlo simulations and finite-size scaling analysis in order to investigate the statistical behavior of the model. It is interesting to find that our coupling mechanism induces a continuous order–disorder phase transition on random regular graphs, but the critical phenomenon disappears on square lattices. Besides, switching from one sign of the spontaneous magnetization to the opposite is also observed near the critical region. Given specific parameters, the system may self-organize to a critical state from any initial conditions. In addition, a mean-field method is developed to study the properties of the phase transition analytically and the solutions are in good agreement with our numerical results qualitatively.


Introduction
Besides physics [1,2], chemistry [3], biology [4] and other natural science, pattern formation is also an important topic in social dynamics. It focuses on the transition from an initial disordered configuration to a final ordered state in social phenomenons [5] such as opinion formation [6][7][8][9], cultural dynamics [10], language dynamics [11,12], etc. In recent years, many studies have been published in modeling social dynamics through methods of statistical physics. Take opinion formation as an example, a binary variable like an Ising spin is often introduced to represent two opposite opinions. The update rules of these models may differ, but the order-disorder transition is always the common theme in social dynamics [5].
Many models have been proposed to model opinion formation, among which the isotropic majority vote model [13] is a well-known non-equilibrium model and much attention has been attracted to it recently. The isotropic majority vote model is defined on a square lattice where each site is occupied by a spin variable. At each time, a randomly chosen spin flips with probability q if it agrees with the majority sign of the spins in its neighborhood. There is a continuous order-disorder phase transition in majority vote model and it has proven the fact that non-equilibrium stochastic systems with up-down symmetry fall in the universality class of the equilibrium Ising model [13]. Since then, many investigations on the majority vote model have been performed and a lot of work focus on the effects of underlying topologies, such as random graphs [14,15], scale-free networks [16], small world networks [17,18] and other complex networks [19][20][21][22][23]. The continuous orderdisorder phase transition is observed in different kinds of network topologies but with varied critical noise q c and critical exponents [24]. Besides, modifications such as multi-state [25][26][27], different agents [28] and the effect of individuals inertia [29] have also been considered by researchers.
Recently, multilayer networks, also known as multiplex networks, i.e. a set of interconnected networks where each layer represents a distinct feature, have been widely explored by many researchers [30][31][32]. The reason why multilayer networks become a hot trend in network science is that many real-world complex systems Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. are composed of several constituents or dynamics and the traditional complex network approach may result in not fully capturing the details of our concerns [33]. Thus it is natural to model a wide variety of complex phenomena through multilayer networks. Of particular interest in the study of multilayer networks is the coevolution between different dynamical processes. Nicosia et al [34] reported the emergence of spontaneous explosive synchronization in a coupled multilayer network, which will not be observed without coupling. In coevolving voter model, the states of the nodes and the network topology coevolve with each other and exhibit fragmentation transitions [35,36]. In the study of coevolutionary games on multilayer networks, different types of interdependencies between layers give rise to fascinating evolutionary outcomes [37][38][39][40][41][42][43].
It should be mentioned that previous studies on majority vote model mainly concern the phase transition and network topology in a single layer network, ignoring a very important fact that in real social systems the noise parameter q which reflects the degree of a person's rational is often heterogeneous. Although Vilela et al [44] considered a binary noise distribution and Vieira et al [45] investigated independent behavior as a second noise, it is worth noticing that the noise parameter q may be dynamical itself and coevolve with the dynamics of spin variables.
In this work, we take the original majority vote model as a two-layer multiplex network in a natural way. We add a reaction-diffusion process in the dynamics of the noise parameters. The dynamics of spin variables and noise parameters are coupled and coevolving with each other. We first introduce the model and then explore the statistical behavior of the model by Monte Carlo simulations. A finite-size scaling analysis is also applied to investigate the critical behavior. In section 4, a mean-field method is employed to study the qualitative dynamical properties of the system. Finally, some concluding remarks are given.

Model
First, we describe our model defined on a two-layer multiplex network. The first layer, which we call the vote layer, is the original majority vote model defined on a complex network. We assign a binary spin variable σ i =±1 to each node and the configuration of the system is denoted by {σ i }. The second layer, which we call the noise layer, is also defined on a complex network, with each site assigned to a noise parameter q i (0q i 0.5).
There is a one-to-one correspondence between a spin variable σ i in vote layer and its noise parameter q i in noise layer. In this way, a given spin σ i flips with probability x 0 and S(0)=0. The summation is over the neighborhood of the site and a ij =1 if site i and j are connected. A schematic diagram of our model is presented in figure 1.
The noise parameter q i reflects the degree of a person's rational. It describes social tension that affects human behaviors. In isotropic majority vote model, the noise parameters are homogeneous, which means q i =q for every site. In our model, however, q i evolves according to a reaction-diffusion equation: where D is the diffusion coefficient. The diffusion term  ( ) D q t r, 2 describes the regional relationship in social environment. The reaction term s ( { }) f q, i depends on ( ) q t r, not only, but the configuration {σ i } in vote layer as well. In our model, the Figure 1. A schematic diagram of our model of N=5 nodes with majority vote dynamics in vote layer and reaction-diffusion process in noise layer. There exists a one-to-one correspondence between a spin variable σ i and a noise parameter q i . reaction term has two mechanisms similar to the self-organized sandpile model [46,47]. One is an autonomous decreasing at a constant rate a, which means the increase of social tension. It is similar to the addition of grains in sandpile model. The other one is the release of social tension with the formation of collective structures in vote layer. It is also similar to the avalanches when the sandpile becomes too steep. This mechanism has taken the second order consequence into account, that is, with the formation of clusters, the government will take actions to reduce the social tension to avoid further civil disturbance. This term is given by bC N i in numerical simulation where C i is the corresponding cluster size of site i in vote layer and N is the number of nodes. Then we have Obviously, the value of r=b/a that reflects the intensity ratio of two mechanisms dominates the evolution of the system with an adequate value of a fixed. If noise parameters in noise layer become too small, the formation of clusters in vote layer will tend to increase it. On the other hand, if noise parameters become too large, the cluster size will also be small and the constant decay will tend to decrease it. Therefore, a constant decay and external flux of noise passing through the noise layer will organize the system itself into a non-equilibrium steady state.
In the following section, we introduce the numerical implementation of our model described above. We have performed the Monte Carlo simulations as follows: (i) Initialize vote layer randomly with a fixed noise parameter q initial for every site.
(ii) We first update the configuration of vote layer and then update noise parameters in noise layer. However, we are dealing with two different dynamical processes that may have different time scales. Here, we first choose a spin variable randomly and update its state according to (1). One Monte Carlo step (MCS) is accomplished after N spins are updated on average. After τ MCS, we then update noise parameters according to (3) simultaneously.
(iii) Quantities of interest are calculated for every update of noise parameters after the non-equilibrium steady state is reached.

Numerical results
We have simulated our model on random k-regular graphs in both vote layer and noise layer with various sizes N. In simulations, we choose a fixed value of a=1×10 −3 and D=1×10 −3 . From the model, it could be found that the value of a and time scale τ have some impacts on the dynamics and statistical properties of the system, but have fewer effects on the final steady state on average. In particular, we are interested in the magnetization M N and mean value of noiseq defined as follows: where á ñ  stands for an ensemble average and r=b/a. We waited 1×10 4 τ MCS to reach the steady state with a time scale τ=5 and the time averages were recorded for the next 1×10 4 τ MCS.
In figure 2, we show the magnetization M N and mean value of noiseq as a function of r=b/a on a random -4 regular graph in both vote layer and noise layer with N ranging from 256 to 8192. As the value of r increases,q that reflects social tension also increases and the order-disorder phase transition occurs.
In figure 3, we investigate statistical behaviors of our model near the critical region. We plot the probability density function of magnetization m for different values of r. We notice that there is a transition from a normal distribution to a bimodal distribution with the decreasing of r. It indicates that there is a phase transition from disordered state to ordered state. Besides, we report a transition from two ordered states in bimodal distribution region, like Ising model on a square lattice slightly below the critical temperature. It is an indication of the existence of a bi-stable state in the system.
We have also examined the effects of time scale τ and diffusion coefficient D in our model. Although the time scale parameter τ does not influence the final average magnetization M N , it influences the dynamics and the transition rate from the two ordered states near critical region. With the increase of τ, the transition becomes more frequent, as shown in figure 4. As for the diffusion coefficient D, it has a great impact on the distribution of noise q after the steady state is reached. With the increase of D, the distribution of noise q in noise layer becomes more homogeneous, as shown in figure 5(a).
Besides, the network topology must have some effects on the dynamics of the coupled system. We have applied random regular graphs with different degree k and square lattice. The dependence of the magnetization M N with the value of r relies heavily on certain topologies in vote layer, as shown in figure 5(b), although topologies in noise layer have little influence on the dependence. The magnetization M N as a function of r for square lattices in both vote layer and noise layer shows a qualitative difference from the case of random regular graphs. With detailed investigations of finite-size scaling analysis in the next section, we study the existence of continuous phase transition in our coupled systems.

Finite-size scaling analysis
Here we study the critical behavior of our model by finite-size scaling analysis. We are interested in susceptibility χ and the Binder cumulant U defined as follows:  where á ñ  stands for an ensemble average. These quantities along with the magnetization M N (r) obey the following finite-size scaling relations: and the critical exponents β, ν and γ satisfy the scaling relation: b n where D eff is the effective dimensionality. In addition, the value of  ( ) r N where χ N (r) has a maximum also scale with the system size N as where c is a constant.
In figure 6, we show results for susceptibility χ and the Binder cumulant U for random -4 regular graphs in both vote layer and noise layer with several sizes N, which confirms the presence of a continuous order-disorder phase transition in our coupled dynamical systems. Further, the critical parameter r c is estimated as the point where the curves in figure 6(b) intersect. We have r c ∼2.09.
In figure 7, we exhibit results for all the relevant critical exponents. In figure 7(a), we plot the dependence of the magnetization M with the system size N. The obtained critical exponent is β/ν=0.223±0.001 according  to (6a). From (8), we obtained the critical exponent 1/ν=0.47±0.03, as shown in figure 7(b). To obtain an estimate of the critical exponent γ/ν, we plot the susceptibility χ versus the system size N. The obtained critical exponent is γ/ν=0.55±0.01 according to (6b). We can also plot the susceptibility at its maximum versus the system size N. It gives us an estimate of γ/ν=0.556±0.006.
To confirm our estimate for r c , we plot  ( ) r N as a function of n -N 1 . The extrapolation gives r c =2.086±0.001, which is in good agreement with the r c obtained from Binder cumulant. Besides, the effective dimensionality D eff is 0.996±0.012.
In previous studies, the noise parameter q is homogeneous and the scaling relations are functions of q, not the value of r in our multilayer systems. Nevertheless, we point out that our critical exponents β/ν and γ/ν are similar to those obtained on random graphs and scale-free networks with degree k=4 [14,16]. In particular, we find that the critical exponent 1/ν is also similar to that for Stauer-Hohnisch-Pittnauer (SHP) networks [23]. Thus, we conclude that our model may belong to the same universality class of the majority vote model on SHP networks.  All the above numerical analysis is done for square lattices in both vote layer and noise layer. It is found that the Binder cumulant curves do not intersect, which is quite different from the case of random regular graphs. It suggests that there is no clear evidence for the critical phenomenon for the coupled square lattice systems. This result gives us an additional example of the effects of network coupling.

Mean field approach
In the following part, we will present a mean-field method to understand the simulation results. The equation for the magnetization s = á ñ m i is given by [48]: Thus (10) can be written as: For a two-dimensional lattice or a random 4-regular graph, the equation above becomes Next we derive the presentation of (3) under thermodynamic limit (  ¥ N ) and without the diffusion term. For simplicity and without losing any qualitative properties, we set the infinite cluster size = -¥ - q q c 2 for q>q c and = ¥ ¥ C for qq c . Thus (3) can be written as follows: x K x n n is Hill function to constraint the output in [0, 1] and | | m is introduced to depict the additional effect of the magnetization.
The fixed points of (16), (17) can be obtained by solving the following equations: Equation (18a) has a solution for m=0, which corresponds to the disordered state when q>q c . From (18b) the corresponding value of q can be easily solved: 1 with r c =1. Notice that the solution for m can be written as: where the critical exponent β=1 for the mean-field theory. The related solution for q can also be solved. From the analysis above, we could find that the parameters in Hill function only have effects on the final steady state of noise q in disordered state. In figure 8, we show the analytical results of magnetization m and noise parameter q as a function of r given K=10 and n=1 in Hill function. For r<1, the noise parameter q<q c =1/6 in noise layer and there exits an ordered state in vote layer. For r>1, the noise parameter q increases with r from 1/6 to 1/2 and the magnetization in vote layer remains disordered. When r=1, the system will evolve to critical point with m=0 and q=q c =1/6 from any given initial states. Comparing figure 8 with figure 2, we find that the theoretical results are consistent well with the simulation results qualitatively, especially for the case of random 4-regular graphs.

Conclusions
In summary, we study the majority vote model coupled with a reaction-diffusion process on a two-layer multiplex network, namely the vote layer and noise layer. In the context of social dynamics, we assumed that social tension is heterogeneous and dynamical. Simple mechanisms similar to the self-organized sandpile model have been introduced: the self-increasing of social tension and a decreasing with the formation with clusters. Different from the traditional isotropic majority vote model, the noise parameter q i in noise layer is itself a function of time which depends on the network topology and the state of the corresponding node in vote layer, and vice versa. Real world systems are usually open and coevolve with the change of the environment. The framework presented in this article may be helpful for the study of similar problems, such as games or other human behavior in social systems and population dynamics in biological systems.
We perform Monte Carlo simulations to investigate the statistical behavior of the model and an orderdisorder phase transition occurs with the increase of model parameter r=b/a. We also observe that there exists a switching from one sign of the spontaneous magnetization to the opposite near critical region. Notice that r=b/a in our model reflects the intensity ratio of the two mechanisms, it is interesting to see that given a specific r (r=2.09 for simulation model on random 4-regular graphs and r = 1 for mean-field model) in noise layer, the system in vote layer may self-organize to the critical point from any given initial states. To confirm the phase transition, we applied a finite-size scaling analysis to study the critical behavior in vote layer. The obtained critical exponents are similar to those on random graphs and scale-free networks and it may belong to the same universality class of the majority vote model on SHP networks. We also stressed that the continuous phase transition is observed only for random regular graphs in this paper. Our results on square lattices show no evidence for the critical phenomenon, which underlines the importance of certain topology on multilayer systems. It is interesting and valuable for further detailed investigation.
To theoretically analyze our model, a mean-field method is also developed to study the properties of our model. The critical parameter is found to be r c =1 and the phase diagram is in good agreement with our numerical results qualitatively obtained on random regular graphs.
Our model reveals the coupling nature of opinion formation in which usually two or more dynamics are intertwining. By a multilayer network approach, the coupled dynamics of social tension and opinion formation coevolve and a steady state is reached. The intensity ratio of the two mechanisms influences the final coevolution outcome. In future studies, more detailed studies about topologies and different mechanisms are necessary. The following two points are interesting for further detailed investigations. One is the reason for the absence of the continuous phase transition in coupled square lattices and the universality class of the coevolving multilayer network model. How will it change with the network topologies and coupling dynamics? The other interesting topic is the possible first-order phase transition of the system. The existence of bi-stable state shows the possibility for the first-order phase transition. With certain external filed like that in the Ising model, the system may exhibit hysteresis or other phenomena related to the discontinuous phase transition. The related studies may deepen our understanding of coupling dynamics on multilayer networks.