The noisy voter model on complex networks

We propose a new analytical method to study stochastic, binary-state models on complex networks. Moving beyond the usual mean-field theories, this alternative approach is based on the introduction of an annealed approximation for uncorrelated networks, allowing to deal with the network structure as parametric heterogeneity. As an illustration, we study the noisy voter model, a modification of the original voter model including random changes of state. The proposed method is able to unfold the dependence of the model not only on the mean degree (the mean-field prediction) but also on more complex averages over the degree distribution. In particular, we find that the degree heterogeneity—variance of the underlying degree distribution—has a strong influence on the location of the critical point of a noise-induced, finite-size transition occurring in the model, on the local ordering of the system, and on the functional form of its temporal correlations. Finally, we show how this latter point opens the possibility of inferring the degree heterogeneity of the underlying network by observing only the aggregate behavior of the system as a whole, an issue of interest for systems where only macroscopic, population level variables can be measured.


Introduction
Stochastic, binary-state models have been used to study the emergence of collective phenomena in a wide variety of systems and fields.][8][9][10][11] In general, these systems are considered to be embedded in a network structure, where the nodes are endowed with a binary-state variable -spin up or down-and the links between nodes represent the interactions or relations between them.7][18][19] This has been shown to be the case, for example, for the critical temperature of the Ising model, [20][21][22] for the epidemic threshold in spreading phenomena, [23][24][25][26] and for the mean return and first-passage times in random walks. 27,28 us, the quantification of the effect of the underlying topology on such systems and dynamics is, from a practical point of view, a matter of prime importance.
A paradigmatic example of this kind of models, with applications in the study of non-equilibrium systems in a wide range of fields, is the voter model. 4,29 ased solely on local, pairwise interactions, this model assumes that, in a single event, a randomly chosen node copies the state of one of its neighbors, also chosen at random.In this paper, however, we are going to focus on the noisy voter model, a variant of the original voter model which, apart from pairwise interactions, includes random changes of state.This variant model has been studied by, at least, four mutually independent strands of research, largely unaware of each other and belonging to different fields.Namely, percolation processes in strongly correlated systems, 30 heterogeneous catalytic chemical reactions, 31,32 herding behavior in financial markets, 33 and probability theory. 34While both the first and the last strands of literature are directly inspired by the voter model, explicitly using terms such as "noisy voter model" or similar, contributions in the contexts of catalytic reactions and financial markets do not refer to the voter model, and use terms such as "catalytic reaction", "herding" or "Kirman model" instead.More recently, the inclusion of random events in the voter model has also been used to reproduce some of the statistical regularities observed in real electoral processes. 35or any finite system, the behavior of the noisy voter model is characterized by the competition between two opposing mechanisms, related to two different types of noise.On the one hand, the pairwise interaction mechanism is related to interfacial fluctuations (internal noise) and tends to order the system, driving it towards a homogeneous configuration -all spins in the same state, whether up or down-.Depending on the dimension of the system, this mechanism leads to a coarsening process or to a metastable partially ordered state, both of them perturbed by finite-size fluctuations (one of which eventually drives the system to full order).In the absence of any other mechanism, as it is the case in the original voter model, the homogeneous configurations become absorbing states of the dynamics. 36On the other hand, the random change mechanism is related to thermal-like fluctuations (external noise) and tends to disorder the system, pulling it from the homogeneous configurations.Therefore, this second mechanism leads to the disappearance of the typical absorbing states of the voter model and to the restoration of ergodicity. 34The main consequence of this competition is the appearance of a noise-induced, finite-size transition between two different behavioral regimes -a mostly ordered regime dominated by pairwise interactions and a mostly disordered regime dominated by noise. 33,37 hile the effect of different network topologies on the behavior of the voter model has been well established, [38][39][40][41] the case of the noisy voter model has received much less attention, most of the corresponding literature focusing only on regular lattices 30,34 or on a fully-connected network. 33,37 inally, the use of a mean-field approach in some recent studies considering more complex topologies [42][43][44] did not allow to find any effect of the network properties -apart from its size and mean degree-on the results of the model.
In this paper, we move beyond the usual mean-field approximations 42,44,45 and propose an alternative analytical approach, based on an annealed approximation for uncorrelated networks and inspired by a recently introduced method to deal with heterogeneity in stochastic interacting particle systems. 46In particular, we approximate the network by a complementary, weighted, fully-connected network whose weights are given by the probabilities of the corresponding nodes being connected in uncorrelated networks of the configuration ensemble [47][48][49][50] -i.e., proportional to the product of their degrees.Furthermore, we present a formulation of the problem in terms of a master equation for the probability distribution of the individual states of the nodes.In this way, we are able to find approximate analytical expressions for the critical point of the transition, for a local order parameter and for the temporal correlations.As opposed to previous mean-field approaches, we find that the degree heterogeneity -variance of the underlying degree distribution-has a significant impact on all of these variables, leading to a larger value of the critical point, a higher level of order and a modification of the functional form of the temporal correlations.As we will show, this latter point opens the possibility of inferring the degree heterogeneity of the underlying network by observing only the aggregate behavior of the system as a whole, an issue of interest for systems where only macroscopic, population level variables can be measured.Finally, these results are confirmed by numerical simulations on different types of networks, allowing for a constant mean degree (k = 8) while leading to different degree distributions.In particular, in order of increasing degree heterogeneity, we focus on Erdös-Rényi random networks, 51 Barabási-Albert scale-free networks 52 and dichotomous networks 16 (whose nodes are assigned one out of two possible degrees, in our case

Model definition
Consider a system composed of N nodes in a given network of interactions.At any point in time, each node i is considered to be in one of two possible states, and is therefore characterized by a binary variable s i = {0, 1}.Moreover, due to the network structure, each node i is also characterized by a certain set of (nearest) neighbors, nn(i), and by its corresponding degree or number of those neighbors, k i .The evolution of the state of each node, s i , occurs stochastically with probabilities that depend on the state of the updating node and on the states of its neighbors.In particular, these probabilities consist of two terms: on the one hand, there are random pairwise interactions between node i and one of its neighbors j ∈ nn(i), after which i copies the state of j; and, on the other hand, there are random changes of state, playing the role of a noise.The transition rates for each node i can be written as where the noise parameter a regulates the rate at which random changes of state take place, and the interaction parameter h does so with the interaction-driven changes of state.Defined in this way, the noisy voter model becomes the network-embedded equivalent of the Kirman model in its original, extensive formulation. 33,37,42 Frthermore, note that, in the limit case of a = 0 and with an appropriate time rescaling, we recover the transition rates of the original voter model.Let us stress here that, even if the model appears to have two parameters, one of them can always be used as a rescaling of the time variable, so that there is only one relevant parameter: the ratio between the two introduced coefficients, a/h.Indeed, only one parameter is introduced in the previous literature in the context of the noisy voter model. 30,34,44 O the contrary, prior works dealing with the Kirman model usually keep both parameters. 33,37,42 Fr consistency with one and the other strands of literature, we are going to consider both parameters explicitly in our analytical approach, while we keep the interaction parameter fixed as h = 1 for our numerical results -allowing the noise parameter a to vary.
In order to characterize the global state of the system we introduce the global variable n, defined as the total number of and taking values n ∈ 0, 1, ..., N. Note that this variable does not take into account any aspect of the network structure.It should be observed that, for a = 0, there are no absorbing states in the model -the probability to move from one state to any other is strictly positive-and therefore the Markov chain is said to be ergodic: in the steady state, averages over time are equivalent to ensemble averages.In practice, the smaller a is, the longer the time needed for both statistics to be actually equivalent.Thus, in the limit case of a = 0 the time needed becomes infinite, and we recover the voter model behavior: non-ergodicity with two absorbing states, at n = 0 and n = N.Moreover, we are going to use the notation x for ensemble averages with random initial conditions, while we leave f (k) for averages over the degree distribution, i.e., Similarly, we will differentiate between the variance of a variable x over realizations, noted as σ 2 [x], and the variance of the degree distribution, labeled as σ 2 k .Note, nonetheless, that for the numerical steady state values to be presented in the following sections, averages are performed both over time and over an ensemble of realizations with random initial conditions, assuming an initial transient of N time units.

General formulation
The stochastic evolution of the system can be formalized as a Markov process.In particular, we can write a general master equation for the N-node probability distribution P(s 1 , . . ., s N ) (see Supplementary Information) and use it to derive general equations for the time evolution of the first-order moments and the second-order cross-moments of the individual nodes' state variables s i , where q i j = r + i + r − i + r + j + r − j and δ stands for the Kronecker delta (see Appendices B and C for details).In general, if the transition rates depend on the individual state variables s i , these equations involve higher order moments and they cannot be solved without a suitable approximation. 46However, for the transition rates of the noisy voter model, due to their particular form, both equations become independent of higher order moments.
For the first-order moments, introducing the transition rates (1) into equation ( 4), we obtain an equation directly solvable in the steady state, when the influence of the initial conditions has completely vanished and thus s i st is independent of i.In this way, we find, for the steady state average individual variables s i and, by definition, for the steady state average global variable n, respectively, the expected results given the symmetry of the system.In the case of the second-order cross-moments, when we introduce the transition rates (1) into equation ( 5), we obtain which, even if independent of higher order moments, cannot be solved in the absence of an explicit knowledge of the network connections -the adjacency matrix.This is due to the presence of sums over neighbors ∑ m∈nn(i) where the terms are not independent of the particular pair of nodes m, i.In order to find the corresponding steady state solution, we introduce below an approximation of the network allowing us to write the previous equation in terms of sums over the whole system.

Annealed approximation for uncorrelated networks
Given a complex network with adjacency matrix A i j and degree sequence {k i }, we can use an annealed graph approach [53][54][55] to define a complementary, weighted, fully-connected network with a new adjacency matrix Ãi j and whose structural properties resemble those of the initial network. 50In particular, we assume that the weights of this new adjacency matrix are given by the probabilities of the corresponding nodes being connected, that is, Ãi j = p i j , where p i j is the probability of node i, with degree k i , being connected to node j, with degree k j .For uncorrelated networks of the configuration ensemble, i.e., random networks with a given degree sequence {k i } and with a structural cutoff at k i < √ Nk, we can approximate the probability of two nodes i, j being connected [47][48][49]56 by In this way, we can approximate the sums over the neighbors of a given node i as sums over the whole network, where f j is a function which can depend on the characteristics of node j (k j and/or s j ).Note that this approximation preserves the initial degree sequence, as it is obvious from and, therefore, the total number of links is also conserved.

Noise-induced, finite-size transition
As shown in the previous literature about the Kirman model, 33,37 in the fully-connected case, the system is characterized by the existence of a finite-size transition between a bimodal and a unimodal behavior, depending on the relative magnitude of the noise and the interaction parameters.For a < h/N the steady state probability distribution of n is found to be bimodal with maxima at the extremes or fully ordered configurations, n = 0 and n = N, meaning that, at any point in time, the most likely outcome of a static observation is to find a large majority of nodes in the same state, whether 0 or 1, with different observations leading to different predominant options (see 57 for an explanation in terms of an effective potential).On the contrary, for a > h/N the distribution of n becomes unimodal with a peak at n = N/2, meaning that, at any point in time, the most likely outcome of an observation is to find the system equally split between both options.Given the ergodicity of the model for a = 0, these probability distributions can also be understood in terms of the fractional time spent by the system with each value of n.In this manner, in the bimodal regime, stochastic realizations of the process will tend to be temporarily absorbed in the proximity of the fully ordered configurations with random switches between them, while realizations in the unimodal regime will spend most of the time with the system more or less equally divided among the two possible individual states, 0 and 1.At the critical point marking the transition between these two behaviors, a c = h/N, the distribution of n becomes uniform, meaning that any share of nodes between the two options is equally likely.Note that this transition is a finite-size effect, since the value of the critical point decreases for increasing system size and vanishes in the thermodynamic limit (N → ∞).
The existence of the referred transition when the system is embedded in a network topology has also been reported in the literature both for the Kirman model 42,43 and in the context of the noisy voter model. 44The above described phenomenology can thus also be observed in different network topologies.As an example, we show in Fig. 1 two realizations of the dynamics for a Barabási-Albert scale-free network corresponding, respectively, to the bimodal [panel a)] and the unimodal regime [panel b)].A mean-field approach has been proposed in the literature, 42,44 leading to an analytical solution for the critical point which does not depend on any property of the network other than its size, a c = h/N, the transition still being a finite-size effect.
Both the analytical and numerical results to be presented here suggest, on the contrary, that the critical point does depend on the network, while they confirm the finite-size character of the transition.As a quantitative description of the transition we are going to use the variance of n: bearing in mind that the variance of a discrete uniform distribution between 0 and N is N(N + 2)/12, we can identify the critical point of the transition as the relationship between the model parameters which leads the steady state variance of n to take the value σ 2 st [n] = N(N + 2)/12.Although it is not necessarily the case, numerical results confirm that the distributions obtained in this manner are indeed uniform.The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.

Variance of n
Introducing the annealed approximation for uncorrelated networks into the equation for the second-order cross-moments of the individual variables s i , equation (8), we can replace the sums over sets of neighbors by sums over the whole system.If we then rewrite this equation in terms of the covariance matrix σ i j , defined as we can use the relation to find an equation for the variance of n, by simply summing over i and j.Finally, after some algebra (see Appendix D for details), we find, in the steady state, under the necessary and sufficient condition that which is generally true and always true for h > 0 and k ≥ 2. Note that the only approximation used in the derivation of equation ( 14) is the estimation of the adjacency matrix involved in the annealed approximation for uncorrelated networks.

5/28
The behavior of the variance σ 2 st [n] as a function of the noise parameter a is shown in Fig. 2 for the three types of networks studied.As we can observe, despite a small but systematic overestimation for intermediate values of the noise parameter -attributable only to the annealed approximation for uncorrelated networks, the only one involved in its derivation-, the main features of the numerical steady state variance are correctly captured by the analytical expression in equation (14).In particular, both its dependence on a and the impact of the underlying network structure are well described by our approach.On the contrary, the mean-field solution proposed in the previous literature, 42 and included in Fig. 2 for comparison, fails to reproduce the behavior of the variance of n for large a and is, by definition, unable to explain its dependence on the network topology.It is, nonetheless, a good approximation for the Erdös-Rényi random network and for values of the noise parameter a 10 −1 .  .Dashed line: Mean-field approximation (see 42 ).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.
Regarding the limiting behavior of the system when a → 0 and when a → ∞, we can observe, for both the numerical and the analytical results presented in Fig. 2, that the influence of the network on the steady state variance of n vanishes in both limits, where we recover the expected behaviors.Notably, in the limit of a → 0 the variance tends to N 2 /4 for all networks, and we progressively recover the voter model behavior; while in the limit of a → ∞ the variance tends to N/4 regardless of the topology, as it corresponds to a purely noisy system composed by N independent units adopting, randomly, values 0 or 1 (equivalent, as well, to a one dimensional random walk confined to the segment [0, N]).
Concerning the impact of the network structure, we can observe in Fig. 2 that for any finite value of the noise parameter, 0 < a < ∞, a larger degree heterogeneity of the underlying topology, measured as the variance of the corresponding degree distribution, leads to a larger steady state variance of n.This behavior is further confirmed by the results to be presented in the next subsection, where we show the steady state variance of n as a function of the variance of the underlying degree distribution σ 2 k , respectively, for two different values of the noise parameter a.As we can observe, even if the numerical results are systematically overestimated, our analytical approach [equation (14)] is able to capture the general features of this dependence and represents a significant improvement from the mean-field prediction of no network impact.
In order to study the bimodal-unimodal transition by using the behavior of the variance of n illustrated in Fig. 2, the variance value corresponding to a uniform distribution is included as a horizontal line, so that the critical a value for each network can be easily identified at the corresponding intersection (marked by vertical dashed lines).Note that values of the variance of n above (below) the uniform distribution line correspond to the system being in the bimodal (unimodal) phase.A first observation is that the referred transition still occurs when the noisy voter model is embedded in a network topology, thus confirming the results reported in the previous literature. 42,44 owever, as opposed to these previous studies, we can observe in Fig. 2 a clear dependence of the critical point on the underlying topology, an effect which seems to be correctly captured by our approach while it goes completely unnoticed, by definition, from a mean-field perspective.The particular features of this dependence will become clear by means of a first-order approximation of the steady state variance σ 2 st [n] with respect to the system size N, allowing us to characterize the asymptotic behavior of the system for both small and large a as well as to find an explicit expression for the critical point a c .

Asymptotic behavior of the variance of n
Given that equation (14) does not allow for an intuitive analytical understanding of the network influence on the steady state variance of n, nor does it allow for an explicit analytical solution for the critical point a c , we develop here a first-order approximation with respect to the system size N, which will also give a relevant insight regarding the asymptotic behavior of the system for both small and large a.In fact, the result of this approximation strongly depends on the relationship between the system size N and the noise parameter a, and we are thus led to consider two different approximation regimes.
In particular, when the noise parameter a is of order O(N −1 ) or smaller, then the product aN is, at most, of order O(N 0 ), and a first-order approximation of equation ( 14) with respect to the system size N leads to corresponding to the asymptotic behavior of the variance of n for small a and large N. On the contrary, when a is of order O(N 0 ) or larger, the product aN is, at least, of order O(N), and the first-order approximation of equation ( 14) becomes corresponding to the asymptotic behavior of the variance of n for large a and large N (see Appendix E for details).
For a more precise characterization of the ranges of validity of these two asymptotic approximations with respect to the noise parameter a, we present in Fig. 3 the variance of n as a function of a for the numerical results and the three corresponding analytical expressions presented so far: the analytical result in equation ( 14), the asymptotic expression for small a in equation ( 16) and the asymptotic expression for large a in equation (17).Note the use a Barabási-Albert scale-free network as an example.Furthermore, we also show in this figure the crossover point a * between both approximations, that we define as the value of a that minimizes the distance between the logarithmic values of both functions (17) and (16).
Noticing that, for both asymptotic approximations, the variance σ 2 st [n] becomes an explicit function of the variance of the underlying degree distribution σ 2 k , we present in Fig. 4 a comparison between these analytical functional relationships and the corresponding numerical results for two different values of the noise parameter a.In particular, taking into account the ranges of validity of the asymptotic approximations characterized above (see Fig. 3), we chose values of the noise parameter respectively before [panel a)] and after [panel b)] the crossover point a * , and both of them in the region of a leading to significant differences between network types (see Fig. 2).
As we can observe, each asymptotic approximation accurately fits the analytical result in equation ( 14) within its respective range of validity, while it becomes clearly inaccurate out of this range.Therefore, we can use these approximations instead of equation ( 14) to better understand the behavior of the system.In this way, we can conclude that, regarding its impact on the results of the model, the most relevant property of the underlying network is not its mean degree, but the variance of its degree distribution relative to the square of its mean degree, σ 2 k /k 2 , a normalized measure of its degree heterogeneity.The results presented in Fig. 4 show that this analysis significantly outperforms the mean-field prediction of no network impact, particularly for networks with large levels of degree heterogeneity.Note, nonetheless, that both asymptotic approximations are subject to the same inaccuracies in reproducing the numerical results as the original analytical expression, i.e., the inaccuracies caused by the annealed approximation for uncorrelated networks: a systematic overestimation of the numerical results and an inability to explain the results for topologies with large structural correlations.

Critical point
As described above, the critical point of the bimodal-unimodal transition can be defined as the relationship between the model parameters a and h leading the steady state variance of n to take the value σ Analytical results [see equation (14)].Dotted line: asymptotic approximation for small a [see equation (16)].Dash-dotted line: asymptotic approximation for large a [see equation (17)].Dashed line: Crossover point between both asymptotic approximations (a * = 0.014157).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.
uniform distribution between 0 and N. A numerical solution for the critical point a c can thus be found by applying this definition to the analytical expression for the variance of n given in equation (14).However, for a fully analytical description of the critical point, we have to use one of the asymptotic approximations presented above, algebraically solvable for a c .In particular, bearing in mind that the value of the critical point of a fully-connected system is of order O(N −1 ) and that the change due to the network structure appears to be of order O(N 0 ) (see Fig. 2), then we can expect the value of the critical point to be still of order O(N −1 ) and we can therefore use the small a asymptotic approximation in equation ( 16) to find to the first-order in N (see Appendix F for details).Both this expression and the mean-field approximation previously proposed in the literature 42,44 are contrasted with numerical results in Fig. 5, where we present the values of the critical point a c for different types of networks as a function of the variance of the corresponding degree distributions σ 2 k .As before, we notice in Fig. 5 a systematic overestimation of the numerical results by our analytical approach, whose origin lies, again, in the annealed approximation for uncorrelated networks.While both equation (18) and the mean-field approximation are able to capture the finite-size character of the transition -the fact that a c → 0 when N → ∞-, only our approach is able to reproduce the influence of the underlying network structure on the critical point.In particular, we observe a numerical behavior approximately consistent with a linear relationship between the value of the critical point and the variance of the underlying degree distribution, as predicted by equation (18).A quantitative assessment of the significance of this dependence can be obtained by observing the shift between the critical points corresponding to the Erdös-Rényi random network and the dichotomous network, the latter being almost a factor of 3 larger than the former.While the persistence of the bimodal-unimodal, finite-size transition in different network topologies had already been reported, 42,44 to the best of our knowledge, no dependence of the critical point on the characteristics of the underlying network has been documented so far for the noisy voter model nor in the context of the Kirman model (see 16 for a similar effect in a different model).k for two values of the noise parameter a.In order to keep all parameters constant except the variance of the degree distribution, a different network type is used for each point (in order of increasing σ 2 k : Erdös-Rényi random network, Barabási-Albert scale-free network and dichotomous network).Circles with error bars: Numerical results (averages over 20 networks, 10 realizations per network and 50000 time steps per realization).Solid line and squares: Analytical results [see equation (14)].Dotted line: asymptotic approximation for small a [see equation (16)].Dash-dotted line: asymptotic approximation for large a [see equation (17)].Dashed line: Mean-field approximation (see 42 ).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.

Local order
We can characterize the local order of the system with an order parameter ρ defined as the interface density or density of active links, that is, the fraction of links connecting nodes in different states, where A i j are the elements of the adjacency matrix.Larger values of ρ imply a larger disorder, corresponding ρ = 1/2 to a random distribution of states, while ρ = 0 corresponds to full order.Furthermore, note that, as opposed to n, the order parameter does take into account the structure of connections between nodes.
While it has not been studied before in the context of the Kirman model, the interface density ρ is commonly used to describe the time evolution of the voter model. 40In the absence of noise, the voter model is characterized by the existence of two absorbing states (n = 0 and n = N), both of them corresponding to full order (ρ = 0).Therefore, the focus is on how the In order to keep all parameters constant except the variance of the degree distribution, a different network type is used for each point (in order of increasing σ 2 k : Erdös-Rényi random network, Barabási-Albert scale-free network and dichotomous network).Symbols: Numerical results (averages over 20 networks, 10 realizations per network and 50000 time steps per realization).Solid line: Analytical results [see equation (18)].Dashed line: Mean-field approximation (see 42,44 ).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8. system approaches these absorbing ordered states.In the presence of noise, on the contrary, the system has no absorbing states, i.e., it is always active.Thus, the focus is not anymore on how it reaches any final configuration, but rather on characterizing its behavior once the influence of the initial condition has vanished, that is, in the steady state.In the context of the noisy voter model, it has been recently shown that, after a short initial transient, the average interface density reaches a plateau at a certain value ρ st , with ρ st > 0 for any non-zero value of the noise and ρ st = 1/2 in the infinite noise limit. 44Moreover, a mean-field pair-approximation has been used to find an analytical solution for ρ st as a function of the level of noise and the mean degree of the underlying network.This analytical solution has been shown to be a good approximation for large values of the noise parameter, the small noise region not having been considered.
Let us start the description of our results by emphasizing that individual realizations of the interface density ρ remain always active for any non-zero value of the noise, as it was also the case for the variable n (see Fig. 1).As an example, we show in Fig. 6 two realizations of the dynamics for a Barabási-Albert scale-free network corresponding, respectively, to the bimodal [panel a)] and the unimodal regime [panel b)].While in the first of them (a < a c ) the system fluctuates near full order, with sporadic excursions of different duration and amplitude towards disorder; in the second (a > a c ), the system fluctuates around a high level of disorder, with some large excursions towards full order.
Introducing the annealed approximation for uncorrelated networks described above into the definition of the order parameter given in equation (19), and focusing on the steady state average value, we obtain In this way, an explicit solution for the steady state average interface density can be found by expressing it in terms of the analytical results presented so far, namely, in terms of the variance σ 2 st [n] (see Appendix G for details), This expression can be contrasted with numerical results in Fig. 7, where we present the steady state average interface density  1).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.
ρ st as a function of the noise parameter a for different types of networks.The mean-field pair-approximation result derived in 44 is also included for comparison.
As we can observe in Fig. 7, our approach correctly captures the behavior of the system for both small (a 10 −3 ) and very large values (a 3) of the noise parameter: both the asymptotic convergence towards ρ st = 0 for small a (voter model result for finite systems) and the convergence towards ρ st = 1/2 for large a (full disorder) are well reproduced.On the contrary, our analytical approach fails to reproduce the numerical results for intermediate values of the noise parameter (10 −3 a 3).The origin of this discrepancy lies in the annealed network approximation: when replacing the original network by a weighted fully-connected topology, all track of local effects is lost -precisely those measured by the order parameter.The fact that this discrepancy is only present for intermediate values of a can be explained, on the one hand, by the lack of locally ordered structures in the fully disordered, large a regime and, on the other hand, by the development of a global order -more and more independent of local effects-for decreasing values of a. Thus, an accurate fit of the numerical results for any value of a can only be expected for topologies where local effects are absent or negligible.In the Supplementary Figure S1 presented in Appendix I we show that the approximation successfully fits the results in a fully-connected network.The good accuracy of the results presented above for the variance of n suggests that the discrepancy between analytical and numerical results appears only when the annealed network approximation is used to derive a relationship between ρ st and σ 2 st [n], and not in the derivation of the latter, for which only global correlations are relevant.Apart from the functional dependence of the interface density on the noise parameter, our approach is also able to capture the influence of the network, which becomes significant for a 10 −2 .In particular, we find that a larger variance of the degree distribution of the corresponding network leads to a smaller interface density, i. e., to a higher level of order.
Even if the mean-field pair-approximation fits the numerical results remarkably well for large and intermediate values of the noise parameter (a 10 −2 ), it is completely unable to reproduce the behavior of the system for small a, and it fails to explain the influence of any network property other than the mean degree.While the pair-approximation allows to capture the short-range order characteristic of intermediate values of a, the assumptions implicit in the derivation of the mean-field result 45 do not allow to reproduce the long-range order characteristic of the small a region.Note that the limiting case of a = 0 (voter   21)].Dashed line: Mean-field pair-approximation (see 44 ).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8. model) is a singular point of the mean-field pair-approximation, 44 leading to the existence of two different solutions: a non-zero solution linked to the result displayed in Fig. 7 -correct in the infinite size limit-, and a zero solution -correct for finite systems.

Inference of network properties from the autocorrelation of n
As explained above, from any initial condition, the system quickly reaches a dynamic steady state, whose active character can be clearly observed in Fig. 1.In order to characterize the dynamic nature of this steady state, let us now focus on the steady state autocorrelation function of n, defined as where τ plays the role of a time-lag.In the fully-connected case, it has been shown in the previous literature 37 that the autocorrelation decays exponentially, with an exponent proportional to the noise parameter, In the case of different network topologies, the mean-field prediction is that no influence of the network is to be expected and, therefore, the same exponential decay as in the fully-connected case is to be found.In contrast with this prediction, both the analytical and numerical results to be presented here show that the network does have a significant impact on the functional form of the steady state autocorrelation of n.
Introducing the annealed approximation for uncorrelated networks described above into the equation for the time evolution of the first-order moments (6), integrating it with carefully chosen initial conditions, and making use of the above reported analytical results (see Appendix H for details), we can find where S 1 is defined as 12/28 This expression can be contrasted with numerical results in Fig. 8, where we present the autocorrelation function, normalized by the variance, for the two extreme cases of a network with no degree heterogeneity (regular 2D lattice) and a highly heterogeneous degree distribution (dichotomous network).Note the logarithmic scale in the y-axis.
Dichotomous analytical Dichotomous numerical Regular 2D analytical Regular 2D numerical It is important to note that, in the case of no degree heterogeneity, the new variable S 1 becomes S 1 = σ 2 st [n].This can be understood by applying, ∀i, k i = k in the averages over the degree distribution in equation ( 14), and introducing this result into the definition of S 1 , equation (24).Thus, for networks with no degree heterogeneity, the steady state autocorrelation function behaves as in the fully-connected case, and as predicted by the mean-field approximation, This single exponential decay is confirmed by the numerical results presented in Fig. 8 for the regular 2D lattice.
On the contrary, for networks with non-zero degree heterogeneity, in general, S 1 = σ 2 st [n], and thus the autocorrelation function consists of two different exponential decay components [see equation (23)].When h > 0, the exponential e −(2a+h)τ decays faster than e −2aτ .Therefore, for long time-lags, we expect the normalized autocorrelation function of any network to be parallel to e −2aτ in log-linear scale, with a vertical shift proportional to its degree heterogeneity and due to the initial deviation from the single exponential behavior.This description is confirmed by the numerical results presented in Fig. 8 for the dichotomous network.
Neither equation ( 14) nor the asymptotic approximate expressions ( 16) and ( 17) allow to infer, for a given system, the values of the two model parameters, a and h, and the normalized variance of the underlying degree distribution, σ 2 k /k 2 , by measuring only the steady state variance of n, σ 2 st [n].Thus, it is impossible, by using only these relationships, to conclude if the fluctuations observed in a given system have a contribution due to the degree heterogeneity of the network, without a prior knowledge of the model parameters a and h.On the contrary, the particular functional form of the autocorrelation function K st [n](τ) -with two exponential decay components whose exponents are different functions of a and hdoes allow for the values of a, h and σ 2 k /k 2 to be inferred from equation (23), in combination with equation (16) or equation (17), by measuring only the temporal correlations of the aggregated variable n, and assuming we can also know the system size N.Note that the use of equation (16) or equation ( 17) can be determined by self-consistency, depending on the value obtained for a.As an example, a simple fit of equation ( 23) to the numerical results presented in Fig. 8 for the dichotomous network leads, in combination with equation ( 16), to the fitted parameter values a = 0.0099, h = 0.94 and σ 2 k /k 2 = 2.539, remarkably close to the actual values used for computing the numerical results, a = 0.01, h = 1 and σ 2 k /k 2 = 2.625.In this way, we are able to infer some information about the underlying network -its normalized level of degree heterogeneity-by studying only the aggregate behavior of the system as a whole.

Discussion
In this paper, we have proposed a new analytical method to study stochastic, binary-state models of interacting units on complex networks.Moving beyond the usual mean-field theories, 42,44,45 this alternative approach builds on a recent study considering heterogeneity in stochastic interacting particle systems 46 and proposes an annealed approximation for uncorrelated networks accounting for the network structure as parametric heterogeneity.
Using the noisy voter model as an example, we have been able to unfold the dependence of the model not only on the mean degree of the underlying topology (the mean-field prediction) but also on more complex averages over the degree distribution.
In particular, we have shown that the degree heterogeneity -i.e., the variance of the underlying degree distribution-has a substantial influence on the location of the critical point of the noise-induced, finite-size transition characterizing the model.This shift of the transition might have important practical implications in real systems, since it suggests that different behavioral regimes can be achieved by introducing changes in the underlying network of interactions.Furthermore, we have studied the influence of the network on the local ordering of the system, finding that a larger degree heterogeneity leads to a higher average level of order in the steady state.Interestingly, we have also found the heterogeneity of the underlying degree distribution to play a relevant role in determining the functional form of the temporal correlations of the system.Finally, we have shown how this latter effect can be used to infer some information about the underlying network -its normalized level of degree heterogeneity-by studying only the aggregate behavior of the system as a whole, an issue of interest for systems where macroscopic, population level variables are easier to measure than their microscopic, individual level counterparts.
Numerical simulations on different types of networks have been used to validate our analytical results, finding a remarkably good agreement for all the properties studied except for the local order, for which a significant discrepancy is found for intermediate levels of noise.The origin of this discrepancy has been shown to lie in the annealed network approximation, whose validity is restricted to global properties or situations where local effects are negligible.The generally good agreement found is all the more remarkable considering that, while the uncorrelated network assumption is essential for the proposed analytical method, we did not impose any particular structural constraint to avoid correlations in the networks used for the numerical simulations. 48,58

A Master equation
We derive here a general master equation for the N-node probability distribution P(s 1 , . . ., s N ), where the individual node variables are binary and take the values s i = {0, 1}.Recalling that r + i is the rate at which node i changes its state from s i = 0 to s i = 1 and r − i the rate at which it does so in the opposite direction, we can directly write differential equations for the probability of node i to be in state s i = 0 and for its probability to be in state s i = 1, respectively, Introducing here the individual-node step operators E +1 i and E −1 i , whose effect over an arbitrary function of the state of node i, f (s i ), is defined as we can rewrite equations (A1) as Multiplying these two equations, respectively, by (1 − s i ) and s i , we can gather them in a single differential equation, and noticing that (1 we can rearrange terms as Finally, we find the master equation for the N-node probability distribution P(s 1 , . . ., s N ) by simply adding up the contribution of every single node

B Equation for the time evolution of the first-order moments s i
We show, in this section, how to obtain a general equation for the time evolution of the first-order moments s i [equation ( 4) in the main text].Let us start by using the definition of the step operators in equation (A2) and the binary character of each individual node state variable, s i = {0, 1}, to derive, for a given function of the state of node i, f (s i ), four relations which will ease later calculations.While the function f might also depend on the other variables, f = f (s 1 , ..., s i , ..., s N ), we restrict our attention, without loss of generality, to the case f (s i ).For the first two relations, we have that where the sums are over the two possible values of s i .Looking at the master equation (A6), one can understand that these two relations translate the fact that any increase in the probability of a given node being in a given state must be accompanied by a corresponding decrease in the probability of the complementary state.Regarding the other two relations, we can write (B4)

15/28
Let us also introduce, for clarity, the notation ∑ {s} to refer to the sum over all the possible combinations of states of all the individual nodes' variables, and ∑ {s} j to indicate the sum over all the possible combinations of states of all the variables except s j , Note that these two definitions are related by which allows us to split the sum over all possible configurations of the system into a sum over the values of one of the variables and a sum over the configurations of the rest of the system.By using the notation in (B5), the average of a given function of the states of the nodes, f (s 1 , . . ., s N ), can be written as Using this expression and the master equation in (A6) we derive an equation for the time evolution of the average value of the state of node i, d s i dt Separating the terms with j = i and those with j = i, we find If we now use the relation (B7) to extract, from the general sum over {s}, the sum over the values of s i for the terms with j = i, while we extract the sum over the values of s j for the terms with j = i, we obtain where we can easily identify relations (B1) and (B2) for the terms with j = i, and relations (B3) and (B4) for the terms with j = i.In this way, we can write which, after combining the sums together again, becomes and we finally find the equation for the time evolution of the first-order moments presented in the main text, C Equation for the time evolution of the second-order cross-moments s i s j In order to find a general equation for the time evolution of the second-order cross-moments s i s j [equation (5) in the main text] we proceed in a similar way as we did in the previous section for the first-order moments.Taking into account the master equation (A6) and using the definition of the average value in (B8), we can write for the second-order cross-moments, For the terms of the sum with k = i, j, we can use relation (B7) to write where the equality follows from an application of relations (B1) and (B2).Similarly, we can use relations (B3) and (B4) to transform, in equation (C1), the terms with k = i = j as and, equivalently, the terms with k = j = i as Note that, for both expressions (C3) and (C4), we have assumed that i = j.In order to study the other case, when i = j, we simply need to notice that, being the possible values of the variables s i = {0, 1}, then s 2 i = s i , and therefore where we have used the result (B14) for the first-order moments derived in the previous section.Thus, we can write an equation for the second-order cross-moments as where q i j = r + i + r − i + r + j + r − j .Finally, using the Kronecker delta, we obtain the expression presented in the main text,

D Variance of n
We derive here an analytical expression for the steady state variance of n [equation (14) in the main text].Let us start by introducing the transition rates of the noisy voter model [equation (1) in the main text] into the equation for the time evolution of the second-order cross-moments obtained in the previous section, equation (C7), s m s i − 2(2a + h) s i s j s m s i . (D1) Applying now the annealed approximation for uncorrelated networks described in the main text [see equation (10)], we can replace the sums over sets of neighbors by sums over the whole system, finding Bearing in mind the definition of the covariance matrix, σ i j = s i s j − s i s j , we can find an equation for its time evolution from equation (6) in the main text and equation (D2), which can be written in terms of only the covariance matrix and the first moments, In the steady state, and using also the steady state solution of the first order moments s i st = 1/2 [equation (7) in the main text], we find Note that, for the sake of notational simplicity, we have dropped the subindex st for the steady state solution of the covariance matrix.Recalling now the relation between the variance of n and the covariance matrix [equation (13) in the main text], we can find an equation for the steady state variance of n by simply summing equation (D5) over i and j, Let us introduce now the set of variables S x , with x ∈ {0, 1, 2, . ..}, and defined as In this way, we can rewrite the steady state variance of n in terms of one of these new variables, S 0 , In order to find an equation for this new variable S 0 , we could use again the equation for the covariance matrix in (D5), multiplying it by k j and summing over i and j, obtaining a solution in terms of the variable S 1 .We could then proceed similarly and find an equation for S 1 as a function of S 2 , for S 3 as a function of S 4 , and so forth.In general, for any x, we have where the overbar notation is used for averages over the degree distribution [see equation (3) in the main text].From equation (D9) we can obtain an expression for the variable S x in terms of only S 1 and S x+1 , By inverting equation (D10), we can write all variables S x+1 in terms of the preceding ones, which has the general form It is easy to see that this recurrence relation has the solution where the choice of S 1 instead of S 0 in the first term allows us to write all the variables S x+1 in terms of only one of them, S 1 .
Note that this choice is required by the presence of a term with S 1 inside B x .Thus, we can write the solution for our original recurrence relation in (D11) as If we now rewrite equation (D14) as we find that the left hand side of this equation vanishes in the limit of x → ∞, where we have used the definition of the variables S x given in equation (D7).A necessary and sufficient condition for the last equality in equation (D16) to hold is that which is generally true and always true for h > 0 and k ≥ 2. Thus, in the x → ∞ limit, we can equate the right hand side of equation (D15) to zero, and find, in this way, a solution for S 1 , Regarding the sums in equation (D19), we can use the sum of the geometric series where the condition of convergence is exactly the same as presented before in equation (D17), and thus generally true and always true for h > 0 and k ≥ 2. In this way, applying the result (D20) to equation (D19) we have where the denominator can be rewritten as thereby finding a final expression for S 1 , If we now go back to the equation for the steady state variance σ 2 st [n] as a function of S 0 , equation (D8), and we use equation (D10) to find an expression for S 0 as a function of S 1 , then we can write an equation for the steady state variance as a function of S 1 , Finally, introducing here what we found for S 1 in equation (D23), we arrive to the final expression for the steady state variance of the global variable n as presented in the main text,

E Asymptotic approximations for the variance of n
We develop here a first-order approximation for the steady state variance of n with respect to the system size N.Given the dependence of the result of this approximation on the relationship between the system size N and the noise parameter a, we are forced to consider two different asymptotic approximation regimes: one for small a [corresponding to equation ( 16) in the main text] and the other for large a [corresponding to equation (17) in the main text].
Let us start by noticing that the structural constraint imposed by the annealed approximation for uncorrelated networks on the degrees of the network, k i < √ Nk, allows us to write equation (D26) as

21/28
In this way, we notice that, depending on the order of the product aN, the approximation of the third term in equation (E1) will lead to different results.In particular, when the noise parameter a is of order O(N −1 ) or smaller, then the product aN is, at most, of order O(N 0 ), and we can continue with the approximation as which, to the first order in N, becomes Using now the definition of the variance of the degree distribution, σ 2 k = k 2 − k 2 , we find the approximation presented in the main text for the steady state variance of n for small a and to the first order in N, Note that the remaining terms are at most of order O(N 3/2 ).
On the contrary, when a is of order O(N 0 ) or larger, then the product aN is, at least, of order O(N), and we can approximate equation (E1) as order in N we have and, finally, we find the approximation presented in the main text for the steady state variance of n for large a and to the first order in N, where the remaining terms are at most of order O(N 1/2 ).

F Critical point approximation
In this section, we derive an analytical approximation for the critical point of the bimodal-unimodal transition [equation (18) in the main text], which can be defined as the relationship between the model parameters a and h leading the steady state variance of n to take the value σ 2 st [n] = N(N + 2)/12, corresponding to a uniform distribution between 0 and N.In particular, bearing in mind that the critical value a c of a fully-connected system is of order O(N −1 ) and that the change due to the network structure appears to be of order O(N 0 ) (see Fig. 2 in the main text), then we can expect the value of the critical point to be still of order O(N −1 ), and we can therefore use the small a asymptotic approximation in equation (E4), The solution of this equation leads to the, for large N, leads to the value of the critical point discussed in the main text, consistent with the assumption of a critical value of order O(N −1 ).Note that assuming, instead, the critical value to be of order O(N 0 ), and using therefore the large a asymptotic approximation in equation (E7), leads again to an a c of order O(N −1 ), inconsistent with the initial assumption.

G Order parameter: the interface density ρ
We obtain, this section, an analytical expression for the order parameter ρ [equation (21) in the main text].ρ is defined as the interface density or density of active links, that is, the fraction of links connecting nodes in different states.In terms of the connectivity matrix A i j , ρ = Restricting our attention to the steady state average value of equation (G2), Nk 2 s i st + s j st − s i s j st , (G3) we can use the steady state mean solution found before for the individual node variables s i , s i st = 1/2, and the definition of the covariance matrix in the steady state, σ i j = s i s j st − 1/4, in order to write where we can identify the variable S 1 [see equation (D7)], Finally, reversing the relation (D25) between the variance of n and the variable S 1 , we can write the steady state average interface density ρ in terms of the variance of n, as it appears in the main text.

H Autocorrelation function of n
We derive here an analytical expression for the steady state autocorrelation function of n [equations ( 23) and ( 24) in the main text], defined as where τ plays the role of a time-lag.As far as the second point in time, t + τ, is concerned, we assume that the system was at n(t) at time t, and hence we can treat n(t) as an initial condition,  (H12) Given that we assume the state of the system at t to be our initial condition, b(t) can be written as and thus we have, for the autocorrelation function, 1 − e −(2a+h)τ ∑

(H14)
Using now the value found before for the steady state solution of the first-order moments, s i st = 1/2, and the definition of the covariance matrix in the steady state, σ i j = s i s j st − s i + e −(2a+h)τ ∑ i j )

= 2 •Figure 1 .
Figure 1.Fraction of nodes in state 1 on a Barabási-Albert scale-free network.Single realizations.The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.

Figure 2 .
Figure 2. Steady state variance of n as a function of the noise parameter a, for three different types of networks: Erdös-Rényi random network, Barabási-Albert scale-free network and dichotomous network.Symbols: Numerical results (averages over 20 networks, 10 realizations per network and 50000 time steps per realization).Solid lines: Analytical results [see equation 14].Dash-dotted lines: Analytical results for the critical points [see equation 18].Dashed line: Mean-field approximation (see42 ).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.

Figure 3 .
Figure 3. Steady state variance of n as a function of the noise parameter a for a Barabási-Albert scale-free network.Symbols: Numerical results (averages over 20 networks, 10 realizations per network and 50000 time steps per realization).Solid line:Analytical results [see equation(14)].Dotted line: asymptotic approximation for small a [see equation(16)].Dash-dotted line: asymptotic approximation for large a [see equation(17)].Dashed line: Crossover point between both asymptotic approximations (a * = 0.014157).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.

Figure 4 .
Figure 4. Steady state variance of n as a function of the variance of the degree distribution σ 2k for two values of the noise parameter a.In order to keep all parameters constant except the variance of the degree distribution, a different network type is used for each point (in order of increasing σ 2 k : Erdös-Rényi random network, Barabási-Albert scale-free network and dichotomous network).Circles with error bars: Numerical results (averages over 20 networks, 10 realizations per network and 50000 time steps per realization).Solid line and squares: Analytical results [see equation(14)].Dotted line: asymptotic approximation for small a [see equation(16)].Dash-dotted line: asymptotic approximation for large a [see equation(17)].Dashed line: Mean-field approximation (see42 ).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.

Figure 5 .
Figure 5. Critical value of the noise parameter a as a function of the variance of the degree distribution of the underlying network, σ 2k .In order to keep all parameters constant except the variance of the degree distribution, a different network type is used for each point (in order of increasing σ 2 k : Erdös-Rényi random network, Barabási-Albert scale-free network and dichotomous network).Symbols: Numerical results (averages over 20 networks, 10 realizations per network and 50000 time steps per realization).Solid line: Analytical results [see equation(18)].Dashed line: Mean-field approximation (see42,44 ).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.

= 2 •Figure 6 .
Figure6.Interface density on a Barabási-Albert scale-free network.Single realizations (the same realizations shown in Fig.1).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.

Figure 7 .
Figure 7. Steady state of the average interface density as a function of the noise parameter a in a linear-logarithmic scale and for three different types of networks: Erdös-Rényi random network, Barabási-Albert scale-free network and dichotomous network.Symbols: Numerical results (averages over 20 networks, 10 realizations per network and 50000 time steps per realization).Solid lines: Analytical results [see equation (21)].Dashed line: Mean-field pair-approximation (see44 ).The interaction parameter is fixed as h = 1, the system size as N = 2500 and the mean degree as k = 8.

Figure 8 .
Figure 8. Autocorrelation function of n in log-linear scale for a dichotomous network and a regular 2D lattice.Symbols: Numerical results (averages over 10 networks, 2 realizations per network and 200000 time steps per realization).Solid lines:Analytical results [see equation(23)].Parameter values are fixed as a = 0.01, h = 1, the system size as N = 2500 and the mean degree as k = 8.

)
Note that the remaining terms are now one order of N smaller than in the previous approximation [equation (E2)].To the first 22/28

2 (
[s i (1 − s j ) + (1 − s i )s j ] (s i + s j − s i s j ) annealed approximation for uncorrelated networks described in the main text [see equation(10)], we findρ = ∑ i j k i k j Nk(s i + s j − s i s j ) s i + s j − s i s j ) .