# Inverse Percolation to Quantify Robustness in Multiplex Networks

Authors: Edwin Montes-Orozco; Roman-Anselmo Mora-Gutiérrez; Bibiana Obregón-Quintana; Sergio-G. de-los-Cobos-Silva; Eric A. Rincón-García; Pedro Lara-Velázquez; Miguel A. Gutiérrez-Andrade
Journal: Complexity
 (2020)

---

## Abstract

Inverse percolation is known as the problem of finding the minimum set of nodes whose elimination of their links causes the rupture of the network. Inverse percolation has been widely used in various studies of single-layer networks. However, the use and generalization of multiplex networks have been little considered. In this work, we propose a methodology based on inverse percolation to quantify the robustness of multiplex networks. Specifically, we present a modified version of the mathematical model for the multiplex-vertex separator problem (m-VSP). By solving the m-VSP, we can find nodes that cause the rupture of the mutually connected giant component (MCGC) and the large viable cluster (LVC) when their links are removed from the network. The methodology presented in this work was tested in a set of benchmark networks, and as case study, we present an analysis using a set of multiplex social networks modeled with information about the main characteristics of the best universities in the world and the universities in Mexico. The results show that the methodology presented in this work can work in different models and types of 2- and 3-layer multiplex networks without dividing the entire multiplex network into single-layer as some techniques described in the specific literature. Furthermore, thanks to the fact that the technique does not require the calculation of some structural measure or centrality metric, and it is easy to scale for networks of different sizes.


---

## Body

## 1. Introduction

### 1.1. Complex Networks and Robustness

A complex network is a network with nontrivial topological characteristics, which do not occur in simple networks, such as degree distributions, high local cohesiveness (measured through the clustering coefficients, etc.), community structures, and hierarchical structures [1].In complex networks, robustness is a measure of the capacity of complex systems versus specific external or internal events to preserve a relatively normal behavior. The robustness of the network is understood as its ability to keep most nodes together after a disturbance [1].Robustness in complex networks can be considered as an inverse percolation process [2] because the percolation theory models the process of placing nodes in a random n-dimensional network with some probability p and thus predict the sudden formation of a giant component (GC). Therefore, the inverse percolation process can be considered as the time when the rupture (separation) of the GC occurs.


### 1.2. Multiplex Networks

In recent years, the study of multilayer networks has been emphasized, thanks to the fact that most real systems have structures with multiple types of links or interactions between nodes [3], for example, multimodal transport systems, biological systems, social networks, and numerous modes of communication [4].On the other hand, multiplex networks are a particular class of multilayer networks, which were introduced to better model complex real-world systems [5–7]. The main characteristic of multiplex networks is that all the nodes in each layer are replicated in other layers, and there is a direct link between each replica node to denote the relationship.Formally, letGP=Gα,C∀α∈1,…,M be a multiplex network where:(i)
Gα=Xα,Eα is a monoplex network called layer α, where Xα and Eα are the set of nodes (it is important to mention that, in multiplex networks, each node belongs to all the M layers.) and the links in layer α, respectively.(ii)
G=Eαβ⊆Xα×Xβ;α,β∈1,…,M,α≠β is the set of interconnections between nodes in different layers. The elements of Ğ are called cross layers, and the elements of each Eα are called intralayer connections of GP.For these types of networks, a question that has been the subject of considerable recent research is as follows: if we work with a multilayer or multiplex network, how robust is the communication with the failure or elimination of certain nodes? In response to this question, it is known that the elimination of certain nodes in a single-layer can cause cascading failures.Furthermore, as in the case of single-layer networks, if the nodes are removed randomly, the effect on the transmission is very slight, but if the nodes to be removed are chosen carefully, they can cause the entire system to collapse.


### 1.3. Approach

In this work, we present a methodology to quantify robustness in multiplex networks based on the inverse percolation process applied to the mutually connected giant component (MCGC) and the large viable cluster (LVC). We analyze this methodology through a set of instance networks, and as a case study, we present a set of networks modeled with information about the main characteristics of 377 universities that belong to the 500 best in the world (QSRanking networks) and 60 universities that belong to the Mexican Higher Education System (MHES networks).It is worth mentioning that the QSRanking networks were modeled using the information available in the QS Intelligence Unit [8], where the main characteristics studied are academic reputation, employer reputation, college students, college appointments, international teachers, and international students.On the other hand, the MHES networks were modeled using information available in the EXECUM [9] repository of the National Autonomous University of Mexico (UNAM), where some of the main characteristics are number of professors with recognition from the National System of Researchers (SNI) of the Mexican National Council of Science and Technology (CONACYT), professors with Masters or Ph.D., registered students, accredited academic programs, publications in SCOPUS, publications in indexed journals, postgraduate programs incorporated into CONACYT’s National Registry of Quality Graduate Programs (PNPC), among others [10].The idea of modeling 2- and 3-layer multiplex networks is based on the advantage of analyzing the most important nodes based on the information of 2 or 3 years. Thus, it obtains an identification based on a multicriteria analysis that reduces the disadvantages of using only one layer.


## 1.1. Complex Networks and Robustness

A complex network is a network with nontrivial topological characteristics, which do not occur in simple networks, such as degree distributions, high local cohesiveness (measured through the clustering coefficients, etc.), community structures, and hierarchical structures [1].In complex networks, robustness is a measure of the capacity of complex systems versus specific external or internal events to preserve a relatively normal behavior. The robustness of the network is understood as its ability to keep most nodes together after a disturbance [1].Robustness in complex networks can be considered as an inverse percolation process [2] because the percolation theory models the process of placing nodes in a random n-dimensional network with some probability p and thus predict the sudden formation of a giant component (GC). Therefore, the inverse percolation process can be considered as the time when the rupture (separation) of the GC occurs.


## 1.2. Multiplex Networks

In recent years, the study of multilayer networks has been emphasized, thanks to the fact that most real systems have structures with multiple types of links or interactions between nodes [3], for example, multimodal transport systems, biological systems, social networks, and numerous modes of communication [4].On the other hand, multiplex networks are a particular class of multilayer networks, which were introduced to better model complex real-world systems [5–7]. The main characteristic of multiplex networks is that all the nodes in each layer are replicated in other layers, and there is a direct link between each replica node to denote the relationship.Formally, letGP=Gα,C∀α∈1,…,M be a multiplex network where:(i)
Gα=Xα,Eα is a monoplex network called layer α, where Xα and Eα are the set of nodes (it is important to mention that, in multiplex networks, each node belongs to all the M layers.) and the links in layer α, respectively.(ii)
G=Eαβ⊆Xα×Xβ;α,β∈1,…,M,α≠β is the set of interconnections between nodes in different layers. The elements of Ğ are called cross layers, and the elements of each Eα are called intralayer connections of GP.For these types of networks, a question that has been the subject of considerable recent research is as follows: if we work with a multilayer or multiplex network, how robust is the communication with the failure or elimination of certain nodes? In response to this question, it is known that the elimination of certain nodes in a single-layer can cause cascading failures.Furthermore, as in the case of single-layer networks, if the nodes are removed randomly, the effect on the transmission is very slight, but if the nodes to be removed are chosen carefully, they can cause the entire system to collapse.


## 1.3. Approach

In this work, we present a methodology to quantify robustness in multiplex networks based on the inverse percolation process applied to the mutually connected giant component (MCGC) and the large viable cluster (LVC). We analyze this methodology through a set of instance networks, and as a case study, we present a set of networks modeled with information about the main characteristics of 377 universities that belong to the 500 best in the world (QSRanking networks) and 60 universities that belong to the Mexican Higher Education System (MHES networks).It is worth mentioning that the QSRanking networks were modeled using the information available in the QS Intelligence Unit [8], where the main characteristics studied are academic reputation, employer reputation, college students, college appointments, international teachers, and international students.On the other hand, the MHES networks were modeled using information available in the EXECUM [9] repository of the National Autonomous University of Mexico (UNAM), where some of the main characteristics are number of professors with recognition from the National System of Researchers (SNI) of the Mexican National Council of Science and Technology (CONACYT), professors with Masters or Ph.D., registered students, accredited academic programs, publications in SCOPUS, publications in indexed journals, postgraduate programs incorporated into CONACYT’s National Registry of Quality Graduate Programs (PNPC), among others [10].The idea of modeling 2- and 3-layer multiplex networks is based on the advantage of analyzing the most important nodes based on the information of 2 or 3 years. Thus, it obtains an identification based on a multicriteria analysis that reduces the disadvantages of using only one layer.


## 2. Related Work

Currently, there is no consensus about the definition of robustness in complex networks; however, there are some works that analyze these concepts [11–14].In some works, such as [15–17], the authors define several measures to quantify robustness in single-layer networks, which are based on the identification of the most vulnerable nodes in the network, and estimate the importance of a disturbance or failure.For example, in the work of Ren et al. [18], the authors describe that node removal algorithms are inefficient in certain conditions and propose a simple and efficient link removal strategy called Hierarchical Power Iterative Normalized Cut (HPI-Ncut). On the other hand, Jianhua et al. [19], show the design of two optimal mixed-integer nonlinear programming schemes to investigate the optimal robustness of free-scale single-layer networks.On the other hand, Morone and Makse [20] present a methodology based on the mapping of the optimal percolation problem in random networks, intending to identify the minimum set of influential nodes (spreaders).In the work of Montes-Orozco et al. [21], the authors present a methodology based on inverse percolation that finds the set of nodes that disconnect the network in order to quantify the robustness in single-layer networks.For multilayer networks, Shekhtman et al. [22] present a survey about the vulnerability and resilience of systems with two or more interdependent networks; specifically, they show the consequences of a disturbance or attack on a particular node of a specific network and describe that this is the case when some nodes have an essential role in at least two networks (layers).Furthermore, there are some works about cascade failures in multilayer networks [16, 23, 24]; however, they only focus on targeted attacks and not in system failures. For example, Kleineberg et al. [24] describe how an attack can bring cascading failures inside multilayer networks, and in the study developed by Schneider et al. [16], the authors describe how the attacks to certain node into a specific layer affect the nodes connected to it in other layers.In the work of Baxter et al. [25], the authors develop a theory of sparse multiplex networks, which makes it possible to find the mutually connected giant component (MCGC) in a 2-layer multiplex network with arbitrary correlations between connections of different types. Here, the authors found that the correlations between the overlapping and nonoverlapping bonds markedly change the phase diagram of the system, leading to multiple hybrid phase transitions.Bianconi [26] posed the following question: are there metrics that can be used to design robust multilayer networks? Some authors have tried to find a measure to answer this problem [27]; however, the question remains open [4, 28]. On the other hand, in the work of Vaknin et al. [28], the authors present a study based on the propagation of attacks located in multiplex networks, where some experiments are executed using different values for the assortativity index and the percolation threshold.In the work of Cozzo et al. [4], the authors analyze the triadic relationships for the generalization of transitivity between layers and clustering coefficient in multiplex networks. Also, the authors show how the structure of the multiplex networks adds a new fundamental intransitivity parameter.Stella [29] presents a study about the strength of the lexicon in the face of the progressive failure of words with multiplex network attacks. Here, the author indicates that the MCGC is robust to many targeted attacks that, instead, dramatically disrupted the large viable cluster (LVC).Osat et al. [30] present an analysis based on optimal percolation, where the authors describe that approximating the optimal percolation solution in a multiplex network with valid solutions for monolayer networks extracted from the multiplex can have serious consequences in the characterization of the true robustness of the system.Finally, Montes-Orozco et al. [31] present a methodology to identify the spreader nodes of COVID-19 through the inverse percolation approach in 5-layer multiplex networks using the analysis of the relationship between sociocultural and economic characteristics with the number of infections and deaths caused by the COVID-19 virus in different countries.


## 3. Materials and Methods

In this section, we present the main characteristics of benchmark networks, the strategy that we used to model the MHES networks, and the development of the multiplex-vertex separator problem (m-VSP).


### 3.1. Benchmark Networks

It is important to mention that because robustness in multiplex networks is a new topic, there is no reference set to compare the results. Therefore, to verify if m-VSP works correctly for several network models, we use four 3-layer multiplex networks with different characteristics of the reference set given by De Domenico et al. [32].The networks selected for this work are CKM Physicians Innovation Multiplex Social (CKM) [33], Hepatitis C Virus Multiplex Genetic (Hepatitis C) [34], Krackhardt High Tech Multiplex Social (Krackhardt) [35], and Lazega Law Firm Multiplex Social (Lazega) [36].Next, we present Table1, which shows the number of links and nodes for each layer in each benchmark network, where the identifiers are l1 for layer 1, l2 for layer 2, and l3 for layer 3.Table 1
Main characteristics for benchmark networks.
NetworkNodesl1l2l3l1, l2 and l3CKM449498423246Hepatitis C80432106Krackhardt145792021Lazega71739972671At this point, it is important to mention that, although benchmark networks have few links, we use them to verify that our methodology works in different models of multiplex networks. Therefore, we present the numerical results for the main structural metrics calculated in all the benchmark networks.In order to facilitate the understanding of the results about the structural metrics of networks, we briefly present some fundamental concepts:(i)
Clustering coefficient [37]: it quantifies how much a node is interconnected with its neighbors, where two nodes are neighbors if there is a link that joins them. For nondirected graphs, the clustering coefficient for a node ni is calculated in the following way:(1)Ci=2E′kiki−1:E′⊂E,whereki is the degree of the node ni and E′ is the subset of links that connect any pair of nodes in the neighborhood of node ni. For this paper, we say that a node na belongs to the neighborhood of node ni if there is a link between ni and na.(ii)
Average degree [38]: it is a measure that divides the number of links by the number of nodes.(iii)
Average path length [39]: it defines the average number of steps that must be traveled through the shortest path for all possible pairs of nodes.(iv)
Diameter of the network [40]: it is given by the longest route of all the shortest routes between any pair of nodes.(v)
Density [41]: measures the connectivity of the network at the global level, and it is calculated as(2)current  connectionspotential connections,where “potential connections” (PC) are calculated as(3)PC=NN−12,withN as the total number of nodes in the network and “current connections” are the number of links that exists in the network.In Table2, we present the numerical values for each characteristic: components, diameter, path length, size of MCGC, and size of LVC. For each characteristic, there are three columns, l1, l2, and l3 that include the information of layers 1, 2, and 3, respectively.Table 2
Characteristics of benchmark multiplex networks.
NetworkComponentsDiameterPath lengthSize of MCGCSize of LVCl1l2l3l1l2l3l1l2l3l1, l2, and l3l1, l2, and l3CKM35192279143.13.43.9105104Hepatitis C2562103442221.38080Krackhardt1113341.31.632121Lazega1313541.82.21.77153The information provided in Tables1 and 2 describes the main characteristics of the selected networks.For example, according to Table1, the Hepatitis C network has 106 nodes in each layer and 80, 43, and 2 links in layers 1, 2, and 3, respectively. Besides, according to Table 2, this network has short diameters (4, 4, and 2) and short average path lengths (2, 2, and 1.3).On the other hand, in order to verify that our methodology is easily scalable for large multiplex networks, in the following section, we show the modeling of 17 2- and 3-layer multiplex networks, in which each has a greater number of links in every single layer.


### 3.2. QSRanking and MHES Networks

As mentioned above, we present a case study analysis of the QSRanking and MHES networks. For these networks, each layer represents a different year, where a node represents each university and the links between the nodes are created using the following strategy based on the Mahalanobis distance.Mahalanobis distance measures a way to determine the similarity between two variables. In contrast with the Euclidean distance, Mahalanobis distance takes into account the correlation between the random variables.Formally, the Mahalanobis distance between two random variables with the same probability distributionx and y with covariance matrix Σ is defined as(4)dmx,y=x−yTΣ−1x−y.In our case, the Mahalanobis distance can be considered as follows: based on the different characteristics, we need to measure the similarity between two universities and thus relate those that are similar. For example, assuming that, for each university, only 2 characteristics are considered (number of students and number of professors), and with these data, we can build a vectorxi=x1i,x2iT for each university i.If we consider that the first characteristic (students) has values between 100 and 200, while the second (teachers) is between 15 and 25, then we use the Euclidean distance:(5)dex1,x2=x11−x122+x21−x222.Since the differences in the first characteristic are greater than in the second, less importance is given to the second characteristic. For this reason, we prefer to add the statistics of the data to the distance measure, weighting according to its variance: the variables with less variance will have more importance than those with greater variance. In this way, we intend to equalize the importance of all the characteristics in the final result.Therefore, the expression is as follows:(6)d2x1,x2=x11−x12σ12+x21−x22σ22,where σi is the standard deviation of the component ii of the measurement vectors.Otherwise, in vectorial notation,(7)dex1,x2=x1−x2TS−1x1−x2,where S is a matrix whose elements on the diagonal are sii=σi2.However, the expression above has a problem that most of the characteristics are not independent, that is, each one depends on the other since it is more likely that when a university has more students, it will have a greater number of professors. So, to incorporate the dependency between the two variables, we can substitute the diagonal matrixS by the covariance matrix Σ:(8)dmx1,x2=x1−x2TΣ−1x1−x2, which is the Mahalanobis distance.Then, based on the above, the relationships of the universities are given by the number of characteristics in which the universities are similar, and their quantification is obtained as follows:(1)
The Mahalanobis distance [42] between each pair of universities is calculated.(2)
The median of Mahalanobis distances is calculated.(3)
For each pair of universities with a lower distance than the median, a link is added.It is essential to mention that the information provided by ExECUM is divided into two sets: aggregated and disaggregated data. We use both sets of information to analyze the MHES.In order to better understand the results shown in the following sections, we present the identifiers used for each QSRaking and MHES network:(1)
Each monoplex QSRanking and MHES network represents a different year, where the QSRanking networks have an index∗QS.(2)
For the MHES networks, an index 1 is used to represent the networks that were generated using disaggregated data. For networks generated with aggregated data, we used an index 2. For example, 2015.1 and 2015.2 are the networks that represent the relationship between the universities in 2015, according to the disaggregated and aggregated data, respectively.(3)
The multiplex networks are formed by a set of monoplex networks. For example, for MHES multiplex networks, 2013-2015.2 is the multiplex network formed by the 2013.2, 2014.2, and 2015.2 MHES monoplex networks, and for QSRanking networks, 2018-2020.QS is the multiplex network formed by 2018.QS, 2019.QS, and 2020.QS QSRanking monoplex networks.In Table3, we present the number of nodes and links for each layer.Table 3
Main characteristics for MHES and QSRanking networks.
NetworkNodesl1l2l3l1, l2, and l32013-2014.112141144—602014-2015.111441210—602015-2016.112101036—602016-2017.110361018—602013–2015.1121410221210602014–2016.1114412101036602015–2017.1121010361018602013-2014.211301154—602014-2015.211541172—602015-2016.21172920—602016-2017.29201008—602013–2015.2113011541172602014–2016.211541172920602015–2017.211729201008602018-2019.QS964811426—3772019-2020.QS1142610668—3772018-2019.QS96481142610668377As mentioned above, the QSRanking and MHES networks have a large number of links in each layer regarding the number of nodes. Therefore, these networks can help us verify how efficient our methodology is.Now, in order to see the relationships between universities as shown in Figure1, the 2-layer multiplex network is modeled with the disaggregated data for the years 2013 and 2014, and Figure 2 shows the 3-layer network modeled with the aggregated data for the years 2015, 2016, and 2017.Figure 1
Graph of the MHES network with disaggregated data for the years 2013 and 2014.Figure 2
Graph of the MHES network with nonstandardized disaggregated data for the years 2015, 2016, and 2017.Here, it is important to mention that, in each multiplex network, the nodes represent universities in the MHES, and the intralayer links represent their relationships for a specific year. The interlayer links represent the connection of replica nodes in different years (same nodes in other layers).For example, in Figure2, the links in layers l1, l2, and l3 represent the relationships of each pair of universities in the years 2015, 2016, and 2017, respectively, while the links that join l1 with l2 and l2 with l3 are the identity relationship.Since the same universities in different layers have different relationships, finding which nodes are essential to remain the multiplex networks connected is not a trivial problem.


### 3.3. Inverse Percolation through the Adaptation of Vertex Separator Problem (m-VSP)

As mentioned above, the robustness in a network can be quantified, finding the nodes that cause the rupture of the GC (inverse percolation process). In multiplex networks, we can identify the GC as the MCGC or the LVC, which are defined as [43](1)
Mutually connected giant component: each nodei is in the MCGC if it has at least one neighbor j that belongs to the MCGC and if all its replica nodes in each interdependent network are also in the MCGC.From this definition, Bianconi et al. [43] deduce that if a node i in a certain layer of a connected network of networks is in the MCGC, then all its replica nodes are also in the MCGC.(2)
Large viable component: viability can be seen as the existence of several paths, connecting any pair of nodes in each of the layers of the entire network. By definition, an LVC is a set of nodes interconnected with each other at the same time on all individual layers [25, 44] of maximal size.This makes an LVC at most the size of the intersection of the set of connected nodes at each network layer.On the other hand, the methodology presented in this work is defined as an adaptation of the methodology denoted in [21], which is based on the resolution of the vertex separator problem (VSP) [45]. Specifically, an m-VSP consists in finding a minimum set of nodes C that, being deleted from the network, it produces a disconnection of the MCGC and the LVC of the multiplex network into two connected components (A, B), in such a way that A and B are maximized. Formally, the m-VSP can be summarized as:(i)
Instance: a MCGC or a LVC of a multiplex network GP =Gα,C.(ii)
Problem: to find a partition of the MCGC or the LVC that results into three disjoint setsA, B and C, and A and B nonempty, such thatP1) There is no interlayer or intralayer links between the elements belonging to each set A, B, or CP2) A and B are maximizedP3) C is minimizedThus, if we solve the m-VSP on a multiplex network, we obtain nodes that cause the rupture of the MCGC or LVC. Therefore, we can assure that the resolution of the m-VSP is an inverse percolation process and finds the set of relevant nodes for network communication.It is important to mention that we can quantify the robustness of the networks by analyzing the percentage of nodes that belong to the separator setC of the MCGC and LVC. For example, for MCGC, a high percentage (more than 70%) of nodes in C indicates that the network is robust, while a low percentage (less than 30%) indicates that the network is not very robust, and for LVC, a high percentage (more than 60%) of nodes in C indicates that the network is robust, while a low percentage (less than 40%) indicates that the network is not very robust.


### 3.4. Simulated Annealing to Solve m-VSP

Because the original VSP is considered NP-hard [46, 47], the adaptation from multiplex networks is NP-hard too [48, 49]; therefore, we can solve the adaptation of VSP using a heuristic method.In this work, we used a simulated annealing algorithm (SA) [50] developed in C language. SA is a technique based on local search, and it is necessary to establish the structure of neighborhoods. In the development of this work, two neighborhoods were used: Hamming distance 1 (Hamming distance is defined as the number of elements that have to be changed to transform a solution into another valid solution) and Hamming distance 2. For the first type of neighborhood, it is added or deleted one node of the current solution with the following steps:(1)
A random number between 0 and 1 is generated(2)
If random≤0.5, one node is removed; otherwise, one node is added(3)
To remove an element, it is chosen randomly, giving the same probability to each element that belongs to the current solution(4)
To add an element, it is chosen randomly, giving the same probability to each element that does not belong to the current solutionFor the second neighborhood, two nodes are removed, or two nodes are added, and if it is the case that two nodes cannot be removed or added, one is deleted, and a different one is added. With this, the distance of Hamming 2 is maintained.Finally, the algorithm finished when the best solution it is not updated, when 50 temperature updates have elapsed, orTk=Tf.SA requires 4 control parameters, and in order to obtain a good performance, it is important to find the effective parameter settings for the technique. Then, using the differential evolution algorithm (DE) [51]; we obtain the following values (for more information, see [10]):(i)
Initial temperatureTi  = 1500(ii)
Final temperatureTf  = 0.001(iii)
Cooling programTkγ, with γ  = 0.95(iv)
Times that a new neighbor is generatedLk = 20


## 3.1. Benchmark Networks

It is important to mention that because robustness in multiplex networks is a new topic, there is no reference set to compare the results. Therefore, to verify if m-VSP works correctly for several network models, we use four 3-layer multiplex networks with different characteristics of the reference set given by De Domenico et al. [32].The networks selected for this work are CKM Physicians Innovation Multiplex Social (CKM) [33], Hepatitis C Virus Multiplex Genetic (Hepatitis C) [34], Krackhardt High Tech Multiplex Social (Krackhardt) [35], and Lazega Law Firm Multiplex Social (Lazega) [36].Next, we present Table1, which shows the number of links and nodes for each layer in each benchmark network, where the identifiers are l1 for layer 1, l2 for layer 2, and l3 for layer 3.Table 1
Main characteristics for benchmark networks.
NetworkNodesl1l2l3l1, l2 and l3CKM449498423246Hepatitis C80432106Krackhardt145792021Lazega71739972671At this point, it is important to mention that, although benchmark networks have few links, we use them to verify that our methodology works in different models of multiplex networks. Therefore, we present the numerical results for the main structural metrics calculated in all the benchmark networks.In order to facilitate the understanding of the results about the structural metrics of networks, we briefly present some fundamental concepts:(i)
Clustering coefficient [37]: it quantifies how much a node is interconnected with its neighbors, where two nodes are neighbors if there is a link that joins them. For nondirected graphs, the clustering coefficient for a node ni is calculated in the following way:(1)Ci=2E′kiki−1:E′⊂E,whereki is the degree of the node ni and E′ is the subset of links that connect any pair of nodes in the neighborhood of node ni. For this paper, we say that a node na belongs to the neighborhood of node ni if there is a link between ni and na.(ii)
Average degree [38]: it is a measure that divides the number of links by the number of nodes.(iii)
Average path length [39]: it defines the average number of steps that must be traveled through the shortest path for all possible pairs of nodes.(iv)
Diameter of the network [40]: it is given by the longest route of all the shortest routes between any pair of nodes.(v)
Density [41]: measures the connectivity of the network at the global level, and it is calculated as(2)current  connectionspotential connections,where “potential connections” (PC) are calculated as(3)PC=NN−12,withN as the total number of nodes in the network and “current connections” are the number of links that exists in the network.In Table2, we present the numerical values for each characteristic: components, diameter, path length, size of MCGC, and size of LVC. For each characteristic, there are three columns, l1, l2, and l3 that include the information of layers 1, 2, and 3, respectively.Table 2
Characteristics of benchmark multiplex networks.
NetworkComponentsDiameterPath lengthSize of MCGCSize of LVCl1l2l3l1l2l3l1l2l3l1, l2, and l3l1, l2, and l3CKM35192279143.13.43.9105104Hepatitis C2562103442221.38080Krackhardt1113341.31.632121Lazega1313541.82.21.77153The information provided in Tables1 and 2 describes the main characteristics of the selected networks.For example, according to Table1, the Hepatitis C network has 106 nodes in each layer and 80, 43, and 2 links in layers 1, 2, and 3, respectively. Besides, according to Table 2, this network has short diameters (4, 4, and 2) and short average path lengths (2, 2, and 1.3).On the other hand, in order to verify that our methodology is easily scalable for large multiplex networks, in the following section, we show the modeling of 17 2- and 3-layer multiplex networks, in which each has a greater number of links in every single layer.


## 3.2. QSRanking and MHES Networks

As mentioned above, we present a case study analysis of the QSRanking and MHES networks. For these networks, each layer represents a different year, where a node represents each university and the links between the nodes are created using the following strategy based on the Mahalanobis distance.Mahalanobis distance measures a way to determine the similarity between two variables. In contrast with the Euclidean distance, Mahalanobis distance takes into account the correlation between the random variables.Formally, the Mahalanobis distance between two random variables with the same probability distributionx and y with covariance matrix Σ is defined as(4)dmx,y=x−yTΣ−1x−y.In our case, the Mahalanobis distance can be considered as follows: based on the different characteristics, we need to measure the similarity between two universities and thus relate those that are similar. For example, assuming that, for each university, only 2 characteristics are considered (number of students and number of professors), and with these data, we can build a vectorxi=x1i,x2iT for each university i.If we consider that the first characteristic (students) has values between 100 and 200, while the second (teachers) is between 15 and 25, then we use the Euclidean distance:(5)dex1,x2=x11−x122+x21−x222.Since the differences in the first characteristic are greater than in the second, less importance is given to the second characteristic. For this reason, we prefer to add the statistics of the data to the distance measure, weighting according to its variance: the variables with less variance will have more importance than those with greater variance. In this way, we intend to equalize the importance of all the characteristics in the final result.Therefore, the expression is as follows:(6)d2x1,x2=x11−x12σ12+x21−x22σ22,where σi is the standard deviation of the component ii of the measurement vectors.Otherwise, in vectorial notation,(7)dex1,x2=x1−x2TS−1x1−x2,where S is a matrix whose elements on the diagonal are sii=σi2.However, the expression above has a problem that most of the characteristics are not independent, that is, each one depends on the other since it is more likely that when a university has more students, it will have a greater number of professors. So, to incorporate the dependency between the two variables, we can substitute the diagonal matrixS by the covariance matrix Σ:(8)dmx1,x2=x1−x2TΣ−1x1−x2, which is the Mahalanobis distance.Then, based on the above, the relationships of the universities are given by the number of characteristics in which the universities are similar, and their quantification is obtained as follows:(1)
The Mahalanobis distance [42] between each pair of universities is calculated.(2)
The median of Mahalanobis distances is calculated.(3)
For each pair of universities with a lower distance than the median, a link is added.It is essential to mention that the information provided by ExECUM is divided into two sets: aggregated and disaggregated data. We use both sets of information to analyze the MHES.In order to better understand the results shown in the following sections, we present the identifiers used for each QSRaking and MHES network:(1)
Each monoplex QSRanking and MHES network represents a different year, where the QSRanking networks have an index∗QS.(2)
For the MHES networks, an index 1 is used to represent the networks that were generated using disaggregated data. For networks generated with aggregated data, we used an index 2. For example, 2015.1 and 2015.2 are the networks that represent the relationship between the universities in 2015, according to the disaggregated and aggregated data, respectively.(3)
The multiplex networks are formed by a set of monoplex networks. For example, for MHES multiplex networks, 2013-2015.2 is the multiplex network formed by the 2013.2, 2014.2, and 2015.2 MHES monoplex networks, and for QSRanking networks, 2018-2020.QS is the multiplex network formed by 2018.QS, 2019.QS, and 2020.QS QSRanking monoplex networks.In Table3, we present the number of nodes and links for each layer.Table 3
Main characteristics for MHES and QSRanking networks.
NetworkNodesl1l2l3l1, l2, and l32013-2014.112141144—602014-2015.111441210—602015-2016.112101036—602016-2017.110361018—602013–2015.1121410221210602014–2016.1114412101036602015–2017.1121010361018602013-2014.211301154—602014-2015.211541172—602015-2016.21172920—602016-2017.29201008—602013–2015.2113011541172602014–2016.211541172920602015–2017.211729201008602018-2019.QS964811426—3772019-2020.QS1142610668—3772018-2019.QS96481142610668377As mentioned above, the QSRanking and MHES networks have a large number of links in each layer regarding the number of nodes. Therefore, these networks can help us verify how efficient our methodology is.Now, in order to see the relationships between universities as shown in Figure1, the 2-layer multiplex network is modeled with the disaggregated data for the years 2013 and 2014, and Figure 2 shows the 3-layer network modeled with the aggregated data for the years 2015, 2016, and 2017.Figure 1
Graph of the MHES network with disaggregated data for the years 2013 and 2014.Figure 2
Graph of the MHES network with nonstandardized disaggregated data for the years 2015, 2016, and 2017.Here, it is important to mention that, in each multiplex network, the nodes represent universities in the MHES, and the intralayer links represent their relationships for a specific year. The interlayer links represent the connection of replica nodes in different years (same nodes in other layers).For example, in Figure2, the links in layers l1, l2, and l3 represent the relationships of each pair of universities in the years 2015, 2016, and 2017, respectively, while the links that join l1 with l2 and l2 with l3 are the identity relationship.Since the same universities in different layers have different relationships, finding which nodes are essential to remain the multiplex networks connected is not a trivial problem.


## 3.3. Inverse Percolation through the Adaptation of Vertex Separator Problem (m-VSP)

As mentioned above, the robustness in a network can be quantified, finding the nodes that cause the rupture of the GC (inverse percolation process). In multiplex networks, we can identify the GC as the MCGC or the LVC, which are defined as [43](1)
Mutually connected giant component: each nodei is in the MCGC if it has at least one neighbor j that belongs to the MCGC and if all its replica nodes in each interdependent network are also in the MCGC.From this definition, Bianconi et al. [43] deduce that if a node i in a certain layer of a connected network of networks is in the MCGC, then all its replica nodes are also in the MCGC.(2)
Large viable component: viability can be seen as the existence of several paths, connecting any pair of nodes in each of the layers of the entire network. By definition, an LVC is a set of nodes interconnected with each other at the same time on all individual layers [25, 44] of maximal size.This makes an LVC at most the size of the intersection of the set of connected nodes at each network layer.On the other hand, the methodology presented in this work is defined as an adaptation of the methodology denoted in [21], which is based on the resolution of the vertex separator problem (VSP) [45]. Specifically, an m-VSP consists in finding a minimum set of nodes C that, being deleted from the network, it produces a disconnection of the MCGC and the LVC of the multiplex network into two connected components (A, B), in such a way that A and B are maximized. Formally, the m-VSP can be summarized as:(i)
Instance: a MCGC or a LVC of a multiplex network GP =Gα,C.(ii)
Problem: to find a partition of the MCGC or the LVC that results into three disjoint setsA, B and C, and A and B nonempty, such thatP1) There is no interlayer or intralayer links between the elements belonging to each set A, B, or CP2) A and B are maximizedP3) C is minimizedThus, if we solve the m-VSP on a multiplex network, we obtain nodes that cause the rupture of the MCGC or LVC. Therefore, we can assure that the resolution of the m-VSP is an inverse percolation process and finds the set of relevant nodes for network communication.It is important to mention that we can quantify the robustness of the networks by analyzing the percentage of nodes that belong to the separator setC of the MCGC and LVC. For example, for MCGC, a high percentage (more than 70%) of nodes in C indicates that the network is robust, while a low percentage (less than 30%) indicates that the network is not very robust, and for LVC, a high percentage (more than 60%) of nodes in C indicates that the network is robust, while a low percentage (less than 40%) indicates that the network is not very robust.


## 3.4. Simulated Annealing to Solve m-VSP

Because the original VSP is considered NP-hard [46, 47], the adaptation from multiplex networks is NP-hard too [48, 49]; therefore, we can solve the adaptation of VSP using a heuristic method.In this work, we used a simulated annealing algorithm (SA) [50] developed in C language. SA is a technique based on local search, and it is necessary to establish the structure of neighborhoods. In the development of this work, two neighborhoods were used: Hamming distance 1 (Hamming distance is defined as the number of elements that have to be changed to transform a solution into another valid solution) and Hamming distance 2. For the first type of neighborhood, it is added or deleted one node of the current solution with the following steps:(1)
A random number between 0 and 1 is generated(2)
If random≤0.5, one node is removed; otherwise, one node is added(3)
To remove an element, it is chosen randomly, giving the same probability to each element that belongs to the current solution(4)
To add an element, it is chosen randomly, giving the same probability to each element that does not belong to the current solutionFor the second neighborhood, two nodes are removed, or two nodes are added, and if it is the case that two nodes cannot be removed or added, one is deleted, and a different one is added. With this, the distance of Hamming 2 is maintained.Finally, the algorithm finished when the best solution it is not updated, when 50 temperature updates have elapsed, orTk=Tf.SA requires 4 control parameters, and in order to obtain a good performance, it is important to find the effective parameter settings for the technique. Then, using the differential evolution algorithm (DE) [51]; we obtain the following values (for more information, see [10]):(i)
Initial temperatureTi  = 1500(ii)
Final temperatureTf  = 0.001(iii)
Cooling programTkγ, with γ  = 0.95(iv)
Times that a new neighbor is generatedLk = 20


## 4. Results and Discussion

In this section, we show the numerical results for the quantification of robustness in multiplex networks based on the m-VSP.


### 4.1. Benchmark Networks

To verify if the proposed methodology can be applied in different types and models of networks, in this section, we present an analysis of quantification of robustness in different 3-layer multiplex networks.At this point, it is important to mention that the robustness analysis in multiplex networks is an emerging research area, and we could not find an appropriate benchmark set in the specific literature.Therefore, we decided to use a set of benchmark networks (Tables1 and 2) that have been used for different purposes, such as the identification of communities and the analysis of their structural properties [52–54].For example, there are some methodologies proposed in the specific literature, such as [16, 22–24], in which, the authors do not carry out analysis or experiments in benchmark networks since they focus only on their own case studies.In Table4, we present the numerical values obtained by applying the m-VSP to the MCGC, while in Table 5, we show the numerical values obtained by applying the m-VSP to the LVC. In these tables, in the first column, we show the identifier for each network; in the second column, the values for C; in the third column, the values for A and B; and in the fourth column, the percentage of nodes in MCGC that belong to C.Table 4
Numerical results for quantification of robustness in benchmark networks (MCGC).
NetworkCA and B% of nodes of MCGC inCCKM5023 and 2847.61Hepatitis1—1.25Krackhardt66 and 128.57Lazega3921 and 254.92Table 5
Numerical results for quantification of robustness in benchmark networks (LVC).
NetworkCA and B% of nodes of LVC inCCKM1919 and 1718.26Hepatitis1—1.25Krackhardt66 and 128.57Lazega2119 and 12%29.57It is important to remember that ifC has a value greater than 70% of the total nodes of the MCGC or the 60% of the LVC, it is considered robust. On the other hand, when C has a value lower than 30% for MCGC or 40% for LVC, we can affirm that the network is not very robust.Based on the information of Table4, we can see that, from Krackhardt and Lazega, the percentage of nodes of the LVC in set C are moderately robust since each network has a C of 28.57% and 54.92%, respectively.In the specific case of the Hepatitis network, there is only one node in the separator set; and the subsetsA and B cannot be formed since each node that belongs to the network remains as an isolated component, and it is very prone to collapse.Based on the information of Table5, we can see that, from Krackhardt and Lazega, the percentage of nodes of the LVC in set C are not robust since each network has a C of 28% and 29.57%, respectively.As in the case of the analysis using the MCGC, for the Hepatitis network, there is only one node in the separator set, and the subsetsA and B cannot be formed since each node that belongs to the network remains as an isolated component. These results are produced because the network has a star structure, and it is very prone to collapse. On the other hand, the CKM network has a separator set C of 19 nodes that are around 18% of the total nodes. Therefore, both networks are not very robust.Therefore, based on the above, we can see that the m-VSP can find the nodes important to maintaining the communication and the correct operation of the systems.


### 4.2. QSRanking and MHES Networks

Based on the previous section results, we can see that the proposed methodology works for different types of network models (applied to its MCGC and LVC). Therefore, we present the structural metrics analysis and quantification of robustness for the QSRanking and MHES networks.


#### 4.2.1. Structural Metrics

In Table6, we present the numerical results for the structural metrics in the monoplex HESM networks modeled using the disaggregated and aggregated data and in the QSRanking monoplex networks, where we can see that, from QSRanking networks, the average degree is around a quarter part of the total number of nodes. For the MHES, the average degree value for networks modeled with aggregated data is smaller than for networks modeled using disaggregated data; however, the maximum degree in most networks (except 2016.2) is around half of the total network nodes.Table 6
Numerical values for structural metrics in QSRanking MHES monoplex networks.
NetworkAverage degreeMax. degreeClusteringPath length2013.118.276300.6321.8822014.118.508300.6111.8852015.118.4300.68881.8742016.116.103300.5741.9662017.116.237300.5742.1132013.215.733300.6182.1062014.214.533300.6292.3252015.214.678300.6342.6352016.25.786130.7163.7572017.215.017270.7222.562018.QS91.3581090.8581.132019.QS198.7131100.8631.282020.QS96.1081090.8651.26We can see that all networks have a high clustering coefficient (between 0.6 and 0.9) and have a small path length between any pair of nodes (between 2 and 3). Therefore, we can deduce that RankingQS and MHES networks follow the small-world model [28], which implies that they comply with the main characteristics of social networks (SN) [36].


#### 4.2.2. Quantification of Robustness

In this section, we present the numerical results obtained by solving the m-VSP to the MCGC and the LVC for each QSRanking and MHES network; Table7 presents the cardinality value for sets A, B, and C and the percentage of nodes of MCGC that belong to C. Here, it is important to mention that, in the QSRanking and MHES networks, the MCGC is formed by all the nodes of the networks.Table 7
Numerical results for m-VSP for QSRanking and MHES networks (MCGC).
NetworkCA and B% of nodes of MCGC inC2013-2014.11428 and 1823.32014-2015.11324 and 2321.72015-2016.11527 and 18252016-2017.11925 and 1631.72013–2015.12418 and 18402014–2016.12721 and 12452015–2017.12419 and 17402013-2014.21827 and 15302014-2015.21427 and 1933.32015-2016.21525 and 20252016-2017.21929 and 1631.72013–2015.22619 and 1543.32014–2016.22420 and 16402015–2017.22618 and 1643.32018-2019.QS112137 and 128342019-2020.QS133137 and 10735.32018–2020.QS174108 and 9546.1According to the results shown in Table7, we can assure that QSRanking and MHES networks are moderately robust; since for all QSRanking networks, C is between 34% and 46%, and for all MHES networks, C is between 20% and 46% of the total nodes belonging to the MCGC. For example, for the networks 2013.1-2014.1 and 2018-2020.QS, we have 23% and 46% of the nodes in the separator set, respectively.Now, we present the numerical results obtained by solving the m-VSP to the LVC for each QSRanking and MHES network; Table8 presents the cardinality value for sets A, B, and C and the percentage of nodes of LVC that belong to C. Here, it is important to mention that the size of the LVC for each network is provided in Table 8.Table 8
Size of LVC for each MHES and QSRanking network.
NetworkSize of LVCNetworkSize of LVC2013-2014.1532015-2016.2492014-2015.1472016-2017.2502015-2016.1492013-2015.2512016-2017.1512014-2016.2462013–2015.1532015-2017.2522014–2016.1472018.QS3192015–2017.1492019.QS3242013-2014.2512020.QS3282014-2015.247——Based on Table8, we can see that, for the MHES and QSRanking networks, the LVC is smaller compared to the MCGC. This is because the MCGC is composed of the union, and the LVC is composed at most by the union of all the giant components of every single layer.According to the results shown in Table9, we can assure that QSRanking and MHES networks are not robust; since for all QSRanking networks, C is between 28 and 33.8% and for all MHES networks, C is between 16.9% and 36.7% of the total nodes belonging to the LVC. For example, for the networks 2016.1-2017.1 and 2019-2020.QS, we have 31.3% and 35.2% of the nodes in the separator set, respectively.Table 9
Numerical results for m-VSP for QSRanking and MHES networks (LVC).
NetworkCA and B% of nodes of LVC inC2013.1-2014.1926 and 2516.92014.1-2015.11126 and 2323.42015.1-2016.11228 and 20242016.1-2017.11623 and 2131.32013.1–2015.11822 and 2033.92014.1–2016.11225 and 2325.52015.1–2017.11822 and 2036.72013.2-2014.21723 and 2033.32014.2-2015.21327 and 1927.62015.2-2016.21225 and 2324.42016.2-2017.21329 and 18262013.2–2015.21822 and 2035.32014.2–2016.21622 and 2234.82015.2–2017.21623 and 2130.72018.QS-2019.QS108142 and 12733.82019.QS-2020.QS114137 and 12635.22018.QS–2020.QS92149 and 13628Furthermore, we can verify that the methodology proposed in this work can find the set of nodes that cause the rupture of the MCGC and the LVC for different models of 2-layer and 3-layer multiplex networks since the sum ofA, B, and C is equivalent to the total number of nodes that belong to the MCGC or the LVC according to the case.It is important to mention that, for MHES networks, there is no difference in the networks modeled with both disaggregated and aggregated data, so they have the same numerical results when m-VSP is applied.The above means that both types of networks are equivalent and imply that it is possible to use a reduced set of information to analyze the behavior of the MHES; therefore, we can assure that the main characteristics of the complex system are properly structured when we modeled the QSRanking and MHES networks.Then, we can verify that the methodology proposed in this work can find the set of nodes that cause the rupture of the MCGC or the LVC for different models of 2-layer and 3-layer multiplex networks.In order to analyze the elements that cause the rupture of the MCGC and the LVC in the QSRanking and MHES networks, based on the results presented in [8, 28] about the classification of universities, we perform an analysis.In Table10, we show the percentage of the 377 universities that we consider to model the QSRanking networks that are located between positions 1 and 100; 101 and 200; 201 and 300; 301 and 400; and 401 and 500 in the main university rankings (like the QS ranking).Table 10
Percentage of universities in 100-universities ranges (QSRanking).
Network1–100 (%)101–200 (%)201–300 (%)301–400401–500 (%)2018-2019.QS20.729.328.415.26.42019-2020.QS18.429.328.717.5%6.12018–2020.QS19.527.53014%9Based on the information presented in Table10, we can see that most of the universities that cause the rupture of the network are between positions 100 and 300 of the QS ranking, and their main characteristics are as follows: average reputation, high employer reputation, average number of citations, few international faculties, and a greater number of national students.Now, in Table11, we show the percentage of universities with high, average, and low values in the university rankings that cause the disconnection of the MHES.Table 11
Percentage of universities with high, average, and low-ranking value (MHES).
NetworkHigh ranking (%)Average ranking (%)Low ranking (%)2013.1-2014.135.728.635.72014.1-2015.130.423.346.32015.1-2016.14033.326.62016.1-2017.136.626.636.62013.1–2015.133.32937.52014.1–2016.129.640.729.62015.1–2017.1332541.62013.1-2014.233.322.244.42014.1-2015.235.728.635.72015.1-2016.24033.326.72016.1-2017.242.126.331.62013.1–2015.234.630.734.62014.1–2016.233.329.134.62015.1–2017.238.430.830.8Based on the results of Table11, we can see that, for all MHES networks, the percentage of universities with high-, average-, and low-ranking values is divided equally, that is, a similar amount of the three types of universities are necessary to cause the rupture of the LVC.From a social network perspective, we can deduce that the links between the different types of universities are essential for the correct performance and development of the QSRanking and MHES networks.Finally, it is important to mention that the proposed methodology obtains a set of nodes that cause the rupture of the MCGC and the LVC without using structural network information or centrality measures.


## 4.1. Benchmark Networks

To verify if the proposed methodology can be applied in different types and models of networks, in this section, we present an analysis of quantification of robustness in different 3-layer multiplex networks.At this point, it is important to mention that the robustness analysis in multiplex networks is an emerging research area, and we could not find an appropriate benchmark set in the specific literature.Therefore, we decided to use a set of benchmark networks (Tables1 and 2) that have been used for different purposes, such as the identification of communities and the analysis of their structural properties [52–54].For example, there are some methodologies proposed in the specific literature, such as [16, 22–24], in which, the authors do not carry out analysis or experiments in benchmark networks since they focus only on their own case studies.In Table4, we present the numerical values obtained by applying the m-VSP to the MCGC, while in Table 5, we show the numerical values obtained by applying the m-VSP to the LVC. In these tables, in the first column, we show the identifier for each network; in the second column, the values for C; in the third column, the values for A and B; and in the fourth column, the percentage of nodes in MCGC that belong to C.Table 4
Numerical results for quantification of robustness in benchmark networks (MCGC).
NetworkCA and B% of nodes of MCGC inCCKM5023 and 2847.61Hepatitis1—1.25Krackhardt66 and 128.57Lazega3921 and 254.92Table 5
Numerical results for quantification of robustness in benchmark networks (LVC).
NetworkCA and B% of nodes of LVC inCCKM1919 and 1718.26Hepatitis1—1.25Krackhardt66 and 128.57Lazega2119 and 12%29.57It is important to remember that ifC has a value greater than 70% of the total nodes of the MCGC or the 60% of the LVC, it is considered robust. On the other hand, when C has a value lower than 30% for MCGC or 40% for LVC, we can affirm that the network is not very robust.Based on the information of Table4, we can see that, from Krackhardt and Lazega, the percentage of nodes of the LVC in set C are moderately robust since each network has a C of 28.57% and 54.92%, respectively.In the specific case of the Hepatitis network, there is only one node in the separator set; and the subsetsA and B cannot be formed since each node that belongs to the network remains as an isolated component, and it is very prone to collapse.Based on the information of Table5, we can see that, from Krackhardt and Lazega, the percentage of nodes of the LVC in set C are not robust since each network has a C of 28% and 29.57%, respectively.As in the case of the analysis using the MCGC, for the Hepatitis network, there is only one node in the separator set, and the subsetsA and B cannot be formed since each node that belongs to the network remains as an isolated component. These results are produced because the network has a star structure, and it is very prone to collapse. On the other hand, the CKM network has a separator set C of 19 nodes that are around 18% of the total nodes. Therefore, both networks are not very robust.Therefore, based on the above, we can see that the m-VSP can find the nodes important to maintaining the communication and the correct operation of the systems.


## 4.2. QSRanking and MHES Networks

Based on the previous section results, we can see that the proposed methodology works for different types of network models (applied to its MCGC and LVC). Therefore, we present the structural metrics analysis and quantification of robustness for the QSRanking and MHES networks.


### 4.2.1. Structural Metrics

In Table6, we present the numerical results for the structural metrics in the monoplex HESM networks modeled using the disaggregated and aggregated data and in the QSRanking monoplex networks, where we can see that, from QSRanking networks, the average degree is around a quarter part of the total number of nodes. For the MHES, the average degree value for networks modeled with aggregated data is smaller than for networks modeled using disaggregated data; however, the maximum degree in most networks (except 2016.2) is around half of the total network nodes.Table 6
Numerical values for structural metrics in QSRanking MHES monoplex networks.
NetworkAverage degreeMax. degreeClusteringPath length2013.118.276300.6321.8822014.118.508300.6111.8852015.118.4300.68881.8742016.116.103300.5741.9662017.116.237300.5742.1132013.215.733300.6182.1062014.214.533300.6292.3252015.214.678300.6342.6352016.25.786130.7163.7572017.215.017270.7222.562018.QS91.3581090.8581.132019.QS198.7131100.8631.282020.QS96.1081090.8651.26We can see that all networks have a high clustering coefficient (between 0.6 and 0.9) and have a small path length between any pair of nodes (between 2 and 3). Therefore, we can deduce that RankingQS and MHES networks follow the small-world model [28], which implies that they comply with the main characteristics of social networks (SN) [36].


### 4.2.2. Quantification of Robustness

In this section, we present the numerical results obtained by solving the m-VSP to the MCGC and the LVC for each QSRanking and MHES network; Table7 presents the cardinality value for sets A, B, and C and the percentage of nodes of MCGC that belong to C. Here, it is important to mention that, in the QSRanking and MHES networks, the MCGC is formed by all the nodes of the networks.Table 7
Numerical results for m-VSP for QSRanking and MHES networks (MCGC).
NetworkCA and B% of nodes of MCGC inC2013-2014.11428 and 1823.32014-2015.11324 and 2321.72015-2016.11527 and 18252016-2017.11925 and 1631.72013–2015.12418 and 18402014–2016.12721 and 12452015–2017.12419 and 17402013-2014.21827 and 15302014-2015.21427 and 1933.32015-2016.21525 and 20252016-2017.21929 and 1631.72013–2015.22619 and 1543.32014–2016.22420 and 16402015–2017.22618 and 1643.32018-2019.QS112137 and 128342019-2020.QS133137 and 10735.32018–2020.QS174108 and 9546.1According to the results shown in Table7, we can assure that QSRanking and MHES networks are moderately robust; since for all QSRanking networks, C is between 34% and 46%, and for all MHES networks, C is between 20% and 46% of the total nodes belonging to the MCGC. For example, for the networks 2013.1-2014.1 and 2018-2020.QS, we have 23% and 46% of the nodes in the separator set, respectively.Now, we present the numerical results obtained by solving the m-VSP to the LVC for each QSRanking and MHES network; Table8 presents the cardinality value for sets A, B, and C and the percentage of nodes of LVC that belong to C. Here, it is important to mention that the size of the LVC for each network is provided in Table 8.Table 8
Size of LVC for each MHES and QSRanking network.
NetworkSize of LVCNetworkSize of LVC2013-2014.1532015-2016.2492014-2015.1472016-2017.2502015-2016.1492013-2015.2512016-2017.1512014-2016.2462013–2015.1532015-2017.2522014–2016.1472018.QS3192015–2017.1492019.QS3242013-2014.2512020.QS3282014-2015.247——Based on Table8, we can see that, for the MHES and QSRanking networks, the LVC is smaller compared to the MCGC. This is because the MCGC is composed of the union, and the LVC is composed at most by the union of all the giant components of every single layer.According to the results shown in Table9, we can assure that QSRanking and MHES networks are not robust; since for all QSRanking networks, C is between 28 and 33.8% and for all MHES networks, C is between 16.9% and 36.7% of the total nodes belonging to the LVC. For example, for the networks 2016.1-2017.1 and 2019-2020.QS, we have 31.3% and 35.2% of the nodes in the separator set, respectively.Table 9
Numerical results for m-VSP for QSRanking and MHES networks (LVC).
NetworkCA and B% of nodes of LVC inC2013.1-2014.1926 and 2516.92014.1-2015.11126 and 2323.42015.1-2016.11228 and 20242016.1-2017.11623 and 2131.32013.1–2015.11822 and 2033.92014.1–2016.11225 and 2325.52015.1–2017.11822 and 2036.72013.2-2014.21723 and 2033.32014.2-2015.21327 and 1927.62015.2-2016.21225 and 2324.42016.2-2017.21329 and 18262013.2–2015.21822 and 2035.32014.2–2016.21622 and 2234.82015.2–2017.21623 and 2130.72018.QS-2019.QS108142 and 12733.82019.QS-2020.QS114137 and 12635.22018.QS–2020.QS92149 and 13628Furthermore, we can verify that the methodology proposed in this work can find the set of nodes that cause the rupture of the MCGC and the LVC for different models of 2-layer and 3-layer multiplex networks since the sum ofA, B, and C is equivalent to the total number of nodes that belong to the MCGC or the LVC according to the case.It is important to mention that, for MHES networks, there is no difference in the networks modeled with both disaggregated and aggregated data, so they have the same numerical results when m-VSP is applied.The above means that both types of networks are equivalent and imply that it is possible to use a reduced set of information to analyze the behavior of the MHES; therefore, we can assure that the main characteristics of the complex system are properly structured when we modeled the QSRanking and MHES networks.Then, we can verify that the methodology proposed in this work can find the set of nodes that cause the rupture of the MCGC or the LVC for different models of 2-layer and 3-layer multiplex networks.In order to analyze the elements that cause the rupture of the MCGC and the LVC in the QSRanking and MHES networks, based on the results presented in [8, 28] about the classification of universities, we perform an analysis.In Table10, we show the percentage of the 377 universities that we consider to model the QSRanking networks that are located between positions 1 and 100; 101 and 200; 201 and 300; 301 and 400; and 401 and 500 in the main university rankings (like the QS ranking).Table 10
Percentage of universities in 100-universities ranges (QSRanking).
Network1–100 (%)101–200 (%)201–300 (%)301–400401–500 (%)2018-2019.QS20.729.328.415.26.42019-2020.QS18.429.328.717.5%6.12018–2020.QS19.527.53014%9Based on the information presented in Table10, we can see that most of the universities that cause the rupture of the network are between positions 100 and 300 of the QS ranking, and their main characteristics are as follows: average reputation, high employer reputation, average number of citations, few international faculties, and a greater number of national students.Now, in Table11, we show the percentage of universities with high, average, and low values in the university rankings that cause the disconnection of the MHES.Table 11
Percentage of universities with high, average, and low-ranking value (MHES).
NetworkHigh ranking (%)Average ranking (%)Low ranking (%)2013.1-2014.135.728.635.72014.1-2015.130.423.346.32015.1-2016.14033.326.62016.1-2017.136.626.636.62013.1–2015.133.32937.52014.1–2016.129.640.729.62015.1–2017.1332541.62013.1-2014.233.322.244.42014.1-2015.235.728.635.72015.1-2016.24033.326.72016.1-2017.242.126.331.62013.1–2015.234.630.734.62014.1–2016.233.329.134.62015.1–2017.238.430.830.8Based on the results of Table11, we can see that, for all MHES networks, the percentage of universities with high-, average-, and low-ranking values is divided equally, that is, a similar amount of the three types of universities are necessary to cause the rupture of the LVC.From a social network perspective, we can deduce that the links between the different types of universities are essential for the correct performance and development of the QSRanking and MHES networks.Finally, it is important to mention that the proposed methodology obtains a set of nodes that cause the rupture of the MCGC and the LVC without using structural network information or centrality measures.


## 4.2.1. Structural Metrics

In Table6, we present the numerical results for the structural metrics in the monoplex HESM networks modeled using the disaggregated and aggregated data and in the QSRanking monoplex networks, where we can see that, from QSRanking networks, the average degree is around a quarter part of the total number of nodes. For the MHES, the average degree value for networks modeled with aggregated data is smaller than for networks modeled using disaggregated data; however, the maximum degree in most networks (except 2016.2) is around half of the total network nodes.Table 6
Numerical values for structural metrics in QSRanking MHES monoplex networks.
NetworkAverage degreeMax. degreeClusteringPath length2013.118.276300.6321.8822014.118.508300.6111.8852015.118.4300.68881.8742016.116.103300.5741.9662017.116.237300.5742.1132013.215.733300.6182.1062014.214.533300.6292.3252015.214.678300.6342.6352016.25.786130.7163.7572017.215.017270.7222.562018.QS91.3581090.8581.132019.QS198.7131100.8631.282020.QS96.1081090.8651.26We can see that all networks have a high clustering coefficient (between 0.6 and 0.9) and have a small path length between any pair of nodes (between 2 and 3). Therefore, we can deduce that RankingQS and MHES networks follow the small-world model [28], which implies that they comply with the main characteristics of social networks (SN) [36].


## 4.2.2. Quantification of Robustness

In this section, we present the numerical results obtained by solving the m-VSP to the MCGC and the LVC for each QSRanking and MHES network; Table7 presents the cardinality value for sets A, B, and C and the percentage of nodes of MCGC that belong to C. Here, it is important to mention that, in the QSRanking and MHES networks, the MCGC is formed by all the nodes of the networks.Table 7
Numerical results for m-VSP for QSRanking and MHES networks (MCGC).
NetworkCA and B% of nodes of MCGC inC2013-2014.11428 and 1823.32014-2015.11324 and 2321.72015-2016.11527 and 18252016-2017.11925 and 1631.72013–2015.12418 and 18402014–2016.12721 and 12452015–2017.12419 and 17402013-2014.21827 and 15302014-2015.21427 and 1933.32015-2016.21525 and 20252016-2017.21929 and 1631.72013–2015.22619 and 1543.32014–2016.22420 and 16402015–2017.22618 and 1643.32018-2019.QS112137 and 128342019-2020.QS133137 and 10735.32018–2020.QS174108 and 9546.1According to the results shown in Table7, we can assure that QSRanking and MHES networks are moderately robust; since for all QSRanking networks, C is between 34% and 46%, and for all MHES networks, C is between 20% and 46% of the total nodes belonging to the MCGC. For example, for the networks 2013.1-2014.1 and 2018-2020.QS, we have 23% and 46% of the nodes in the separator set, respectively.Now, we present the numerical results obtained by solving the m-VSP to the LVC for each QSRanking and MHES network; Table8 presents the cardinality value for sets A, B, and C and the percentage of nodes of LVC that belong to C. Here, it is important to mention that the size of the LVC for each network is provided in Table 8.Table 8
Size of LVC for each MHES and QSRanking network.
NetworkSize of LVCNetworkSize of LVC2013-2014.1532015-2016.2492014-2015.1472016-2017.2502015-2016.1492013-2015.2512016-2017.1512014-2016.2462013–2015.1532015-2017.2522014–2016.1472018.QS3192015–2017.1492019.QS3242013-2014.2512020.QS3282014-2015.247——Based on Table8, we can see that, for the MHES and QSRanking networks, the LVC is smaller compared to the MCGC. This is because the MCGC is composed of the union, and the LVC is composed at most by the union of all the giant components of every single layer.According to the results shown in Table9, we can assure that QSRanking and MHES networks are not robust; since for all QSRanking networks, C is between 28 and 33.8% and for all MHES networks, C is between 16.9% and 36.7% of the total nodes belonging to the LVC. For example, for the networks 2016.1-2017.1 and 2019-2020.QS, we have 31.3% and 35.2% of the nodes in the separator set, respectively.Table 9
Numerical results for m-VSP for QSRanking and MHES networks (LVC).
NetworkCA and B% of nodes of LVC inC2013.1-2014.1926 and 2516.92014.1-2015.11126 and 2323.42015.1-2016.11228 and 20242016.1-2017.11623 and 2131.32013.1–2015.11822 and 2033.92014.1–2016.11225 and 2325.52015.1–2017.11822 and 2036.72013.2-2014.21723 and 2033.32014.2-2015.21327 and 1927.62015.2-2016.21225 and 2324.42016.2-2017.21329 and 18262013.2–2015.21822 and 2035.32014.2–2016.21622 and 2234.82015.2–2017.21623 and 2130.72018.QS-2019.QS108142 and 12733.82019.QS-2020.QS114137 and 12635.22018.QS–2020.QS92149 and 13628Furthermore, we can verify that the methodology proposed in this work can find the set of nodes that cause the rupture of the MCGC and the LVC for different models of 2-layer and 3-layer multiplex networks since the sum ofA, B, and C is equivalent to the total number of nodes that belong to the MCGC or the LVC according to the case.It is important to mention that, for MHES networks, there is no difference in the networks modeled with both disaggregated and aggregated data, so they have the same numerical results when m-VSP is applied.The above means that both types of networks are equivalent and imply that it is possible to use a reduced set of information to analyze the behavior of the MHES; therefore, we can assure that the main characteristics of the complex system are properly structured when we modeled the QSRanking and MHES networks.Then, we can verify that the methodology proposed in this work can find the set of nodes that cause the rupture of the MCGC or the LVC for different models of 2-layer and 3-layer multiplex networks.In order to analyze the elements that cause the rupture of the MCGC and the LVC in the QSRanking and MHES networks, based on the results presented in [8, 28] about the classification of universities, we perform an analysis.In Table10, we show the percentage of the 377 universities that we consider to model the QSRanking networks that are located between positions 1 and 100; 101 and 200; 201 and 300; 301 and 400; and 401 and 500 in the main university rankings (like the QS ranking).Table 10
Percentage of universities in 100-universities ranges (QSRanking).
Network1–100 (%)101–200 (%)201–300 (%)301–400401–500 (%)2018-2019.QS20.729.328.415.26.42019-2020.QS18.429.328.717.5%6.12018–2020.QS19.527.53014%9Based on the information presented in Table10, we can see that most of the universities that cause the rupture of the network are between positions 100 and 300 of the QS ranking, and their main characteristics are as follows: average reputation, high employer reputation, average number of citations, few international faculties, and a greater number of national students.Now, in Table11, we show the percentage of universities with high, average, and low values in the university rankings that cause the disconnection of the MHES.Table 11
Percentage of universities with high, average, and low-ranking value (MHES).
NetworkHigh ranking (%)Average ranking (%)Low ranking (%)2013.1-2014.135.728.635.72014.1-2015.130.423.346.32015.1-2016.14033.326.62016.1-2017.136.626.636.62013.1–2015.133.32937.52014.1–2016.129.640.729.62015.1–2017.1332541.62013.1-2014.233.322.244.42014.1-2015.235.728.635.72015.1-2016.24033.326.72016.1-2017.242.126.331.62013.1–2015.234.630.734.62014.1–2016.233.329.134.62015.1–2017.238.430.830.8Based on the results of Table11, we can see that, for all MHES networks, the percentage of universities with high-, average-, and low-ranking values is divided equally, that is, a similar amount of the three types of universities are necessary to cause the rupture of the LVC.From a social network perspective, we can deduce that the links between the different types of universities are essential for the correct performance and development of the QSRanking and MHES networks.Finally, it is important to mention that the proposed methodology obtains a set of nodes that cause the rupture of the MCGC and the LVC without using structural network information or centrality measures.


## 5. Limitations of the Study and Discussion

The analysis developed in this work is based on the analysis of the rupture of the MCGC or the LVC belonging to each multiplex network. The idea of using the two ways is because the MCGC can give us an idea of the aggregated network; on the other hand, the LVC can give us the information about the nodes that are very important in the behavior of each layer.The idea of modeling two or three different types of layers for the MHES and QSRanking networks is to be able to obtain an identification of the best universities, based on a multicriteria analysis that reduces the disadvantages of using only one set of nodes; thus, the technique is capable of identifying those universities that, based on different characteristics, are the most important for the operation of the Mexican and international higher education system.We apply the multiplex network approach using the adaptation of the vertex separator problem to identify those nodes based on their main characteristics, when their links are eliminated from all the layers of the network, causing the rupture of the system.Therefore, a limitation of the study is that the drastic change of one or more variables by year changes the dynamics of the system.


## 6. Conclusions

This paper presents the adaptation of the vertex separator problem (VSP) for multiplex networks, denoted as m-VSP, which using inverse percolation obtains the nodes that, when removed from the system, cause the rupture of the mutually connected giant component (MCGC) or the large viable cluster (LVC).Also, based on the tests performed in this work, we can deduce that the main advantage of using m-VSP is that no structural calculations or centrality metrics are needed.Numerical results show that the m-VSP can work on different types of multiplex networks and that m-VSP can get a set of nodes that cause the rupture of the MCGC or the LVC.For QSRanking networks, the results show that most universities that cause the rupture of the LVC are between positions 101 and 300 in the ranking, and few universities are between positions 0-100 or 401-500.On the other hand, in the specific case of MHES networks, the results show that most universities that cause the rupture of the LVC are divided into three types: universities that obtain high, average, or low values in the university rankings.As future work, we propose the implementation of intelligence-based techniques of self-adaptive swarms and/or exact methods for solving the mathematical model m-VSP.



---

*Source: 8796360-2020-10-31.xml*