Consensus ranking for multi-objective interventions in multiplex networks

In complex networks, a small subset of nodes often has disproportionate influence on the behavior of the system [1–4], and controlling such nodes allows us to steer the network to desired states [5–7]. For example, vaccinating a small fraction of carefully selected nodes can suppress large-scale disease outbreaks [8], or in a social network information is often disseminated by a small set of influencers [9]. Network centrality measures—such as degree, PageRank, or betweenness centrality—rank nodes based on their importance with respect to some process of interest. Therefore, high-centrality nodes provide effective targets to influence the behavior of a system: a node with high eigenvector centrality is an ideal target for vaccination [10], or a node in a social network with high closeness centrality is an effective influencer [1]. Nodes, however, typically participate in multiple networks simultaneously forming a multilayer or multiplex network [11–13]. For example, a person typically participates in a number of friendship and professional social networks, and these layers are often qualitatively different in nature, including affiliative and competitive interactions [14, 15]. Perturbing a node to influence one layer inadvertently affects the other layers that the node is active in, which can have drastic unintended consequences [5, 16–18]. Therefore, it is desirable to identify targets in a multiplex network that maximize our influence in certain layers, yet minimize any unwanted impact on others. More generally, we consider the scenario where the structure, the dynamics, and our intervention goals are qualitatively different for each layer of a multiplex network, and we aim to find nodes that effectively balance the trade-offs between the multiple objectives.

Recently, a number of centrality measures were introduced to rank nodes in multilayer or multiplex networks [19–22]. These centralities typically assume qualitatively similar layers with similar dynamics; and therefore they may provide useful intervention targets if the structure, the dynamics, and our intervention goals are all similar for each layer. For example, the family of 'versatility centralities' extend common single-layer centralities to multilayer networks; as such, PageRank versatility is useful to identify influential scientists in a multiplex co-authorship network [23–25]; or airports with high betweenness versatility are prone to congestion in the multiplex air traffic network [26]. Such multiplex centralities, however, do not address the more complex, and therefore less studied, scenario where the dynamics on the layers or our objectives are qualitatively different. For example, social and political multiplex networks are often comprised of cooperative (friendship, collaboration, alliance) and competitive (hostility, fighting, distrust) layers [14, 15, 27]; possible interventions might aim to reduce hostility while maintaining cooperation. In fact, this is a practical issue that caretakers of captive nonhuman primate groups face: sometimes unusually aggressive animals are removed with the goal of reducing overall levels of conflict, and it is desirable that this removal does not negatively affect the cohesion of the group [28–30]. For another example, consider the recent work modeling brain activity as a two-layer multiplex network with synchronization and transport dynamics. A possible strategy to affect the behavior of the system is to target nodes that are influential in both layers [31]. Or another recent work studied a two-layer SIS model where one layer represents the spread of a disease and the other the spread of awareness of that disease [32]. Although both layers are characterized by SIS dynamics, a possible intervention would have different goals for the two layers: blocking the propagation of the disease and promoting awareness; and the best spreaders and the best blockers are known to have different properties [33]. In these example multiplex systems, the function of each layer and/or our corresponding intervention goals are qualitatively different and ranking nodes in such networks has so far remained an uninvestigated problem.

In this paper, we explore the use of rank aggregation methods to identify target nodes for multi-objective interventions in multiplex networks. Instead of identifying influential nodes based on a single integrated multiplex centrality measure, we rank the nodes in each layer separately based on a centrality that is relevant to the function of the particular layer and our objectives. We then use rank aggregation methods to find a consensus ranking; high-ranking nodes are expected to balance all the objectives simultaneously. An expanding number of rank aggregation methods exist, and although these were originally studied by economists in the context of social choice theory [34], many also found applications in other disciplines, from computer science to biology [35, 36]. In this paper, we focus on two widely-used methods, Borda count and Kemeny aggregation [37, 38]. In section 2, we introduce consensus ranking in multiplex networks, and we demonstrate several of its properties using model and real networks. As an example application, in section 3, we study degree-based node removal of multiplex networks, such that we aim to destroy the largest connected component of some layers, while maintaining connectivity in the rest. In case of multiple objectives, optimality only exists in the Pareto-sense; therefore we rely on a generalized version of the Borda count algorithm that allows us to assign varying weight to the different objectives. With the weighted Borda count, we systematically study the trade-off between the competing goals, using analytically solvable model networks and a set of real networks.

A multiplex network consists of L layers, where each layer is a network ${{ \mathcal G }}_{\alpha }$ (α = 1, 2, ..., L) with the same set of nodes {v₁, v₂, ..., v_N}. These layers represent distinct types of interactions or relationships between the nodes, and there is no restriction on the structure of the individual layers, e.g. some layers might be directed, while others undirected, or they can be weighted or unweighted. To quantify the importance of a node i in an individual layer α, we calculate the node's centrality c_α,i as if layer α is in isolation. Nodes with high-centrality are important for the functioning of that given layer, and targeting such nodes is an effective strategy to influence or monitor the behavior of that particular layer. The type of centrality measure that works best depends on the function of the layer, for example, betweenness centrality is used for layers representing a transportation network, or eigenvector centrality is useful for layers governed by diffusion-like dynamics. However, nodes in a multiplex network do not perform a single task, but they simultaneously participate in all layers, prompting the question: How do we identify nodes that are important in all or most layers? This is particularly challenging if layers perform qualitatively different functions or our criteria for ranking nodes is different for each layer. To overcome these difficulties, instead of directly combining the centralities of each layer, we first determine the node ranks in each layer and then we use rank aggregation methods to identify a consensus ranking. An expanding number of rank aggregation methods exist [34–36], here we focus on two widely-used and intuitive methods: the Borda count (BC) and the Kemeny aggregation (KA) [37, 38].

We define the rank of node i in layer α as ${r}_{\alpha ,i}=N-{b}_{\alpha ,i}$, where

$$\begin{eqnarray}&&{b}_{\alpha ,i}=| \{j\in \{1,2,\,\ldots ,\,N\}:\,{c}_{\alpha ,i}\gt {c}_{\alpha ,j}\}| \end{eqnarray} \tag{ 1 }$$

is the number of nodes that have smaller centrality than node i. We denote the list of node ranks in layer α as ${r}_{\alpha }=({r}_{\alpha ,1},{r}_{\alpha ,2},\,\ldots ,\,{r}_{\alpha ,N})$. If there are no ties in the rankings (i.e. ${c}_{\alpha ,i}\ne {c}_{\alpha ,j}$ if $i\ne j$) r_α provides an ordinal ranking, the top rank being 1. If ties exist, the ranks of tied nodes are set to be equal and assigned the worst value, e.g. if nodes i and j have equal centralities c_α,i = c_α,j and are top ranked, we assign r_α,i = r_α,j = 2 to both nodes.

The BC algorithm is perhaps the most straight-forward rank aggregation method [38], it works by assigning a score b_i to each node i equal to the average number of nodes that i outranks in the layers of the multiplex:

$$\begin{eqnarray}&&{b}_{i}=\displaystyle \frac{1}{L}\displaystyle \sum _{\alpha =1}^{L}{b}_{\alpha ,i},\end{eqnarray} \tag{ 2 }$$

where b_α,i is defined in equation (1). The BC ranking, ${R}_{\mathrm{BC}}$, is simply a ranking of the nodes based on the scores b_i. The algorithm requires sorting the nodes; therefore its computational complexity is $O(N\mathrm{log}N)$.

The KA algorithm belongs to a larger class of rank aggregation methods that aim to identify a consensus ranking R such that R is at the minimum average distance from the rankings r_α [39, 40]. Formally, it is defined as

$$\begin{eqnarray}&&R={\mathrm{argmin}}_{s}\displaystyle \frac{1}{L}\displaystyle \sum _{\alpha =1}^{L}D({r}_{\alpha },s),\end{eqnarray} \tag{ 3 }$$

where D(r, s) is some measure of distance between two rankings r and s, such as Hamming distance, Cayley distance, etc The KA algorithm uses the Kemeny distance [37], which naturally takes into account ties:

$$\begin{eqnarray}&&D(r,s)=\displaystyle \frac{1}{N(N-1)}\displaystyle \sum _{i\lt j}| r(i,j)-s(i,j)| ,\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}r(i,j)=\left\{\begin{array}{ccc}1, & \mathrm{if} & {r}_{i}\lt {r}_{j},\\ 0, & \mathrm{if} & {r}_{i}={r}_{j},\\ -1, & \mathrm{if} & {r}_{i}\gt {r}_{j},\end{array}\right.\end{eqnarray} \tag{ 5 }$$

and s(i, j) is defined similarly. Kemeny distance is normalized such that if r and s are the same, D(r, s) = 0; and if r and s rank nodes in opposite order, D(r, s) = 1. If there are no ties in the rankings, the Kemeny distance is equivalent to Kendall's distance⁶ .

Unfortunately, identifying the Kemeny consensus is an NP-complete problem for L ≥ 4 [41]; therefore approximate methods have to be used. Over the years, a high number of such algorithms were proposed: a recent survey compared 104 algorithms and combinations of algorithms [42]. In fact, the previously introduced BC algorithm can be considered as a simple heuristic approximation of the KA problem. Here, we implemented a local search algorithm, which was identified as providing optimal trade-off between accuracy and run-time [42]. The algorithm starts from ${R}_{\mathrm{BC}}$ and outputs a new consensus ranking ${R}_{\mathrm{KA}}$ by finding a local minimum of equation (3) using a restricted set of transformations. The computational complexity of our implementation is O(N²); we provide the details of the algorithm in appendix A.

2.1. Model networks

We now compare the two algorithms by applying them to a multiplex network model with tunable pairwise Kemeny distance between its layers. To generate these multiplex networks, we first independently create L layers, and we determine the rank r_α,i of each node i in each layer α with respect to some centrality. We then generate a node label sequence for each layer to ensure that the distance between each pair of layers is D. Here, we use the Erdős–Rényi model (ER) or the scale-free static model (SF) to generate each layer, and we use degree centrality to rank the nodes. We describe the procedure in detail in appendix B.

By choosing parameter D of the model network, we set the strength of consensus: if D = 0.5, the layers are independent and there is no intrinsic consensus; if D = 0, the ranking in each layer is the same and consensus is perfect. We compare the two algorithms by calculating the average distance of the layers from consensus:

$$\begin{eqnarray}&&{D}_{\mathrm{con}}^{\mathrm{BC}/\mathrm{KA}}=\displaystyle \frac{1}{L}\displaystyle \sum _{\alpha =1}^{L}D({r}_{\alpha },{R}_{\mathrm{BC}/\mathrm{KA}}),\end{eqnarray} \tag{ 6 }$$

where ${R}_{\mathrm{BC}/\mathrm{KA}}$ is the consensus ranking found by the BC or KA algorithm. Note that ${D}_{\mathrm{con}}$ corresponds to the cost function of the KA algorithm provided by equation (3), and we chose ${D}_{\mathrm{con}}$ for comparison because BC is sometimes considered as an approximation of KA [42]. The results, however, have to be interpreted with care, finding a ranking that corresponds to lower ${D}_{\mathrm{con}}$ does not necessarily mean that it is better, because ultimately BC and KA define the consensus in different ways. Similarly to the fact that no single best centrality measure exists to rank the nodes, the preferred rank aggregation algorithm also depends on our particular purposes.

Figure 1 shows ${D}_{\mathrm{con}}$ as a function of the distance between layers D. Since our implementation of the KA algorithm works by improving the solution of the BC algorithm, we always find that ${D}_{\mathrm{con}}^{\mathrm{KA}}\leqslant {D}_{\mathrm{con}}^{\mathrm{BC}}$ (the equality holds for L = 2, where both algorithms find the same consensus, corresponding to the global minimum of ${D}_{\mathrm{con}}$). However, this improvement is marginal in terms of ${D}_{\mathrm{con}}$. If strong consensus exists (D ≈ 0) both algorithms closely approximate this consensus. Even if the layers are independent (D ≈ 0.5), we are able to find a consensus with ${D}_{\mathrm{con}}\lt 0.5$, especially for multiplexes with only a few layers. For the latter case, the consensus does not capture an intrinsic property of the multiplex network, rather it identifies nodes that by chance have high rank in many layers; and therefore potentially provide low-cost targets for influencing multiple layers simultaneously.

Zoom In Zoom Out Reset image size

Note that figure 1 shows results for ER layers, we found that using the SF layers produces almost indistinguishable results. In fact, the performance of the algorithms only depends on the rankings r_α, and not on the details of each layer α. The only property of r_α that is specific to ER or SF layers is the number of ties in the ranking, and our results indicate that the performance of the BC and KA algorithms in terms of ${D}_{\mathrm{con}}$ are insensitive to this.

2.2. Real networks

The BC and KA algorithms allow us to identify and analyze consensus ranking in real systems represented by multiplex networks. Here, we investigate three examples:

• 1.  
  Airline network. A multiplex network with 5 layers representing the United States air traffic network in 2013. The nodes are airports and links indicate direct flights between them, and the layers correspond to the 5 largest carriers⁷ .
• 2.  
  Primate social network. A multiplex network with 4 layers containing interaction data from one week of observations of a captive rhesus macaque troupe. The layers represent different interactions between animal pairs, including conflict, signaling of subordination, grooming, and huddling [43, 44].
• 3.  
  Human social network. A multiplex network with 5 layers, where nodes are members of the Department of Computer Science at the Aarhus University and the links indicate various social relationships: Facebook friendship, spending leisure time together, working together, co-authorship, and regularly having lunch together [45].

Table 1 summarizes some basic properties of these networks. The airline network has a clear heavy-tailed degree distribution with hubs that have significantly more connections than average nodes. It is not possible to be as definite about the social networks due to their small size; there are, however, nodes that are connected to a significant fraction of the network, for example, in the fifth layer of the human social network the largest hub is connected to almost half of the other nodes. The primate social network is composed of competitive interactions (conflict and subordination) and cooperative interactions (grooming and huddling)—the existence of antagonistic and affiliative layers is a general property of social and political networks [14, 15, 27]. For the primate network, we find that the adversarial layers are more heterogeneous than the affiliative layers. Interestingly, similar pattern was observed for human social networks obtained from an online game that allowed competition and alliances between players [46].

Table 1.  Properties of the example multiplex networks. We provide the number of nodes (N) and layers (L), the average degree of each layer ($\left\langle k\right\rangle$); and we quantify the degree heterogeneity of each layer with $\kappa =\left\langle {k}^{2}\right\rangle /\left\langle k\right\rangle$ (for Poisson degree distribution $\kappa =1+\left\langle k\right\rangle$ and for scale-free networks $\kappa \to \infty$ as $N\to \infty$) and with the maximum degree ${k}_{\max }$.

	N	L	$\left\langle k\right\rangle$	$\kappa =\langle {k}^{2}\rangle /\left\langle k\right\rangle$	${k}_{\max }$
Airline	575	5	4.63, 2.95, 9.02,	28.7, 39.1, 37.2,	109, 126, 125,
			5.99, 4.11	31.9, 28.3	95, 66
Primate Soc.	100	4	8.62, 4.90,	12.1, 8.9	35, 29,
			4.32, 4.98	5.8, 6.3	13, 12
Human Soc.	61	5	6.33, 4.07, 0.69,	7.8, 9.5, 2.3,	15, 15, 5,
			2.89, 6.36	5.7, 10.8	14, 27

When ranking nodes in a network, we chose a node centrality depending on the function of the system and our goals. Here, for illustration purposes, we calculate the node rankings for each network based on degree centrality. Figures 2(a)–(c) shows the pairwise Kemeny distances between layers for each network. Typically, we find D < 0.5, indicating positively correlated rankings; the correlation is the strongest in the airline network, while rankings in the social networks are less aligned with each other. We observe an interesting pattern in the primate social networks: the distance is low between affiliative layers (grooming and huddling), and also low between the interactions involving social hierarchy (conflict and signaling); but the rankings corresponding to competitive layers are independent from rankings of affiliative layers.

Zoom In Zoom Out Reset image size

We also identify the consensus ranking for each network using the BC and the KA algorithms. Figure 2(d) shows the average distance of the layers from consensus, ${D}_{\mathrm{con}}$, for the original networks and their randomized counterparts, where we shuffled the node labels in each layer, eliminating any correlation between the rankings. For all multiplexes, we find stronger consensus in the real instances than in their randomized versions: the largest difference was observed in the airline network, while smallest in the primate network. As expected, KA always finds a consensus with a lower ${D}_{\mathrm{con}}$ than BC; the difference, however, is marginal with the exception of the airline network.

Overall, we found for both model and real multiplexes that the BC and KA consensus rankings are similar in terms of ${D}_{\mathrm{con}}$, while the BC algorithm is faster and extremely simple to implement. Both algorithms effectively identify the consensus if the layers are strongly correlated, and even if the layers of a multiplex are independent, we can identify a consensus ranking that provides low-cost targets for simultaneous intervention on multiple layers. So far, we only compared the consensus rankings in terms of ${D}_{\mathrm{con}}$, in the next section we demonstrate how consensus ranking can identify targets for a simple degree-based node removal model with multiple objectives.

In this section, we explore an example of using consensus ranking to identify effective targets for multi-objective interventions in multiplex networks. A classic result of network science is that complex networks with heavy-tailed degree distributions are robust against random node removal; while targeted removal of high-degree hubs rapidly breaks down their large-scale connectivity [3]. Complex networks, however, rarely exist in isolation, nodes typically participate in multiple networks simultaneously. Therefore, removing a node from a multiplex network to reduce connectivity in one layer, might also remove connections from other layers that otherwise we would like to preserve. More specifically, here, we consider the problem of removing nodes from a multiplex network with L layers such that for a set of layers ${{ \mathcal L }}_{{\rm{D}}}$ we aim to reduce the size of the largest component, while we want to keep the rest of the layers ${{ \mathcal L }}_{{\rm{K}}}$ intact. We rank the nodes in each layer based on degree centrality, but to reflect our different objectives, we reverse the rankings for the layers that we keep intact:

$$\begin{eqnarray}\begin{array}{rcl}{r}_{\alpha ,i}^{{\rm{D}}} & = & N-| \{j\in \{1,\,\ldots ,\,N\}:\,{k}_{\alpha ,i}\gt {k}_{\alpha ,j}\}| ,\\ {r}_{\alpha ,i}^{{\rm{K}}} & = & N-| \{j\in \{1,\,\ldots ,\,N\}:\,{k}_{\alpha ,i}\lt {k}_{\alpha ,j}\}| ,\end{array}\end{eqnarray} \tag{ 7 }$$

where ${r}_{\alpha ,i}^{{\rm{D}}}$ and ${r}_{\alpha ,i}^{{\rm{K}}}$ are the ranks of node i in layers that we want to destroy and keep, respectively; and k_α,i is the degree of node i in layer α. Given a multiplex network, we calculate these ranks for each layer, and we use the BC and KA algorithms to identify consensus rankings. We then iteratively remove the nodes based on this consensus ranking starting with the highest ranked nodes.

Figure 3 shows the relative size of the largest connected component S(f) as a function of the fraction of nodes removed f for model multiplex networks with L = 4 ER or SF layers (for details about the model networks see appendix B). For the ER example, we find that both the BC and KA consensus-based removal reduces S faster for the layers we aim to destroy, and slower for layers we aim to preserve than random node removal, but BC preforms better at reducing S, while KA is better at preserving connectivity. Similarly for the SF example, BC destroys the targeted layers faster than KA, but does not keep the rest of the layers intact more than random removal would. In either case, we are not as effective as if the layers would be in isolation.

Zoom In Zoom Out Reset image size

Comparing the BC and KA rankings in figure 3, we find that neither method is objectively better than the other, instead they provide a different trade-off between the two objectives: BC preforms better at destroying layers, while KA does better at keeping layers intact. In fact, this is a pattern that we widely observed varying the parameters of the multiplex model networks. Therefore, from hereon in this section, we will focus on the BC ranking and we will study the trade-off by introducing a weighted version of BC. Further reasons for focusing on the BC algorithm are that it is much faster than the KA and the simplicity of the BC algorithm allows us to analytically solve S(f) for the model networks.

The original BC algorithm provides one possible trade-off between the competing objectives of destroying certain layers, while keeping others intact. To explore other possible trade-offs, we introduce a weighted Borda count (wBC) algorithm that allows us to assign varying preference to the different objectives. The weighted Borda score of node i is defined as

$$\begin{eqnarray}&&{b}_{i}=\displaystyle \frac{1}{L}\displaystyle \sum _{\alpha =1}^{L}{w}_{\alpha }{b}_{\alpha ,i},\end{eqnarray} \tag{ 8 }$$

where w_α = w, if we destroy layer α ($\alpha \in {{ \mathcal L }}_{{\rm{D}}}$); w_α = 1 − w, if we keep layer α intact ($\alpha \in {{ \mathcal L }}_{{\rm{K}}}$); and $w\in [0,1]$. The choice w = 0.5 corresponds to the unweighted Borda score (with a multiplier of 1/2); if w = 1, we only care about the layers we aim to destroy; and if w = 0, we only care about the layers we aim to keep intact. In the following, we first derive an analytical solution for S(f), and then we systematically investigate how multiplex network structure affects the trade-off between the different objectives using model and real networks.

3.1. Analytical solution

To analytically solve the size of the giant connected component (GCC) for consensus-based removal, we first calculate the size of the GCC for a general degree-based removal strategy (DBS) on a single-layer network, and then we map the consensus-based process to a degree-based one for each layer. By DBS we mean a node removal process where the probability of removing a node only depends on its degree, i.e. the probability of removing node i is f(k_i). Let s be the probability that a randomly selected link leads to the GCC. Assuming local tree-like structure and uncorrelated networks, we calculate s using the self-consistent equation

$$\begin{eqnarray}&&1-s=\displaystyle \sum _{k=1}^{\infty }\displaystyle \frac{k}{c}\left(1-f(k)\right)p(k){\left(1-s\right)}^{k-1}+\displaystyle \sum _{k=1}^{\infty }\displaystyle \frac{k}{c}f(k)p(k),\end{eqnarray} \tag{ 9 }$$

where p(k) is the degree distribution and c is the average degree of the network. The first term on the right hand side is the probability that following a random link leads to a node that is not removed and not part of the GCC, and the second term is the probability that the node is removed. Once s is determined, we obtain the relative size of the GCCs using the equation

$$\begin{eqnarray}&&1-S=\displaystyle \sum _{k=0}^{\infty }\left(1-f(k)\right)p(k){\left(1-s\right)}^{k}+\displaystyle \sum _{k=0}^{\infty }f(k)p(k).\end{eqnarray} \tag{ 10 }$$

Several node removal processes can be described as a DBS. For example, the choice $f(k)\equiv f$ leads to simple random node removal. Generally the best DBS to destroy a single-layer network removes nodes starting with the highest degree. Formally we express this as

$$\begin{eqnarray}f(k)=\left\{\begin{array}{ccc}1, & \mathrm{if} & k\gt K,\\ \displaystyle \frac{f-{\displaystyle \sum }_{k\gt K}p(k)}{p(K)}, & \mathrm{if} & k=K,\\ 0, & \mathrm{if} & k\lt K,\end{array}\right.\end{eqnarray} \tag{ 11 }$$

meaning that we remove all nodes with degree higher than K and a certain fraction of nodes with degree K, where f is the total fraction of nodes removed, $K=\inf \{k:1-f\leqslant P(k)\}$ and P(k) is the CDF of the degree distribution. Also, the most effective DBS to keep a single-layer network intact is to remove nodes with the lowest degree first:

$$\begin{eqnarray}f(k)=\left\{\begin{array}{ccc}0, & \mathrm{if} & k\gt K,\\ \displaystyle \frac{f-{\displaystyle \sum }_{k\lt K}p(k)}{p(K)}, & \mathrm{if} & k=K,\\ 1, & \mathrm{if} & k\lt K,\end{array}\right.\end{eqnarray} \tag{ 12 }$$

where $K=\inf \{k:f\leqslant P(k)\}$.

To map the consensus-based node removal to a DBS, we first calculate the Borda score of a node with multiplex degree (k₁, ..., k_L) as

$$\begin{eqnarray}&&b({k}_{1},\,\ldots ,\,{k}_{L})=\displaystyle \frac{w}{L}\displaystyle \sum _{\alpha \in {{ \mathcal L }}_{{\rm{D}}}}{P}_{\alpha }({k}_{\alpha }-1)+\displaystyle \frac{1-w}{L}\displaystyle \sum _{\alpha \in {{ \mathcal L }}_{{\rm{K}}}}\left(1-{P}_{\alpha }({k}_{\alpha })\right),\end{eqnarray} \tag{ 13 }$$

where P_α(k) is the CDF of the degree distribution of layer α. Next, we calculate the Borda score distribution, i.e. the probability that a randomly selected node of the multiplex has Borda score b:

$$\begin{eqnarray}&&{p}_{{\rm{B}}}(b)=\displaystyle \sum _{{k}_{1},\,\ldots ,\,{k}_{L}=0}^{\infty }p({k}_{1},\,\ldots ,\,{k}_{L})\delta \left(b-b({k}_{1},\,\ldots ,\,{k}_{L})\right),\end{eqnarray} \tag{ 14 }$$

where p(k₁, ..., k_L) is the multiplex degree distribution. When removing f fraction of nodes based on consensus ranking, the probability of removing a node with Borda score b is

$$\begin{eqnarray}f(b)=\left\{\begin{array}{ccc}1, & \mathrm{if} & b\gt B,\\ \displaystyle \frac{f-{\displaystyle \sum }_{b\gt B}{p}_{{\rm{B}}}(b)}{{p}_{{\rm{B}}}(B)}, & \mathrm{if} & b=B,\\ 0, & \mathrm{if} & b\lt B,\end{array}\right.\end{eqnarray} \tag{ 15 }$$

where $B=\inf \{b:1-f\leqslant {P}_{{\rm{B}}}(b)\}$ and ${P}_{{\rm{B}}}(b)$ is the CDF corresponding to ${p}_{{\rm{B}}}(b)$ provided by equation (14). For each layer α, the consensus-based removal can be formulated as an effective DBS with

$$\begin{eqnarray}&&{f}_{\alpha }({k}_{\alpha })=\displaystyle \frac{{\displaystyle \sum }_{\{{k}_{\beta }:\beta \ne \alpha \}}p({k}_{1},\,\ldots ,\,{k}_{L})f(b({k}_{1},\,\ldots ,\,{k}_{L}))}{{\displaystyle \sum }_{\{{k}_{\beta }:\beta \ne \alpha \}}p({k}_{1},\,\ldots ,\,{k}_{L})},\end{eqnarray} \tag{ 16 }$$

where p(k₁, ..., k_L) is the joint multiplex degree distribution, and the summations are over the degrees in all layers, except α.

Substituting equation (16) into (9) and (10) provides S_α(f), the relative size of the GCC of layer α for consensus-based node removal. We numerically evaluate these equations for a class of multiplex model networks where the layers are either ER or SF networks, for SF networks we use the degree distribution provided in equation (B6). Figure 4 compares the analytical solution of S_α(f) to simulations, showing excellent agreement. We find that there is a non-trivial connection between S_α(f) and f as the node removal process switches back-and-forth between destroying and preserving layers. Note that the methods that we used to numerical evaluate the necessary equations become intractable for increasing number of layers; therefore for the numerical solutions we restrict ourselves to L = 4 ER and L = 2 SF layers.

Zoom In Zoom Out Reset image size

3.2. Model networks

The w parameter of the wBC algorithm allows us to assign different level of importance to the two competing objectives: $w\approx 1$ favors breaking down layers, while w ≈ 0 focuses on keeping layers intact. To systematically investigate this trade-off, we introduce the coefficients

$$\begin{eqnarray}\begin{array}{rcl}{O}_{{\rm{D}}} & = & \displaystyle \frac{1}{| {{ \mathcal L }}_{{\rm{D}}}| }\displaystyle \sum _{\alpha \in {{ \mathcal L }}_{{\rm{D}}}}\displaystyle \frac{{S}_{\alpha }-{S}_{\alpha }^{\min }}{{S}_{\alpha }^{\max }-{S}_{\alpha }^{\min }},\\ {O}_{{\rm{K}}} & = & \displaystyle \frac{1}{| {{ \mathcal L }}_{{\rm{K}}}| }\displaystyle \sum _{\alpha \in {{ \mathcal L }}_{{\rm{K}}}}\displaystyle \frac{{S}_{\alpha }^{\max }-{S}_{\alpha }}{{S}_{\alpha }^{\max }-{S}_{\alpha }^{\min }},\end{array}\end{eqnarray} \tag{ 17 }$$

where S_α is the size of the GCC in layer α after removing f fraction of its nodes using the wBC algorithm; ${S}_{\alpha }^{\min }$ (${S}_{\alpha }^{\max }$) is the size of the GCC after removing f fraction of its highest (lowest) degree nodes, i.e. it is the best we could do if layer α was in isolation. The coefficients ${O}_{{\rm{D}}}$ and ${O}_{{\rm{K}}}\in [0,1]$ measure how well the objectives are achieved relative to the case when we remove nodes from each layer independently. Analytical solution of equation (17) is obtained by first calculcating ${S}_{\alpha }^{\min },{S}_{\alpha }^{\max }$, and S_α using equations (11), (12), and (16).

Figure 5 shows the trade-off curve between ${O}_{{\rm{D}}}$ and ${O}_{{\rm{K}}}$ for model multiplex networks with L layers, where we aim to destroy half of the layers, while keeping the other half intact. If there would be no trade-off, the curve would be a single point at $({O}_{{\rm{K}}},{O}_{{\rm{D}}})=(1,1);$ if all layers would be identical, the trade-off curve would be the diagonal line ${O}_{{\rm{D}}}=1-{O}_{{\rm{K}}}$.

Zoom In Zoom Out Reset image size

Figure 5(a) shows that increasing the number of layers affects the trade-off in two distinct ways: (i) as L increases the number of conflicting objectives also increase, and the trade-off becomes more severe approaching the diagonal line. (ii) For $L\geqslant 4$, we aim to destroy multiple layers simultaneously, which cannot be done as efficiently as if they where in isolation; therefore even if $w=1,{O}_{{\rm{D}}}$ remains less than 1. Increasing L has a similar, but weaker effect on ${O}_{{\rm{K}}}$, since destroying a network is more difficult than keeping it intact. Changing the pairwise Kemeny distance D between the layers (figure 5(b)), we find that strong consensus (D ≈ 0) leads to strong trade-off between ${O}_{{\rm{D}}}$ and ${O}_{{\rm{K}}}$, but allows efficient destruction (${\max }_{w}{O}_{{\rm{D}}}\approx 1$) or preservation (${\max }_{w}{O}_{{\rm{K}}}\approx 1$) of layers. Figure 5(c) shows the effect of degree heterogeneity, we find that destroying SF layers and preserving ER layers entails the least amount of trade-off, while destroying SF and preserving ER is the most difficult. Finally, figure 5(d) shows that the trade-off is the most pronounced at the initial stages of node removal (small f), and the competing objectives are less restrictive in case of large-scale removals (large f).

3.3. Real networks

In the previous section, we used model multiplex networks to understand how basic network properties affect the ${O}_{{\rm{D}}}$–${O}_{{\rm{K}}}$ trade-off. Real networks, however, have more complex structure. To investigate their effect, we use the three example datasets introduced in section 2.2: the human social network, the primate social network, and the airline network. We assign half of the layers in each multiplex to be destroyed and the rest to be kept intact, and we calculate the ${O}_{{\rm{D}}}$–${O}_{{\rm{K}}}$ trade-off after removing a fraction of the nodes (f = 0.4 for the airline network and the human social network, and f = 0.55 for the primate social network). We then compare the trade-off curves to the following randomized null models:

• 1.  
  Full randomization (FR): we randomly rewire all links, in effect replacing each layer with an ER network with the same number of nodes and links.
• 2.  
  Node label randomization (NLR): we shuffle the node labels in each layer, removing any correlation or consensus between layers, but leaving the structure of each layer unchanged.
• 3.  
  Degree preserved randomization (DPR): we rewire links such that the multiplex degree of each node is unchanged. This randomization, therefore, preserves the pairwise Kemeny distances between the layers, and removes all structure within the layers beyond their degree sequence.
• 4.  
  Node label and degree preserved randomization (NLR + DPR): combining NLR and DPR removes correlations both within layers and between layers, and only preserves the degree distributions of the individual layers.

Figure 6 provides the trade-off curves for the three example networks, each network showing a distinct behavior. For the airline network, we aim to destroy two layers and keep the other three intact. In section 2.2, we showed that the airline network has heterogeneous degree distributions and strong consensus between its layers. For model networks, we found that strong consensus leads to significant ${O}_{{\rm{D}}}$–${O}_{{\rm{K}}}$ trade-off, and indeed, we find that for the airline network the trade-off curve is close to the diagonal. We find that the DPR trade-off curve is almost indistinguishable from the original, while the trade-off for FR, NLR, and NLR + DPR is significantly weaker, meaning that both the inter-layer correlations and the degree distributions are necessary to explain the strong trade-off. We also find that ${\max }_{w}{O}_{{\rm{D}}}=1$ for the original and NLR curve, but not for FR, this indicates that the simultaneous destruction of the two layers is aided by the presence of hubs and is uneffected by the inter-layer correlations.

Zoom In Zoom Out Reset image size

In the case of the primate social network, we aim to destroy the two layers that are related to competition (conflict and signaling) and to preserve the affliative layers (grooming and huddling). In fact, it is common practice for the caretakers of captive primate groups to remove unusually aggressive individuals with to goal of reducing overall levels of conflict; and it is desirable that this removal does not negatively affect the cohesion of the group [28–30]. Figure 6(b) shows the ${O}_{{\rm{D}}}$–${O}_{{\rm{K}}}$ curves for the primate network, and we find a very weak trade-off. This is largely explained by the inter-layer correlations, which in contrast to the airline network, reduce the trade-off. To understand this recall figure 2(b), where we showed that Kemeny distance is small between the conflict-signaling and the grooming-huddling layer pairs, while there is no or even negative correlation between the competitive and affiliative layers. This particular structure allows us to simultaneously disrupt the competitive layers without affecting the affiliative layers. Furthermore, the competitive layers have more heterogeneous degree distributions than the affiliative layers (table 1), a property also seen in human social networks [46]. Using model networks, we showed that destroying heterogenous layers and keeping homogeneous layers intact reduces the trade-off (figure 5(c)). Indeed, comparing the FR and DPR + NLR curves shows that the degree distribution of the layers contributes to the weak trade-off, albeit less than the inter-layer correlations.

Finally, for the human social network, our goal is to reduce connectivity in two layers and keep the other three intact. For both the airline network and the primate social network, we found that the inter-layer correlations significantly affect the ${O}_{{\rm{D}}}$–${O}_{{\rm{K}}}$ trade-off; interestingly, in case of the human social network inter-layer correlations have little effect. The original trade-off curve is best approximated by the NLR, which removes correlations between layers, but preserves all structure within each layer; and all other randomizations show stronger trade-off than the original. Therefore, we conclude that the internal structure of the layers beyond the degree distribution, such as community structure, is what reduces trade-off.

In this paper, we explored the use of consensus rankings to identify effective targets for multi-objective interventions in multiplex networks. Our strategy is to calculate a centrality for each layer that is relevant to its specific function and our specific objectives, rank the nodes based on these centralities, and then using rank aggregation methods identify a consensus ranking. As an example process, we studied the degree-based node removal process, where we aimed to destroy the largest connected component in some layers, while keeping the other layers intact. We demonstrated that removing the nodes in order of consensus ranking effectively balances these competing objectives.

The advantage of our method is that it is agnostic to the specific properties of the layers and our goals, making it a widely applicable tool. However, the price of this flexibility is that methods designed for a specific system likely outperform our general approach—future work should explore this possible trade-off. In our work, we extensively investigated how inter-layer structural correlations affect the trade-off between different objectives. We have not, however, explored the scenario when the dynamics of the layers are also directly coupled. Centrality measures that take such coupling into account were only developed for multilayer networks where all layers follow qualitatively similar dynamics [21, 26]. It would be interesting to extend rank aggregation techniques to directly consider coupled dynamics.

We thank Haochen Wu for the airline network dataset. We gratefully acknowledge support from the US Army Research Office MURI Award No. W911NF-13-1-0340, DARPA Award No. W911NF-17-1-0077, and the National Institutes of Health Award No. R24-OD011136.

Many methods exist to approximate the Kemeny consensus, a recent survey compared the performance of 104 algorithms and combinations of algorithms [42]. It found that so-called local search methods provide an optimal trade-off between accuracy and run-time, meaning that algorithms with significantly longer runtime only marginally decreased the cost function provided by equation (3). Local search algorithms start from an initial ranking which can be either random or an approximation of the Kemeny consensus provided by another algorithm. Then this ranking is improved on by a series of local transformations that decrease the cost function. Here, we implement a version of local search based on the simple insert sort algorithm.

As initial ranking we use ${R}_{\mathrm{BC}}$, the output of the BC algorithm. Let node i₁ be the top ranked node in ${R}_{\mathrm{BC}}$, node i₂ the second, and so on; if there are ties in ${R}_{\mathrm{BC}}$, we randomly break them. We then construct a new ranking ${R}_{\mathrm{KA}}$ by starting from an empty ranking and iteratively adding nodes. First, we add node i₁. We then add i₂ such that it has new rank (i) above i₁, (ii) tied with i₁, or (iii) ranked below i₁. The new rank of i₂ is chosen to minimize the cost function. We repeat this step until all nodes are assigned a new rank.

We introduce a method to generate multiplex networks with N nodes, L layers and tunable Kemeny distance D_αβ between each pair of layers α, β ∈ {1,2, ..., L}. We start by generating L independent layers using any single-layer network model of choice. We then calculate the rank r_α,i of each node i in each layer α based on some centrality c_α,i. Note that we are not restricted to use the same single-layer network model or the same centrality for all layers. In the following, we describe a procedure to re-label the nodes in each layer to specify the inter-layer dependency between the ranks of nodes.

Let $P({c}_{1},\,\ldots ,\,{c}_{L})$ be the CDF of the multiplex centrality distribution, i.e. $P({c}_{1},\,\ldots ,\,{c}_{L})$ is the probability that a randomly selected node i has ${c}_{\alpha ,i}\leqslant {c}_{\alpha }$ for every α; and let P_α(c_α) be the CDF of the marginal centrality distribution of layer α. Specifying the multiplex centrality distribution would allow us to control the dependency between layers; we, however, cannot arbitrarily choose P(c₁, ..., c_L), since the marginals P_α(c_α) are determined by the properties of the individual layers. Yet, according to Sklar's theorem, we can always write the multiplex centrality distribution as

$$\begin{eqnarray}&&P({c}_{1},\,\ldots ,\,{c}_{L})=C({P}_{1}({c}_{1}),\,\ldots ,\,{P}_{L}({c}_{L})),\end{eqnarray} \tag{ B1 }$$

where C(u₁, ..., u_L) is an L-variate copula [47]. A copula $C:{[0,1]}^{L}\to [0,1]$ is the CDF of a random vector (u₁, ..., u_L) with uniform margins [48, 49]. The advantage of this representation is that it separates the marginal distributions of each layer, specified by P_α(c_α), and the interdependency structure between layers, specified by C(u₁, ..., u_L). If the marginals P_α(c_α) and P_β(c_β) are continuous, the Kemeny distance D_αβ between the two layers is the same as the Kendall's distance and it is completely defined by the copula

$$\begin{eqnarray}&&{D}_{\alpha \beta }=1-2{\int }_{0}^{1}{{\rm{d}}{u}}_{\alpha }{\int }_{0}^{1}{{\rm{d}}{u}}_{\beta }C({u}_{\alpha },{u}_{\beta })\displaystyle \frac{{\partial }^{2}C({u}_{\alpha },{u}_{\beta })}{\partial {u}_{\alpha }\partial {u}_{\beta }},\end{eqnarray} \tag{ B2 }$$

where C(u_α, u_β) is the two-dimensional marginal of the L-variate copula [49]. For example, for independent layers we have $C({u}_{1},\,\ldots ,\,{u}_{L})={\prod }_{\alpha =1}^{L}{u}_{\alpha }$, and substituting this into equation (B2), we get D = 0.5 for all layer pairs.

To re-label the nodes such that the multiplex centrality distribution P(c₁, ..., c_L) follows equation (B1), we draw a random vector u_i = (u_{1, i}, ..., u_{L, i}) from the L-variate copula for each node i, and from this we obtain a vector of ranks s_i = (s_{1, i}, ..., s_{L, i}) where

$$\begin{eqnarray}&&{s}_{\alpha ,i}=N-| \{j\in \{1,\,\ldots ,\,N\}:\,{u}_{\alpha ,i}\gt {u}_{\alpha ,j}\}| .\end{eqnarray} \tag{ B3 }$$

We re-label each node such that ${r}_{\alpha ,i}$, the rank of i in layer α, becomes equal to s_α,i. For example, if s_i = (1, 3), we re-label nodes such that the highest ranked node in layer 1 and third highest ranked node in layer 2 are labeled i.

Throughout this paper, we use the Gaussian copula

$$\begin{eqnarray}&&{C}_{\rho }({u}_{1},\,\ldots ,\,{u}_{L})={{\rm{\Phi }}}_{\rho }({{\rm{\Phi }}}^{-1}({u}_{1}),\,\ldots ,\,{{\rm{\Phi }}}^{-1}({u}_{L})),\end{eqnarray} \tag{ B4 }$$

where Φ_ρ is the CDF of the L-variate standard normal distribution with correlation matrix $\rho \in {[-1,1]}^{L\times L}$ and Φ⁻¹ is the inverse CDF of the standard normal distribution. Substituting the Gaussian copula into equation (B2), we find the relationship

$$\begin{eqnarray}&&{\rho }_{\alpha \beta }=\sin \left(\displaystyle \frac{\pi }{4}(1-{D}_{\alpha \beta })\right),\end{eqnarray} \tag{ B5 }$$

where ρ_αβ is an element of the correlation matrix ρ. Therefore, choosing the correlation matrix ρ allows us to set the Kemeny distance between layers.

Finally, note that the above procedure does not take into account ties in the ranks, i.e. we assume that ${c}_{\alpha ,i}\ne {c}_{\alpha ,j}$ for any $i\ne j$. Therefore, if there are any ties in the centralities (e.g. if degree centrality is used), we randomly break them. Furthermore, equation (B2) assumes that the marginal centrality distributions are continuous, if this is not the case, the Kemeny distance D_αβ between layers is no longer exactly provided exactly by equations (B2) and (B5); through simulations, however, we found that D_αβ is still well approximated by them.

B.1. Static scale-free model

To generate the scale-free layers, we use the static model [50]. Starting from N unconnected nodes, we assign a weight w_i = i^−α to each node i = 1, ..., N, where $\alpha \in [0,1)$. We then randomly select two nodes i and j with probability proportional to w_i and w_j, respectively, and if there is no link between nodes i and j, we connect them. We repeat this step until L links are added. The resulting network has average degree c = 2L/N and its degree distribution can be written as sum of Poisson distributions

$$\begin{eqnarray}&&p(k)=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{c}_{i}^{k}}{k!}{{\rm{e}}}^{-{c}_{i}},\end{eqnarray} \tag{ B6 }$$

where ${c}_{i}=2{{Lw}}_{i}/({\sum }_{j=0}^{N}{w}_{j})$ is the expected degree of node i. For large N, the degree distribution is approximated as

$$\begin{eqnarray}&&p(k)=\displaystyle \frac{{[c(1-\alpha )]}^{1/\alpha }}{\alpha }\displaystyle \frac{{\rm{\Gamma }}(k-1/\alpha ,c[1-\alpha ])}{{\rm{\Gamma }}(k+1)},\end{eqnarray} \tag{ B7 }$$

where Γ(z) is the gamma function and Γ(z, a) is the upper incomplete gamma function. We refer to this network as scale-free, because the tail of the distribution decays as a power-law, i.e. $p(k)\sim {k}^{-\gamma }$, where γ = 1 + 1/α.

When numerically evaluating the equations developed in section 3.1, we use equation (B6) to represent p(k).