Irreducibility of multilayer network dynamics: the case of the voter model

Real-world interactions often happen at different levels and are therefore properly modelled by means of multilayer networks. Such multilayer approaches [1–3] have been applied to fields ranging from energy infrastructure [4] and transport [5–7], to epidemiology [8]. The multilayer set up can either describe interconnected networks with nodes of the same nature in each layer, but interacting with nodes of different nature in a different layer, or a multiplex structure with nodes of the same nature interacting via a different network in each layer. In any case, a central methodological question is that of multilayer reducibility, that is, when the multilayer framework is really needed to explain new phenomena, or when the system description can be reduced to an appropriately aggregated or reformulated single-layer network. An interesting contribution in this direction has recently come from the study of the structural reducibility of multilayer networks, i.e. of the possibility of aggregating some of the layers of a multi-dimensional network while preserving its distinguishability from the corresponding single-layer aggregated graph [9]. Although some recent works have pointed out that multiplex dynamics can be intrinsically different from their equivalent single-layer counterparts [10–14], little attention has been devoted so far to the problem of reducing a process taking place on a multilayer network to a dynamically equivalent process on an appropriate single-layer network. In this paper we address this general question, by implementing the voter model on multiplex networks and by studying the reducibility of its dynamics as a test case.

The voter model [15] is a nonequilibrium lattice model [16] which gives a standard framework for studying the influence of social imitation on the process of opinion formation [17]. A basic question considered in this context is when and how the system reaches an absorbing state with all the interacting nodes in the same state, or when an active dynamical situation of coexistence of different states prevails. The answer to these questions is known to depend crucially on the network structure and on the update rules employed [18–22]. The voter model has also been instrumental to understand fundamental phenomena in coevolution dynamics in which node states and network structure have coupled dynamical evolution with two different time scales [23–27]. In terms of comparison with real data, a metapopulation voter model has been recently shown to be able to account for voting patterns in the US general elections [28].

Quite frequently, social interactions happen in different concurrent contexts, so that any model of social imitation based on a single-layer representation of social relationships should probably be regarded only as a first-order approximation of a potentially more complex dynamics. In order to better mimic the multidimensional nature of social interactions, we consider here a multiplex voter model, where agents interact at two distinct layers. We assume also that the system exhibits a certain level of multiplexity, meaning that a fraction q of the agents is present on both layers [26]. As a consequence, any change in the state of those agents on either of the two layers is automatically transferred to the other layer, effectively coupling the voter model dynamics taking place on each of the two networks. By taking into account heterogeneity in the participation of agents on layers we aim at reproducing an interesting feature of real-world multilevel social systems [29]. In particular, it has been shown that the percentage of users of an online social network (e.g., Twitter) which have also an account on another online social network (e.g., Facebook or LinkedIn) lies somewhere between 30% and 70%, as reported by the social media matrix periodically published by the Pew Research Center [30]. These empirical findings confirm the necessity to consider systems whose multiplexity is neither 0 nor 1, but somewhere in between. A second important feature of social interactions is the fact that different layers often share a number of common links. This has been studied in details in online social systems. For instance, a significant overlap between the designation of other individuals as friends, correspondents, or trading partners has been found [31, 32]. Overlap is thus an intrinsic feature of real-world multilayer and interdependent networks, and its presence can substantially influence the system—see, for instance, studies of its role in percolation [33–35]. Here, we quantify this property of multilayer networks by means of the edge-overlap parameter introduced in [36].

In this work we address a fundamental question, that is whether the coupling of several voter models into a multiplex dynamics, motivated by the multiplicity of contexts influencing real interactions among individuals, gives rise to qualitatively different phenomena or effects than those observed in the classical single-layer setup. In particular, we want to establish whether multiplexity makes any difference for a simple social dynamics like the voter model, or if instead the multiplex voter model can indeed be reduced to an equivalent voter model dynamics on an appropriately constructed single-layer network.

A preliminary result on the question of the multilayer reducibility exists for a bilayer of uncorrelated networks evolving according to the voter model with adaptive links [26], but in the limit of zero edge overlap. In that model each layer is associated with a network plasticity parameter, that controls the rate at which relations among agents are rearranged, such that if the values of plasticity at the two layers are sufficiently different the system displays a network shattered fragmentation. On the other hand when the layers have the same plasticity one finds that the results coming from a pair approximation in the thermodynamic limit are equivalent to those of an appropriate aggregated single-layer network [26]. We build upon this, in time independent networks, examining both the asymptotic properties of the thermodynamic limit and the characteristic times to reach an absorbing state for finite systems. We also consider different schemes to obtain a possible equivalent aggregated single-layer network and highlight the importance of the degree of multiplexity and edge overlap.

The paper is organized as follows. Section 2 introduces our voter dynamics multilayer model. Section 3 describes our numerical findings. These are compared in section 4 with single-layer theoretical results for two natural methods of collapsing the multilayer into an aggregated single-layer network, and with numerical results obtained by considering a weighted single-layer aggregate. Section 5 analyses an optimum single-layer reduction method. Conclusions are discussed in section 6.

We consider here the case of a two-layer undirected unweighted multiplex network, described by the pair of binary adjacency matrices $\{{A}^{[1]},{A}^{[2]}\}$, where ${A}^{[\alpha ]}\equiv \{{a}_{{ij}}^{[\alpha ]}\}$, and ${a}_{{ij}}^{[\alpha ]}=1$ if and only if node i and node j are connected by a link at layer α, and zero otherwise. On each layer we have N nodes. A parameter of interest in this study is the average edge overlap ω, that is the probability that an edge is present on both layers:

$$\begin{eqnarray}&&\omega =\displaystyle \frac{\displaystyle \sum _{i,j}{a}_{{ij}}^{[1]}{a}_{{ij}}^{[2]}}{2K},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray*}&&K=\displaystyle \sum _{i}\displaystyle \sum _{j\gt i}\displaystyle \frac{{a}_{{ij}}^{[1]}+{a}_{{ij}}^{[2]}}{1+{a}_{{ij}}^{[1]}{a}_{{ij}}^{[2]}}\end{eqnarray*}$$

is the number of edges of the graph obtained by aggregating the two layers into a single one [36]. Notice that ω = 0 only if each edge exists in exactly one of the two layers, but not in the other one, while we have ω = 1 only if all the edges exist in both layers.

Each node i on layer α is associated to a binary state variable ${s}_{i}^{[\alpha ]}(t)$, where ${s}_{i}^{[\alpha ]}(t)$ can be either +1 or −1. Moreover, we assume that a fraction q of the N nodes participates in both layers, requiring that if i is one of these qN nodes then its state at the two layers will be identical at the end of every update. Such nodes participating in both layers are chosen randomly at initialisation. We can think of the parameter q as the interlayer connectivity or the degree of structural multiplexity of the system. The model is illustrated in figure 1.

Zoom In Zoom Out Reset image size

The multiplex voter dynamics consists of a sequence of time steps. During a time step we perform N updates, and each update consists of three elementary operations, as follows: (i) a layer α is selected at random, with uniform probability; (ii) one of the nodes i on layer α is chosen at random and its state ${s}_{i}^{[\alpha ]}$ is updated according to the classical voter model dynamics, that is, ${s}_{i}^{[\alpha ]}$ becomes the same as that of a randomly chosen neighbour of i on layer α; (iii) if the updated node i participates in both layers, then the state of the corresponding node in the other layer β changes as well by setting ${s}_{i}^{[\beta ]}={s}_{i}^{[\alpha ]}$. With this third operation, state changes can propagate across layers: indeed the presence of a fraction of nodes existing in both layers intertwines the voter dynamics on the two layers, so that in general the evolution of the overall multiplex dynamics might differ from the one we would observe on two independent networks of the same size. As a limiting case, the dynamics reduces to that of a classical voter model on a single-layer network only when q = 1.

It is well known that in connected finite-size single-layer networks the voter model dynamics always reaches an absorbing state, where all the nodes have exactly the same state, in a survival time that scales with the system size N[16–18, 20, 21]. In networks of high effective dimensionality (including random networks), when $N\to \infty$ the dynamics sustains an active disordered state in which nodes continue to change their state [19]. Such active state is the one observed asymptotically in large systems, before finite-size fluctuations pull the system towards the absorbing state. The classical order parameter to measure the activity in the voter model is the so-called interface density ρ(t), defined as the fraction of active edges of the network, i.e. of those edges whose endpoints have different states. In the multiplex model we can define:

$$\begin{eqnarray}&&{\rho }^{[\alpha ]}(t)=\displaystyle \frac{1}{2{K}^{[\alpha ]}}\displaystyle \sum _{i}\displaystyle \sum _{j\lt i}{a}_{{ij}}^{[\alpha ]}| {s}_{i}^{[\alpha ]}(t)-{s}_{j}^{[\alpha ]}(t)| ,\end{eqnarray} \tag{ 2 }$$

for each layer α, α = 1, 2, where we denote by ${K}^{[\alpha ]}=\displaystyle \frac{1}{2}{\displaystyle \sum }_{i,j}{a}_{{ij}}^{[\alpha ]}$ the total number of edges on layer α. Notice that ρ = 0 if and only if all the nodes have the same state, while larger values of ρ are associated to active configurations. A second quantity of interest for a finite system is the average time $\langle T\rangle$ to reach an absorbing state of consensus. Such a time can be defined as $\langle T\rangle ={\displaystyle \int }_{0}^{\infty }{P}_{s}(t){\rm{d}}t$ from the survival probability P_s(t), i.e. the probability for a system to be active at time t [37].

For single-layer uncorrelated networks at sufficiently large times t > N, P_s(t) decreases exponentially ${P}_{s}(t)\sim {{\rm{e}}}^{-2t/\tau }$ [38]. Hence, $\langle T\rangle \sim \tau /2$, so that the dependencies of $\langle T\rangle$ on system size and on the moments of the degree distribution are given precisely by the corresponding dependencies of the characteristic time τ. For uncorrelated networks one can find exact expressions for the value of the average interface density in the thermodynamic limit ρ^single, as well as for the characteristic time to reach the absorbing state τ^single [23, 38]. An important feature is that these analytical predictions of ρ^single and τ^single depend only on the values of the first two moments of the degree distribution P(k) of the network, i.e. on ${\mu }_{1}=\langle k\rangle ={\displaystyle \sum }_{k}{kP}(k)$ and ${\mu }_{2}=\langle {k}^{2}\rangle ={\displaystyle \sum }_{k}{k}^{2}P(k)$, and not on any other microscopic property of the network, and read:

$$\begin{eqnarray}&&{\rho }^{{\rm{single}}}=\displaystyle \frac{{\mu }_{1}-2}{3({\mu }_{1}-1)},\end{eqnarray} \tag{ 3 }$$

where ρ^single is the average over surviving runs, and

$$\begin{eqnarray}&&{\tau }^{{\rm{single}}}=\displaystyle \frac{({\mu }_{1}-1){\mu }_{1}^{2}N}{({\mu }_{1}-2){\mu }_{2}}.\end{eqnarray} \tag{ 4 }$$

In the following we will focus on the values ρ(q, ω) and τ(q, ω) of interface density and characteristic time of the multiplex voter model as a function of the average edge overlap of the system ω and of the fraction q of nodes present in both layers.

We studied the voter model on a duplex network made by two random regular graphs with N nodes, each node having a degree equal to μ. We therefore have ${\mu }_{1}=\langle {k}^{[1]}\rangle =\langle {k}^{[2]}\rangle =\mu$, and μ₂ = μ². In our simulations we fixed N = 1000 and μ = 4, and we studied the dynamics of the system by varying the value of the edge overlap ω and the fraction q of nodes participating in both layers. In order to obtain multiplex networks with different values of edge overlap, we started from a duplex consisting of two identical regular random graphs with L = Nμ/2 edges. Then, the desired value of edge overlap was obtained by rewiring the edges in only one of the two layers, say layer 2, and maintaining fixed the degree sequence. Notice that each edge might exists either on one of the two layers or on both. Hence, by rewiring on one layer an edge that exists on both layers, the number of edges present in both layers decreases by one unit, while the number of edges present in just one of the two layers increases by two units. Consequently, if we rewire a fraction r of the L edges of the second layer in such a way that each rewiring decreases the number of edges existing on both layers, we obtain a duplex with (1 + r)L distinct edges in total, of which (1 − r)L (i.e., those that have not been rewired) are present on both layers, while the $2{rL}$ remaining ones are present only on one of the two layers. Thus, the edge overlap in equation (1) can be rewritten as:

$$\begin{eqnarray}&&\omega =\displaystyle \frac{(1-r)L}{(1+r)L}=\displaystyle \frac{1-r}{1+r},\end{eqnarray} \tag{ 5 }$$

so that the fraction r of edges of the original duplex to be rewired in order to obtain a prescribed value of edge overlap ω is:

$$\begin{eqnarray}&&r=\displaystyle \frac{1-\omega }{1+\omega }.\end{eqnarray} \tag{ 6 }$$

The rewiring procedure described above is similar to the one used in [39] to obtain a two-layer duplex from a single-layer network. An important difference is that our procedure is based on the rewiring of single edges, while in [39] the model proceeds by relabeling nodes, and thus effectively rewiring a certain number of edges at each step while maintaining one of their endpoints fixed.

Since the system is finite, the dynamics will eventually converge to the absorbing frozen state corresponding to an interface density equal to zero. Consequently, for each time t, we will evaluate the average value of the interface density only over the surviving runs, i.e. on those realisations of the dynamics which are still active at time t. We associate the asymptotic value of this numerically obtained quantity with the value at the thermodynamic limit, given by equation (3), and refer to it as interface density, while being clear whether we are referring to the numerical or analytical results. In the case of a duplex with two layers, we have of course two values of interface density, one for each layer. If $\langle {k}^{[1]}\rangle =\langle {k}^{[2]}\rangle$ the ensemble averages on the two layers will be equal, and hence the activity of an arbitrary layer is representative of the typical activity of the entire system. We therefore use the asymptotic value of the interface density of an arbitrary layer to reflect the activity of the multiplex, and dispense with the α superscript.

In figure 2(a) we report the values of interface density ρ(q, ω) as a function of q, and for different values of ω. Let us consider first the case ω = 0 of no overlap between the structure of the two layers of the multiplex. When q = 0, i.e. when inter-layer state passing is not allowed, the system effectively corresponds to two identical but independent single-layer voter dynamics, so that ρ(q, ω) is in accordance with the classical analytical predictions for single-layer networks (since μ₁ = 4, equation (3) gives the activity ρ(0, 0) = 0.22). On the other hand, when q = 1, i.e. when all the nodes participate in both layers and each edge exists only in one layer, the system is in all respects identical to a single-layer network with μ₁ = 2 μ (for μ₁ = 8, equation (3) gives the activity ρ(1, 0) = 0.286). For intermediate values of q the dynamics interpolates monotonically between the two extreme cases, i.e. two dynamically indistinguishable voter models on single-layer networks with μ₁ = μ (q = 0), and one voter model on a single-layer network with μ₁ = 2 μ (q = 1). The picture changes completely as soon as ω is large enough. In general, when the overlap is above some limit, then ρ(q, ω) is a non-monotonic concave function of q, with a maximum at a given value of q in [0,1] which depends on ω. Notice that, when ω = 1 and q = 1, i.e. if the two layers are identical and all the nodes participate in both layers, the dynamics is identical to a voter model on a single-layer network with μ₁ = μ. In fact, since all the edges exist on both layers, a node participating in both layers will have only μ distinct neighbours, and will be connected to each of them on both layers. Hence, in a multiplex network with ω = 1, the interface density ρ of the voter model takes the same value 0.22 at q = 0 and at q = 1, while for intermediate values of q the interface density is higher than that of a voter model on each of the two layers.

Zoom In Zoom Out Reset image size

We now consider the characteristic time τ of the multiplex dynamics, where by characteristic time we understand twice the inverse exponent of the multiplex survival probability ${P}_{s}^{m}(t)$. We consider the multiplex as active if at least one of the layers is active, and find that for all q > 0 the survival probability of the multiplex does decay exponentially with some exponent τ = τ(q, ω) (figure 2(b)). The only exception is the q = 0 case of the fully disconnected multiplex, where the survival probabilities of the layers are independent. In this case the survival probability of the multiplex is given by the probability that at least one of the layers is active, ${P}_{s}^{m}(t,q=0,\omega =0)\sim {P}_{s}(t)(2-{P}_{s}(t))$, or ${P}_{s}^{m}(t,0,0)\sim 2{{\rm{e}}}^{-2t/{\tau }^{{\rm{single}}}}-{{\rm{e}}}^{-4t/{\tau }^{{\rm{single}}}}$. This means that for q = 0 the multiplex survival probability does not scale exponentially, and hence τ(q = 0) is not well-defined. Figure 2(c) shows the characteristic times τ(q, ω) obtained (the value at q = 0 is that of an approximate exponential fit). The trend shows a peak at small q, followed by an exponential decay with increasing number of interlayer connections. The value of edge overlap ω controls the rate of decrease. We stress that the peak is not a consequence of the definition of τ for the multiplex, and the consequent rogue value of τ at q = 0: the limit is truly singular, with plots for the actual average time until absorption displaying the same features, as do the plots for the characteristic exponent layer by layer. We also note that the peak is robust with respect to system size. Therefore the slowest finite-size multiplexes are ones where the layers are interconnected by the smallest number of links. Multiplexes with more interlayer connections will stop being active faster, as will wholly disconnected systems. This insight can be understood by realising that in a multiplex with few interlayer links, the layers will most of the time function as completely disconnected networks that half the time will try to settle into the absorbing states of consensus different from one another. Since this situation is now prohibited by the few interlayer links, the multiplex will freeze only when one of the layers switches over and both layers reach the same consensus. It is this switching behaviour that is responsible for the peak in τ for small q. We also notice that this nonlinear effect does not depend on the edge overlap ω. Here, unlike with the behaviour of interface density, the overlap does not change qualitatively the behaviour of τ for increasing q. In fact, the more interconnected the multiplex, the higher the role of edge overlap: tuning up ω causes a decrease in the interface density and results in a longer-lived multiplex, and the effect becomes more pronounced with q.

It was shown in [26] that the interface density of the multiplex voter model dynamics as a function of q and in absence of edge overlap can be rewritten as the interface density of a voter model on a single-layer network having ${\mu }_{1}=\langle k\rangle (1+q)$ (where k is the average degree of each of the two original layers), under some appropriately rescaled time. Here it is important to note that the approximation of [26] treated the interlayer connections as inherently probabilistic: the q parameter was a probability that each node's state gets passed on to the other layer. It was shown that the analytics of such a system, in the particular case ω = 0, are equivalent to those of the voter model on some properly aggregated single-layer network. However, as [26] only considered the thermodynamic limit, it is not obvious whether the aggregate displays a corresponding rescaling of characteristic time. Nevertheless, it should be possible to devise an aggregate that results from the flattening of the multiplex into a single-layer network, so that the resulting network will have a first moment of the degree distribution equal to the expression μ₁ given above.

The main question now is whether such a reducibility is possible in the most general case in which the multiplex has edge overlap $\omega \ne 0$. This is indeed the most interesting case, i.e. when in a multiplex there are correlations between the edges at the different layers [29, 36]. In the following our working hypothesis will be that, if the multiplex can indeed be reduced to a monoplex (in the variables of interest), then there exists an aggregate graph such that the behaviour of the voter model on the multiplex is completely described by equations (3) and (4) evaluated on the corresponding monoplex. Since those equations depend just on the number of nodes and on the first two moments of the degree distribution of the resulting single-layer network, we will now consider two standard aggregation procedures, and derive analytically the values of the first and second moments of the resulting degree distributions. Let us remind that, by definition, the first moment of the degree distribution of a graph is equal to ${\mu }_{1}=\langle k\rangle ={\displaystyle \sum }_{k=1}^{\infty }{kP}(k)$, where P(k) is the degree distribution, so that P(k = κ) is the probability that the degree of a node sampled at random from the graph is equal to κ. The probability P(k) can be also written as:

$$\begin{eqnarray}&&P(k)=\displaystyle \frac{{N}_{k}}{\displaystyle \sum _{{\ell }}{N}_{{\ell }}},\end{eqnarray} \tag{ 7 }$$

where N_k is the number of nodes in the graph having degree equal to k and the normalisation is just the total number of distinct nodes.

Given a duplex formed by two random regular graphs with identical degree μ, with edge overlap ω and where a fraction q of nodes participates in both layers, we can distinguish two classes of nodes. The nodes in the first class exist in only one layer, and we indicate as k_single the number of their neighbours, while nodes in the second class exist in both layers (these are the nodes having an inter-layer link) and we indicate their degree as k_both. The former class has degree:

$$\begin{eqnarray}&&{k}_{{\rm{single}}}=\mu ,\end{eqnarray} \tag{ 8 }$$

whereas nodes present on both layers have degree k_both = 2μ when ω = 0. However, if the edge overlap is not null, i.e. when ω > 0, then the nodes being present in both layers have degree k_both equal to:

$$\begin{eqnarray}&&{k}_{{\rm{both}}}=2\mu -q\mu \omega =\mu (2-q\omega ).\end{eqnarray} \tag{ 9 }$$

In fact, the degree of a node i present in both layers is equal to the sum of its degrees on the two layers (2μ) minus the expected number of its edges which are present on both layers. This number is equal to the probability that a neighbour j of i is also present on both layers (which is equal to q), times the probability that the edge (i, j) is present in both layers (which is equal to ω) multiplied by the number of neighbours of i (i.e., μ). Hence we get the correction qμω. In the particular case in which ω = 1, we get k_both = μ (2 −  q), while for ω = 0 we recover k_both = 2μ.

Let us now consider the two following distinct aggregation procedures. They are illustrated in figure 3 and correspond to the two most standard ways to aggregate a duplex into a single-layer network. The two flattening procedures differ in the total number of nodes and also in the number of nodes of degree k_single and k_both that they produce. This in turn changes the effective system size of the aggregate network, and the first two moments of the degree distribution. In the following sections we compute these quantities for the two methods, and assess how well the aggregates fare in describing the behaviour of the true multiplex. For the sake of completeness, we will also introduce a third aggregation method that results in a weighted single-layer network.

Zoom In Zoom Out Reset image size

4.1. Aggregation method 1

We can obtain a single-layer representation of a multiplex by putting the layers side by side with no effective node overlap. Therefore N_single = 2N(1 − q), N_both = 2Nq, and the total number of nodes in the aggregate is just N_aggr = 2N. The degrees of the nodes participating in just one or both layers are respectively equal to k_single and k_both. By using equation (7) we have:

$$\begin{eqnarray}&&P(k={k}_{{\rm{single}}})=(1-q),\end{eqnarray} \tag{ 10 }$$

and

$$\begin{eqnarray}&&P(k={k}_{{\rm{both}}})=q.\end{eqnarray} \tag{ 11 }$$

The first moment μ₁^aggr of the aggregated graph is then equal to:

$$\begin{eqnarray}{\mu }_{1}^{{\rm{aggr}}} & = & \displaystyle \sum _{k=1}^{\infty }{kP}(k)\\ & = & {k}_{{\rm{single}}}P(k={k}_{{\rm{single}}})+{k}_{{\rm{both}}}P(k={k}_{{\rm{both}}})\\ & = & {\mu }_{1}P(k={k}_{{\rm{single}}})+{\mu }_{1}(2-q\omega )P(k={k}_{{\rm{both}}})\\ & = & (1-q){\mu }_{1}+q(2{\mu }_{1}-q\omega {\mu }_{1})\\ & = & {\mu }_{1}(1+q-\omega {q}^{2}),\end{eqnarray} \tag{ 12 }$$

and the second moment

$$\begin{eqnarray}{\mu }_{2}^{{\rm{aggr}}} & = & {(1-q){\mu }_{1}^{2}+q{\mu }_{1}^{2}(2-q\omega )}^{2}\\ & = & {\mu }_{1}^{2}[1+3q-4\omega {q}^{2}+{\omega }^{2}{q}^{3}].\end{eqnarray} \tag{ 13 }$$

Under no overlap these reduce to

$$\begin{eqnarray}{\mu }_{1}^{{\rm{aggr}}} & = & {\mu }_{1}(1+q),\\ {\mu }_{2}^{{\rm{aggr}}} & = & {\mu }_{1}^{2}(1+3q).\end{eqnarray} \tag{ 14 }$$

Note that this method produces the same effective rescaling of the first moment μ₁ as given by the analytical estimation of the interface density of the multiplex in the thermodynamic limit (see [26]).

4.2. Aggregation method 2

This method reduces the number of effective nodes in the aggregate by treating nodes present in both layers as one node. So, while for the rest of the nodes we still have N_single = 2 N (1 − q), the number of 'multiplex' nodes is now half as much as in method 1, N_both = Nq. The total number of distinct nodes in the aggregate is equal to ${N}_{{\rm{aggr}}}=2N(1-q)+{Nq}=N(2-q)$. The degrees of the single and 'multiplex' nodes are as before equal to k_single and k_both respectively. Therefore

$$\begin{eqnarray}&&P(k={k}_{{\rm{single}}})=\displaystyle \frac{2(1-q)}{(2-q)},\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&P(k={k}_{{\rm{both}}})=\displaystyle \frac{q}{(2-q)}.\end{eqnarray} \tag{ 16 }$$

Consequently, we have for the first moment of the degree distribution ${\mu }_{1}^{{\rm{aggr}}}$:

$$\begin{eqnarray}{\mu }_{1}^{{\rm{aggr}}}=\displaystyle \sum _{k=1}^{\infty }{kP}(k) & = & \displaystyle \frac{2\;N(1-q)\times {\mu }_{1}+{Nq}\times {\mu }_{1}(2-q\omega )}{N(2-q)}\\ & = & {\mu }_{1}\displaystyle \frac{(2-{q}^{2}\omega )}{(2-q)}\end{eqnarray} \tag{ 17 }$$

and for the second moment ${\mu }_{2}^{{\rm{aggr}}}$:

$$\begin{eqnarray}&&{\mu }_{2}^{{\rm{aggr}}}={\mu }_{1}^{2}\displaystyle \frac{({q}^{3}{\omega }^{2}-4{q}^{2}\omega +2q+2)}{(2-q)}.\end{eqnarray} \tag{ 18 }$$

When the layers are uncorrelated and the overlap ω = 0, we have

$$\begin{eqnarray*}{\mu }_{1}^{{\rm{aggr}}} & = & \displaystyle \frac{2{\mu }_{1}}{(2-q)},\\ {\mu }_{2}^{{\rm{aggr}}} & = & 2{\mu }_{1}^{2}\displaystyle \frac{(1+q)}{(2-q)}.\end{eqnarray*}$$

We notice that, when q = 0 (i.e., if we have two non-interacting layers), we get ${\mu }_{1}^{{\rm{aggr}}}={\mu }_{1}$ and ${\mu }_{2}^{{\rm{aggr}}}={\mu }_{1}^{2}={\mu }_{2}$, while for q = 1 we have ${\mu }_{1}^{{\rm{aggr}}}=2{\mu }_{1}$ and ${\mu }_{2}^{{\rm{aggr}}}=4{\mu }_{1}^{2}={(2{\mu }_{1})}^{2}$.

4.3. Aggregation method 3

Up to this point we have considered the most straightforward procedures to create a single-layer network out of a multilayer system. The separate layers in our initial systems were unweighted, and by design so are the aggregates produced by the first two aggregation methods. However, if a link between two nodes (i, j) exists in both layers, then the probability for these two nodes to interact in the multiplex dynamics is twice as large than that associated to a pair of nodes which is connected only in one of the two layers. In this Subsection we explore whether an appropriately weighted single-layer aggregate network can provide a better approximation for the multiplex voter dynamics. In particular, we considered as a third possibility a single-layer weighted network with N(2 − q) nodes in total, as in the aggregate costructed with method 2, where the weight of edge (i, j) is either 1 (if that edge exists in only one of the two layers) or 2 (when the edge exists on both layers). Since the resulting network has edges with different weights, it is not immediate how to express the interface density as a function of the first and second moments of the degree sequence, or of the sequence of node strengths (where the strength of a node is equal to the sum of the weights of the edges incident on it). Hence, we investigated such a situation through numerical simulations. The update is implemented as follows. We randomly select a link with a probability proportional to its weight and, if the states of the two end-nodes are different, one of them is changed. For these homogenous networks this edge-centric evolution is equivalent to the node-centric one by which the voter model on the layers of the multiplex evolve.

We computed the interface density by averaging over surviving runs at some large enough time, when the trend has reached the asymptote. A systematic fitting of the profiles of interface density as a function of ω and q resulting from the simulations allowed to conclude that in the single-layer weighted aggregate graph ρ(q, ω), q > 0, can be well approximated as:

$$\begin{eqnarray}&&\rho (q,\omega )\approx {\rho }_{0}\left[1+\displaystyle \frac{q}{3}-\displaystyle \frac{{q}^{2}}{10}\left(1+\displaystyle \frac{3}{2}\omega +{\omega }^{2}\right)\right],\end{eqnarray} \tag{ 19 }$$

where ${\rho }_{0}=\displaystyle \frac{\mu -2}{3(\mu -1)}$ is the value of the interface density observed in an unweighted single-layer aggregate with average degree equal to μ when we set ω = 0 and q = 0.

4.4. Comparing aggregated to multiplex dynamics

We are now ready to compare the results of the numerical simulations of the multiplex system to the theoretical values given by the corresponding equations for the unweighted aggregates, and with the numerical simulations of the weighted aggregate. Equations (3) and (4) describe the behaviour of the monoplex. We contrast the two unweighted aggregates by putting in those equations the respective effective value for the first and second moments and the total number of nodes. In other words, we take $N\to {N}_{{\rm{aggr}}}$, ${\mu }_{1}\to {\mu }_{1}^{{\rm{aggr}}}$, and ${\mu }_{2}\to {\mu }_{2}^{{\rm{aggr}}}$, where the effective values differ depending on the aggregation method.

We consider first the behaviour in the thermodynamic limit, described by the interface density ρ. The results are reported in figure 4. It is interesting to notice that all the three methods yield aggregates whose qualitative behaviour with ω and q corresponds to the trend observed in the multiplex: for small q the system becomes more active by increasing the number of interlayer edges whereas, after a certain point, the activity can decrease, depending on whether or not the edge overlap is significant enough. However, none of the methods provide a quantitive explanation for a general (ω, q). Method 1 in general does a better job than method 2 at being systematically consistent, albeit predicting substantially higher values than the other two methods for intermediate values of q. As we have seen before, the analytical results for Method 1 (at ω = 0) correspond to a multiplex with probabilistic interlayer connections of intensity q. Therefore, at ω = 0, comparing the numerical results for the multiplex, and the analytical trend for method 1 can be used to gain insight into the difference induced by an alternative method of inteconnecting the layers of the multiplex (in fact, when viewed in this way, the small magnitude of the differences becomes more surprising than their presence). The main quantitative differences arise for a wide range of intermediate overlap values, for medium-to-large interlayer connectivity—precisely the region of parameter space motivated by real-world systems. We note that also the weighted aggregate shows a general qualitative agreement with the multiplex dynamics, but the quantitative discrepancies appear to be systematic, though somehow obscured by finite-size fluctuations. We conclude that considering the numerical simulations of the weighted aggregate instead of the analytical approximations of the two standard unweighted aggregates adds to, but does not substantially aid our understanding of the multiplex dynamics.

Zoom In Zoom Out Reset image size

The conclusion we draw from the analysis of the interface density is that, as long as edge overlap is taken into account, standard unweighted and weighted aggregates can only reproduce some qualitative features of large multiplex systems, but cannot capture quantitatively the behaviour of real multiplex networks.

Consider now the characteristic time of the approach to absorbing states of finite-size systems (figure 5). The difference between the results of the numerical simulations of the multiplex dynamics and the predictions at q = 0 of the unweighted aggregates is due to τ(q = 0, ω) not being defined for the multiplex; yet the aggregate timescale is close. The monotonic decrease for q > 0 exhibited by the numerical simulations is captured only by methods 2 and 3, with coincidence of τ values for q = 1. It is interesting to note that neither of the two unweighted aggregates is able to capture the jump of τ at small q, while the weighted aggregate can somehow reproduce this effect, at least qualitatively. Overall, while Method 1 appears to work best in the thermodynamic limit (this aggregation method results in a much slower system for large q than the observed duplex), Methods 2 and 3 do a better job in capturing the finite-size behaviour. In conclusion, although the three aggregation procedures can come close to describing the qualitative features of long-term duplex activity, they are not sophisticated enough to capture the long timescales associated with sparsely interconnected systems. In the following section we will see whether and how it is possible to devise more complex aggregation procedures in order to reproduce quantitatively the dynamics of the multiplex voter model.

Zoom In Zoom Out Reset image size

In principle, the results of the previous section cannot absolutely rule out the possibility that there exist other aggregation methods reproducing the behaviour of ρ(q, ω) observed in the multiplex voter model. Our hypothesis is that the deviations from the theory found for ω > 0 are due to the additional nonlinearity induced by inter-layer state copying made possible by the presence of a fraction of inter-layer edges. In order to better investigate these nonlinear effects, we inverted equation (3) to compute the effective value of the average degree ${\mu }_{1}^{{\rm{aggr}}}(q,\omega )$ of an ideal aggregated unweighted single-layer network able to reproduce the observed value of ρ(q, ω) for each value of ω and q. In formula:

$$\begin{eqnarray}&&{\mu }_{1}^{{\rm{aggr}}}(q,\omega )=\displaystyle \frac{3\rho (q,\omega )-2}{3\rho (q,\omega )-1}.\end{eqnarray} \tag{ 20 }$$

As made evident by figure 6, for each value of edge overlap, ${\mu }_{1}^{{\rm{aggr}}}(q,\omega )$ can be fitted very well by a quadratic polynomial in q. However, the actual values of the coefficients of the fit depend on ω in a non-trivial way. We can formally write:

$$\begin{eqnarray}&&{\mu }_{1}^{{\rm{aggr}}}(q,\omega )=F(\omega ){q}^{2}+G(\omega )q+H(\omega ).\end{eqnarray} \tag{ 21 }$$

Zoom In Zoom Out Reset image size

The problem now is to find an expression for F(ω), G(ω) and H(ω). We started by making an ansatz for the functional dependence of those coefficients on ω, and then fitting these functions, for each value of q, by using several realisations of ${\mu }_{1}^{{\rm{aggr}}}(q,\omega )$ corresponding to different values of ω.

We found that for all the values of ω the quadratic fit of ${\mu }_{1}^{{\rm{aggr}}}(q,\omega )$ yields H(ω) = μ, so we just focused on the other two coefficients. The values of F(ω) and G(ω) are reported in the panels of figure 7 as black circles. It is then reasonable to assume that F(ω), G(ω) are polynomial functions of ω. We found that in order to accurately reproduce the behaviour of ${\mu }_{1}^{{\rm{aggr}}}(q,\omega )$ in the whole range of ω, both F(ω) and G(ω) should be at least third-order polynomials in ω, as made evident by the plots reported in figure 7.

Zoom In Zoom Out Reset image size

Notice that the predictions of ${\mu }_{1}^{{\rm{aggr}}}$ provided by the two theoretical arguments reported above, based on the linear superpositions of the two layers, contained only linear terms in ω. However, the fit of ${\mu }_{1}^{{\rm{aggr}}}(q,\omega )$ confirms that, when ω > 0, the behaviour of the multiplex voter model is the result of a highly nonlinear combination of the two layers, suggesting that a trivial single-layer equivalent of the multiplex voter model dynamics does not exist, especially when the underlying multiplex network is characterised by a non-negligible overlap.

That simple aggregation procedures do not produce the observed scaling is also evident by examining the timescales on which finite systems approach absorbing states. For q > 0, the characteristic time τ(q, ω) is an exponentially decreasing function of q, $\tau (q,\omega )/N\sim {{\rm{e}}}^{-q\alpha (\omega ,N)}$, α > 0 (figure 2(b)). Figure 6(b) shows that $\alpha (\omega ,N)=-a\omega +b$, a > 0, and that the system-size dependency does not enter into it. Hence, $\tau (q,\omega )/N\sim {{\rm{e}}}^{q(\omega a-b)}$. Thus increasing the edge overlap ω results in longer-lived systems, while adding more interlayer links produces the opposite effect. This was additionally confirmed by examining the behaviour of the rescaled time until absorption $\langle T\rangle$, which showed the same qualitative trend and almost identical a and b coefficients.

Multilayer networks allow to extend the applicability of network theory to more realistic contexts in which nodes are connected through concurrent interaction patterns of different kinds. However, a fundamental open question to answer is whether the added complexity yielded by multilayer networks is really needed to model network phenomena, or if instead there exist simple ways of representing multiplex dynamics through appropriately constructed processes occurring on appropriately constructed single-layer networks. We have investigated here the problem of reducing the multiplex voter model to an equivalent single-layer dynamics. We have considered the predictions about the level of activity of the multiplex voter model in the thermodynamic limit, as measured by the interface density, as well as the time to reach the absorbing state for finite systems. We have found that results for the interface density based on single-layer aggregated graphs, either weighted or unweighted, are accurate only when there is little or no interaction between the layers (q sufficiently small), or when there is full connectivity (q ∼ 1) and the edge overlap ω is either 0 or 1. For the complementary broad range of parameters, any standard aggregation procedure can only give some qualitative information about the interface density, but fails to reproduce the quantitative details. We showed that edge overlap and multiplexity have two opposite effects on the long-term dynamics of the multiplex voter model, and in particular that an increase in the value of edge overlap can counter the action of increasing the fraction of interlayer links, leading to an overall decrease in the interface density of the multiplex. In fact, we showed numerically that any equivalent single-layer representation of the multiplex voter model dynamics entails the construction of an aggregate network which is a highly nonlinear combination of the original layers.

These results will be found at the same time surprising and interesting by all the researchers aiming at modelling social interaction in real-world scenarios. As a matter of fact, it has been recently shown that multilayer social networks are normally truly multiplex, meaning that they are characterised by non-negligible values of edge overlap [30] and by an intermediate level of multiplexity [29, 31, 32]. And as we have shown in this work, standard single-layer voter model approximations are quantitatively inaccurate in explaining the values of the interface density when the edge overlap ω and the degree of multiplexing q are far from their extreme values 0 or 1. The peculiarity of the multiplex dynamics is also reflected in the average timescales of consensus for finite systems, as we found that a multiplex with very few interlayer connections takes the longest time to reach consensus, much more than would probabilistically be needed by two disconnected monoplexes. This nonlinear effect is not captured by any of the simple unweighted aggregation procedures proposed in the paper, but it is qualitatively reproduced by the weighted aggregate. Hence, although the analytical treatment of the weighted aggregate, and its comparison with the multiplex systems, is beyond the scope of this work, we would still like to highlight it as a potentially fruitful direction for future research.

The case of dynamical irreducibility of a multiplex process presented in this work raises the important question of whether other unknown phenomena might be lurking in the multilayer structure of real-world systems. This question, together with the insights about the intrinsically multidimensional nature of the multiplex voter model, represents a stimulus to perform further research along these lines.

This work has been supported by the Spanish MINECO and FEDER under projects INTENSE@COSYP (FIS2012-30634), and by the EU Commission through the project LASAGNE (FP7-ICT-318132). VL also acknowledges support from EPSRC project GALE (EP/K020633/1).