# Multilayer Network Science: from Cells to Societies

Oriol Artime1, Barbara Benigni1, Giulia Bertagnolli1, Valeria d’Andrea1,
Riccardo Gallotti1, Arsham Ghavasieh1, Sebastian Raimondo1, Manlio De
Domenico1,2 1Complex Multilayer Networks Lab, Fondazione Bruno Kessler, Via
Sommarive 18, 38123 Povo (TN), Italy  2Physics and Astronomy Department,
University of Padova, via Marzolo 8, 35131, Padova, Italy

(23 August 2022)

###### Abstract

Networks are convenient mathematical models to represent the structure of
complex systems, from cells to societies. In the past decade, multilayer
network science – the branch of the field dealing with units interacting in
multiple distinct ways, simultaneously – was demonstrated to be an effective
modeling and analytical framework for a wide spectrum of empirical systems,
from biopolymer networks (such as interactome and metabolomes) to neuronal
networks (such as connectomes), from social networks to urban and
transportation networks. In this Element, a decade after the publication of
one of the most seminal papers on this topic, we review the most salient
features of multilayer network science, covering both theoretical aspects and
direct applications to real-world coupled/interdependent systems, from the
point of view of multilayer structure, dynamics, and function. We discuss
potential frontiers for this topic and the corresponding challenges in the
field for the future.

##  I Introduction

What is a complex system? It is a _network_ of actors or units related by
special types of relationships or interactions which, together, form a whole.
Being two proteins within a cell or two individuals within a social group,
relationships and interactions tie together the units in such a way that “the
whole is larger than the sum of its parts”, a concept initially introduced by
the Greek philosopher Aristotle and later exploited by Gestalt psychologists,
at the end of 19th century, to explain human perception beyond the traditional
atomistic view.

In fact, the “whole” exhibits features that each actor or unit, in isolation,
does not (and also, could not) exhibit: therefore, it is usually difficult, if
not impossible, to understand a system from the analysis of its components
alone, like in atomistic or other reductionist theories. The framework
required to study such relationships and interactions is known as Network
Science111We refer the reader to this interesting, non-technical, and recent
introduction to the basic concepts characterizing complex systems
[complexity_explained_2019]..

The foundations of Network Science can be found in the pioneering work by
Leonhard Euler in 1736, when the famous mathematician provided the first
mathematically grounded proof to solve, definitively, the problem of the Seven
Bridges of Königsberg. He mapped the empirical problem of traversing the city
of Königsberg – under the constraint that one should use each one of its seven
bridges only one time – to the abstract problem of performing a special walk
through a graph. Since Euler’s solution, graph theory quickly developed in the
successive two centuries, culminating in the pioneering works by Paul Erdős
and Alfréd Rényi on random graphs and their statistical analysis at the end of
the 50’s of the past century.

Figure 1:  Multilayer representation of a co-authorship network. Nodes
represent authors publishing papers in the journals of the American Physical
Society, links connect two authors if they have published a paper together.
Layers encode distinct sub-topics of physics (e.g., Geophysics or Nuclear
Physics): links within the same layer represent co-authorship of one paper
about the same topic, while links between layers indicate co-authorship of one
paper categorized simultaneously across distinct topics. Figure from
[de2020illustrations] under Creative Commons Attribution-ShareAlike 4.0
International License.

For decades, graph theory has been widely used by social scientists and
(systems) biologists to map connections between individuals and biological
units, respectively, to gain novel insights about the _properties of a system_
, the relevance of a _unit within the system_ and the _organization of units
within the system_. In 1974, François Jacob, Nobel Prize in Physiology or
Medicine in 1965, described biology as the science dealing, effectively, with
systems within systems [trewavas2006brief], well before the age of genomics
and large-scale biology. He recognized that biological systems can be mapped
as well to units of systems at a larger scale: in fact, protein interact with
each other to make the cell function, cells interact with each other to make
tissues and organs, which in turn interact with each other to build an
organism. Finally, at the top of this hierarchical web of interactions,
organisms interact with each other to define a population, like our society.
In the same decade, similar ideas regarding the non-trivial interdependencies
between scales were laid out by the 1977 Nobel Laureate in Physics Phillip
Anderson, in the context of natural sciences [anderson1972more].

Social scientists, as biologists, were among the first to face the existence
of multiple levels (or scales) as well as multiple _layers_ of descriptions
for the units of a social system. In the early 70’s of the 20th century, Wayne
W. Zachary observed the interactions within a group of individuals belonging
to a Karate club, over 3 years [zachary1977information], to understand the
dynamics of conflicts which allowed him to predict the outcome of the group
split, happened later. He annotated interactions across eight distinct
contexts, from “the association in and between academic classes at the
university” to “attendance at intercollegiate karate tournaments held at local
universities”. However, at that time, the mathematical framework required to
study a network with multiple _layers of complexity_ – such as the eight
contexts – was not yet developed, and Zachary opted for an approximation:
aggregating the multiple interactions across contexts into a single
representative number denoting the intensity of the relationship between a
pair of actors.

Figure 2:  A multilayer transportation network, where connections using a
particular means of transport are associated with intra-layer links and inter-
modal exchanges are represented the inter-layer links. Here, the national
public transportation network for Great Britain [gallotti2014anatomy,
gallotti2015multilayer] as rendered by MuxViz [de2015muxviz]. (Figure from
[gallotti2015multilayer]).

In this work, we will better understand the challenges faced by systems
biologists and social scientists between the 70’s and the past decade, while
introducing the basic concepts required to define the framework of _multilayer
network science_ with an interdisciplinary language which should be familiar
to biologists, social scientists, computer scientists, applied mathematicians
and physicists. Therefore, it will become clear that, for instance, Zachary’s
approach was a possible model to study the Karate club network, but likely
neither the most accurate nor the most predictive one. We will discuss under
which conditions a system admits a multilayer representation, providing
examples such as the ones shown in Figs. 1 and 2, where units are individuals
and geographic areas, respectively, and interactions represent co-authorship
of scientific papers and transportation routes, respectively. Another
emblematic example, accounting for the temporal and socio-spatial
interdependence typical of many systems, concerns the organization of
ecological systems [pilosof2017multilayer]. Finally, very recently, multilayer
modeling in systems biology and medicine has been used to integrate
information about biological processes, drug-target, genotype and phenotype to
the subset of the human interactome targeted by SARS-CoV-2, the virus of the
COVID-19 [Verstraete2020] (see Fig. 3).

This work is full of examples like the ones above and we hope to make clear
the broad spectrum of potential interdisciplinary applications of the
multilayer framework.

Figure 3:  Illustration of CovMulNet19, the multilayer network encoding
COVID–19 Genotype-Phenotype-Drug interactions. A schematic map of intra- and
inter-layer interdependencies between diseases, symptoms, drugs, Gene Ontology
terms, human proteins and viral proteins of SARS-CoV-2, the virus of COVID–19.
Figure from [Verstraete2020].

Our ultimate goal is to guide the reader through the potential applications of
multilayer modeling, which nowadays provides a well-established paradigm for
the analysis of systems characterized by multiple levels and multiple layers
of description, including systems whose structure changes over time, with the
aim of providing the reader with the tools required to model and analyze their
systems in terms of coupled layers, as well as with the conditions under which
this approach is plausible or not.

It is worth remarking here that this work should be considered as an extended
introduction to the field but, at the same time, not the most complete one.
For this reason, we point the reader to the first reviews
[kivela2014multilayer, boccaletti2014structure, wang2015evolutionary,
de2016physics, battiston2017new] and recent books [bianconi2018multilayer,
cozzo2018multiplex] on this topic or more specifically on analysis and
visualization of multilayer networks [de2021multilayer], which, taken together
with our work, will provide a more comprehensive view of the field.

In Chap. 2 we will introduce the representation of multilayer networks based
on the tensorial formulation [de2013mathematical], providing the mathematical
ground for the analytical techniques for structure (Chap. 3) and dynamics
(Chap. 4), allowing the reader to find a reference for the analysis of
versatility (or multilayer centrality) and mesoscale organization (or
community detection), as well as for percolation, synchronization, competition
and modeling of intertwined phenomena. Towards the end (Chap. 5), we will
discuss a few selected advances in network science – namely the latent
geometry of a complex network based on network-driven processes and the
statistical theory of information dynamics leading to the formalism of network
density matrices – and will discuss their recent generalization and
application to multilayer networks. Finally (Chap. 6), we will show how
multilayer networks are ubiquitous and can be used for modeling complex
systems, from cells to societies.

##  II Representation of multilayer systems

###  II.1 Tensorial representation of a complex network

One convenient way to represent, mathematically, a complex network is by means
of its adjacency matrix [barrat2008dynamical, newman2018networks,
estrada2012structure, barabasi2016network, latora2017complex]. However, to
deal with multilayer networks, it might be more convenient to introduce first
the more general concept of tensor, a multilinear function which maps objects
defined in a vector space into other objects of the same type, regardless of
the choice of a coordinate system. For instance, a simple scalar $x$ is also a
rank-0 tensor, a vector $x_{i}$ is a rank-1 tensor and a matrix $X_{ij}$ is a
rank-2 tensor. More generally, given a vector space $\mathcal{V}$ with
algebraic dual space222This is the space of all the possible linear
transformations that map an object of $\mathcal{V}$ into a real number. For
instance, think about $\mathcal{V}=\mathbb{R}^{2}$ and the linear functional
$f:\mathbb{R}^{2}\longrightarrow\mathbb{R}$: it follows that $f(x,y)=ax+by$,
with $a,b$ two integer numbers, is an element of $\mathcal{V}^{\star}$.
$\mathcal{V}^{\star}$ over the real numbers $\mathbb{R}$, one can define the
tensor $M$ as the multilinear function

$\displaystyle
M:\mathcal{V}^{\star}\times\mathcal{V}^{\star}\times...\mathcal{V}^{\star}\times\mathcal{V}\times\mathcal{V}\times...\mathcal{V}\longrightarrow\mathbb{R},$
(1)

where the number of products is $m$ for the vector space and $n$ for its dual.
This definition formally characterizes a rank-$mn$ tensor
$M^{i_{1}i_{2}...i_{n}}_{j_{1}j_{2}...j_{m}}$ which is $m$-covariant and
$n$-contravariant: in fact, under a change of basis $B$, $m$ components
transform as the same linear mapping of the change of basis ($B$), whereas $n$
components transform as the inverse one ($B^{-1}$). Therefore, in general,
there are two types of canonical basis: the covariant basis denoted by
$e_{i}(a)$ ($a=1,2,...,m$) which is defined in $\mathcal{V}$, and the
contravariant (or dual) basis denoted by $e^{i}(b)$ ($b=1,2,...,n$) which is
defined in $\mathcal{V}^{\star}$. If the vector space is Euclidean, the
coordinates of the canonical vectors and their duals are the same, whereas
this is not the case in general. In the following, to define an adjacency
matrix, or a rank-2 adjacency tensor, we will work in the Euclidean space but
we will keep the covariant and contravariant notation, since it will allow us
to generalize the results to the case of non-Euclidean spaces. The interested
reader can find more about the tensorial framework in any good linear algebra
textbook, while for the purpose of this work it is sufficient to understand
how we can use tensors in practice in a few key situations.

Let us start by better defining the canonical vectors in the case of networks.
For a graph with $N$ nodes, the canonical covariant vectors $e_{i}(a)$ defined
in the space of nodes $\mathbb{R}^{N}$ are $N$ rank-1 tensors of dimension $N$
with all entries equal to 0 except for the $a$-th entry, which is equal to 1.
Similarly for canonical contravariant vectors. The product of canonical
vectors gives canonical matrices, e.g., $E_{ij}(ab)=e_{i}(a)e_{j}(b)$ is a
rank-2 covariant tensor with all components equal to zero except for the one
corresponding to the $a$-th row and the $b$-th column, equal to 1. Similarly,
we can build contravariant tensors and mixed tensors, i.e., tensors obtained
by the product between the covariant and contravariant vectors.

The careful reader has noticed, at this point, that we have defined the
_outer_ product of two canonical vectors, also known as the Kronecker product,
which gives a rank-2 tensor as a result. This result is general: the outer
product of two tensors $X$ and $Y$ is a new tensor $Z$ with a number of
covariant (contravariant) indices given by the sum of the number of covariant
(contravariant) indices of $X$ and $Y$. Therefore, the outer product of two
tensors is always a tensor of higher order than the original ones: e.g.,
$X_{ij}^{k}Y_{l}^{mn}=Z_{ijl}^{kmn}$.

It is possible to define also an _inner_ product: in this case, we talk about
a contraction because the rank of resulting tensor is reduced by two units.
For instance, this is the case in the product
$X_{ij}^{k}Y_{k}^{mn}=Z_{ij}^{mn}$, where the index $k$ is covariant for $X$
and contravariant for $Y$. This operation corresponds to summing over the
components of $X$ and $Y$ identified by the index $k$. The careful reader has
noticed that we have omitted the summation symbol: this choice – known as
Einstein summation convention – is optional and often adopted for sake of
simplicity. In the following, we will make use of this convention.

At this point, we are ready to define the adjacency tensor of a complex
network in terms of canonical vectors [de2013mathematical] as

$\displaystyle
W^{i}_{j}=\sum_{a,b=1}^{N}w_{ab}e^{i}(a)e_{j}(b)=\sum_{a,b=1}^{N}w_{ab}E^{i}_{j}(ab),$
(2)

where $w_{ab}$ is a real number, usually non-negative, used to encode the
intensity of the interaction between nodes $a$ and $b$, while
$E^{i}_{j}(ab)\in\mathbb{R}^{N\times N}$ are the mixed canonical rank-2
tensors. We might wonder if $W^{i}_{j}$ is a true tensor, or just a matrix. To
this end, it is enough to understand how it transforms under a change of basis

$\displaystyle B^{i}_{j}=\sum_{a=1}^{N}e^{\prime i}(a)e_{j}(a),$ (3)

a linear function which transforms the basis vector set $\\{e^{i}(a)\\}$ into
a second set $\\{e^{\prime i}(a)\\}$. By noting that $w_{ab}$ must be
invariant with respect to the change of basis, we have:

$\displaystyle{W^{\prime}}^{k}_{l}$ $\displaystyle=$
$\displaystyle\sum_{a,b=1}^{N}w_{ab}e^{\prime
k}(a)e^{\prime}_{l}(b)=\sum_{a,b=1}^{N}w_{ab}B^{k}_{i}e^{i}(a)e_{j}(b)(B^{-1})^{j}_{l}$
(4) $\displaystyle=$ $\displaystyle
B^{k}_{i}\left[\sum_{a,b=1}^{N}w_{ab}e^{i}(a)e_{j}(b)\right](B^{-1})^{j}_{l}=B^{k}_{i}W^{i}_{j}(B^{-1})^{j}_{l},$

i.e., the adjacency object $W^{i}_{j}$ transforms like a tensor
[de2015ranking]. This result is important since a tensor is an object with
features that, in general, are not shared by a matrix or, at higher orders, a
hypermatrix. In fact, the components of a tensor can be always arranged into
hypermatrices, while the opposite is not necessarily true.

Since we work in the Euclidean space, we might wonder why we use this notation
and not a simpler one. In general, this is convenient because of the presence
of directed relationships between nodes: to distinguish between incoming and
outgoing directions, it is sufficient to map this information into covariant
and contravariant indices in such a way that the adjacency tensor $W^{i}_{j}$
represents a linear transformation which maps nodes into a function of their
incoming or outgoing flow. For instance, node $a$ is represented by $e_{i}(a)$
in the space of nodes and $W^{i}_{j}e_{i}(a)=w_{j}(a)$ provides a rank-1
tensor encoding the set of nodes linked by $a$, while $W^{i}_{j}u_{i}=s_{j}$,
with $u_{i}$ the rank-1 tensor with all components equal to 1, provides a
rank-1 tensor encoding the outgoing strength of all nodes. Similarly,
$W^{i}_{j}e^{j}(a)=w^{i}(a)$ gives the set of nodes linking to $a$, while
$W^{i}_{j}u^{j}=s^{i}$ gives the incoming strength of nodes.

Before moving to the next section, it is useful to define some tensors used
throughout this work. We have just seen the rank-1 1–tensor in action:
similarly we can define the rank-2 1–tensor $U^{i}_{j}=u^{i}u_{j}$ or higher-
order tensors. Another fundamental tensor is the Kronecker one, defined by
$\delta^{i}_{j}$, with components equal to 1 if $i=j$ and equal to 0
otherwise.

Figure 4:  A system where nodes are characterized by three distinct types of
interactions, encoded by colored layers. Overall, the system is a multilayer
network, because to describe relationships we need to specify more than one
network. Units are physical nodes: each one is a set of _state nodes_ or
_replicas_ , each one encoding the identity of the corresponding physical node
in each layer separately. _Intra-layer_ edges define connectivity within each
layer, whereas _inter-layer_ edges define connectivity across layers.
Reproduced with permission from [de2020illustrations].

###  II.2 Tensorial representation of a multilayer network

In the previous section we have introduced the fundamental procedure required
to build an adjacency tensor to represent a classical network (a monoplex).
Using a similar procedure, we can build a _multilayer adjacency tensor_ to
represent a multilayer network, as the one shown in Fig. 4. A multilayer
system is characterized by $N$ physical nodes interacting in $L$ distinct ways
simultaneously: each type of interaction defines a _layer_. At variance with
single-layer networks, there are more edges sets to encode: as many as the
number ($L$) of layers and, in general, as many as the number ($L(L-1)$) of
directed pairwise connections between layers, since we have to specify which
node $i$ in a layer $\alpha$ is connected to which node $j$ in a layer $\beta$
($i,j=1,2,...,N$, $\alpha,\beta=1,2,...,L$)333This simple observation suggests
that a good candidate for multilayer adjacency tensor should be a rank-4
tensor.. Note that for simplicity, we are indicating with Greek letters the
indices related to layers and with Latin letters the indices related to nodes.

There are different types of multilayer networks, depending on the presence or
absence of links between layers and on the way nodes are defined (see Fig. 5).
In the following, we will mostly deal with the class of systems characterized
by inter-layer connectivity, since it is not possible to define a meaningful
multilayer adjacency tensor for the class of edge-colored multigraphs444Note
that, instead, it is possible to define a valid hypermatrix encoding this
object, and this hypermatrix can be thought of as an array of matrices
[bianconi2013statistical]..

Figure 5:  Multilayer networks include a broad spectrum of possible models.
Edge-colored networks are useful models when inter-layer connectivity is not
well defined: this is the case of a social network where edges can represent
different types of social relationships (e.g., trust, family, business, etc.)
[cardillo2013emergence, nicosia2015measuring, de2015structural]. Conversely,
in interconnected networks inter-layer connectivity is well defined and allows
us to model a variety of systems [de2014navigability, de2015ranking,
osat2017optimal, radicchi2017redundant], including those with
interdependencies where nodes control and/or are controlled by nodes in
another network [rosato2008modelling, buldyrev2010catastrophic,
vespignani2010complex, gao2012networks, radicchi2015percolation]. Reproduced
with permission from [de2020illustrations].

Let us introduce the canonical rank-1 vectors $e^{\alpha}(p)$
($\alpha,p=1,\ldots,L$) in the space of layers $\mathbb{R}^{L}$, and the
corresponding canonical rank-2 tensors
$E^{\alpha}_{\beta}(pq)=e^{\alpha}(p)e_{\beta}(q)$, similarly to what we have
done for monoplexes. It is straightforward to show [de2013mathematical] that
the linear combination of

$\displaystyle
M^{i\alpha}_{j\beta}=\sum_{a,b=1}^{N}\sum_{p,q=1}^{L}w_{ab}(pq)e^{i}(a)e_{j}(b)e^{\alpha}(p)e_{\beta}(q)$
(5)

fully characterizes a multilinear object in the space $\mathbb{R}^{N\times
L\times N\times L}$. This object is, in fact, the desired multilayer adjacency
tensor since, under a change of coordinates, it transforms like a tensor:

$\displaystyle M^{\prime i\alpha}_{j\beta}$ $\displaystyle=$
$\displaystyle\sum_{a,b=1}^{N}\sum_{p,q=1}^{L}w_{ab}(pq)B^{i}_{k}e^{k}(a)(B^{-1})^{l}_{j}e_{l}(b){\tilde{B}}^{\alpha}_{\gamma}e^{\gamma}(p)({\tilde{B}}^{-1})^{\delta}_{\beta}e_{\delta}(q)$
(6) $\displaystyle=$ $\displaystyle
B^{i}_{k}{\tilde{B}}^{\alpha}_{\gamma}M^{k\gamma}_{l\delta}(B^{-1})^{l}_{j}({\tilde{B}}^{-1})^{\delta}_{\beta}.$

By indicating with
$E^{i\alpha}_{j\beta}(ab;pq)=E^{i}_{j}(ab)E^{\alpha}_{\beta}(pq)$ the
canonical rank-4 tensors, we can simply reduce the definition of the
multilayer adjacency tensor to

$\displaystyle
M^{i\alpha}_{j\beta}=\sum_{a,b=1}^{N}\sum_{p,q=1}^{L}w_{ab}(pq)E^{i\alpha}_{j\beta}(ab;pq),$
(7)

where $w_{ab}(pq)$ encodes the intensity of the interaction between node $a$
in layer $p$ and node $b$ in layer $q$. Note that $w_{ab}(pp)$ indicates the
weights of the links in layer $p$.

It is worth noticing that, as for the space of nodes, in the space of layers,
we can define multilayer 1–tensors and Kronecker tensors as
$U^{i\alpha}_{j\beta}=U^{i}_{j}U^{\alpha}_{\beta}$ and
$\delta^{i\alpha}_{j\beta}$, respectively. Another important tensor,
representing a complete multilayer network without self-edges, will be used
later in this work to characterize multilayer triadic closure: for
consistency, we prefer to introduce it here as
$F^{i\alpha}_{j\beta}=U^{i\alpha}_{j\beta}-\delta^{i\alpha}_{j\beta}$.

At this point, the reader should be familiar enough with tensors to note that
different decompositions are possible. Here, we are not referring to
operations like Tucker decomposition – the higher-order generalization of
singular value decomposition (SVD) [tucker1966some] – but to a linear
decomposition to highlight the fundamental components of a multilayer system.
In fact, we can identify four tensors which encode distinct structural
information:

$\displaystyle m^{j\beta}_{i\alpha}$
$\displaystyle=\underbrace{m^{j\beta}_{i\alpha}\delta_{\alpha}^{\beta}\delta_{i}^{j}+m^{j\beta}_{i\alpha}\delta_{\alpha}^{\beta}(1-\delta_{i}^{j})}_{\text{intra-
layer
relationships}}+\underbrace{m^{j\beta}_{i\alpha}(1-\delta_{\alpha}^{\beta})\delta_{i}^{j}+m^{j\beta}_{i\alpha}(1-\delta_{\alpha}^{\beta})(1-\delta_{i}^{j})}_{\text{inter-
layer relationships}}$ (8)
$\displaystyle=\underbrace{m^{i\alpha}_{i\alpha}}_{\text{self-
relationships}}+\underbrace{m_{i\alpha}^{j\alpha}}_{\text{endogenous}}+\underbrace{m_{i\alpha}^{j\beta}}_{\text{exogenous}}+\underbrace{m_{i\alpha}^{i\beta}}_{\text{intertwining}}$
$\displaystyle=\mathbb{S}_{i\alpha}(M)+\mathbb{N}^{j}_{i\alpha}(M)+\mathbb{X}^{j\beta}_{i\alpha}(M)+\mathbb{I}_{i\alpha}^{\beta}(M)\,.$

Here, the components of the tensor are indicated by $m_{i\alpha}^{j\beta}$
($i,j=1,2,\dots,N$ and $\alpha,\beta=1,2,\dots,L$), while $\delta^{j}_{i}$ and
$\delta_{\alpha}^{\beta}$ indicate the Kronecker delta function in the space
of nodes and layers, respectively. The four tensors encode the following
relationships:

  * •

Intra-layer interactions:

    * –

self-interactions ($\mathbb{S}$): from a node to itself;

    * –

endogeneous interactions ($\mathbb{N}$): between distinct nodes belonging to
the same layer;

  * •

Inter-layer interactions:

    * –

exogenous interactions ($\mathbb{X}$): between distinct nodes belonging to
distinct layers;

    * –

intertwining ($\mathbb{I}$): from a node to its replicas in other layers.

Equation (8) characterizes the “structural
$\mathbb{S}\mathbb{N}\mathbb{X}\mathbb{I}$ decomposition” of the multilayer
adjacency tensor $M$: the models for interconnected systems described in Fig.
5 can be characterized by which $\mathbb{S}\mathbb{N}\mathbb{X}\mathbb{I}$
components contribute to the tensor. In particular, we identify:

  * •

Interconnected multiplex networks: type $\mathbb{S}\mathbb{N}\mathbb{I}$;

  * •

Interdependent networks: type $\mathbb{S}\mathbb{N}\mathbb{X}$;

  * •

General multilayer networks: $\mathbb{S}\mathbb{N}\mathbb{X}\mathbb{I}$.

As previously mentioned, edge-colored networks do not admit a meaningful
representation in terms of multilayer adjacency tensor, but according to this
classification they would define type $\mathbb{S}\mathbb{N}$.

Figure 6:  A rank-2 tensor defined in $\mathbb{R}^{4\times 4}$ is flattened
into a rank-1 tensor defined in $\mathbb{R}^{1\times 16}$, without loss of
information. Reproduced with permission from [de2020illustrations].

From an operational perspective, working with tensors might be complicated and
cumbersome. Nevertheless, from a theoretical perspective, the tensorial
formulation allows us to write complex equations in a very compact and handy
way, and it can be used to guide our intuition about generalizing existing
network descriptors for multilayer analysis, as we will see in the next
chapters. One widely adopted approach is based on the operation known as
_matricization_ (or _flattening_) [kolda2009tensor], which maps a complex
object like a high-order tensor into a lower-order object while preserving the
information content (see Fig. 6 for an emblematic example). In the case of the
multilayer adjacency tensor $M^{i\alpha}_{j\beta}$, which is defined in the
space $\mathbb{R}^{N\times L\times N\times L}$, this operation corresponds to
flatten entries into a rank-2 object defined in the space
$\mathbb{R}^{NL\times NL}$, as shown in Fig. 7. It is manifest that there is
no loss of information, although care is needed when dealing with this new
object, which is known in the literature as _supra-adjacency matrix_
[gomez2013diffusion, de2013mathematical, sole2013spectral].

Figure 7:  As in Fig. 6 but applied to a multilayer adjacency tensor
$M^{i\alpha}_{j\beta}$ in $\mathbb{R}^{2\times 4\times 2\times 4}$ flattened
into a supra-adjacency matrix in $\mathbb{R}^{8\times 8}$. Colors highlight
the diagonal blocks where adjacency matrices corresponding to layers are
encoded. Off-diagonal blocks encode inter-layer connectivity. Reproduced with
permission from [de2020illustrations].

Working with a supra-adjacency matrix comes with several computational
advantages and notational disadvantages. For practical purposes, the
flattening allows for a visual inspection of the multilayer adjacency tensor,
as shown in Fig. 8. The supra-adjacency matrix is, in fact, a block matrix
with intra-layer connectivity encoded into diagonal blocks and inter-layer
connectivity encoded into off-diagonal ones (see also Fig. 8). However, it is
worth remarking that this is a matter of convention: in fact, this arrangement
of blocks is not unique and other arrangements are also valid, although some
ones are less convenient than others for the design of algorithms. This non-
uniqueness of the supra-adjacency matrix makes it less suitable for
theoretical calculations but still an alternative to higher-order tensors.
Fig. 9 shows the supra-adjacency matrix representation for three widely used
multilayer network models.

Figure 8:  A multilayer network, consisting of 3 layers, flattened into a
supra-adjacency matrix [gomez2013diffusion, de2013mathematical]. (A) The
interconnected multiplex network (left) with $N=4$ nodes and $L=3$ layers, the
latter ones are color-coded. Connectivity is weighted (see edge thickness) and
directed (see arrows): the adjacency matrix of each layer is shown on the
right-hand side of this panel. (B) Result of the flattening procedure: the
original rank–4 tensor is now encoded into a rank–2 supra-adjacency matrix
with a block structure. Intra-layer connectivity is encoded into diagonal
blocks (according to colors), whereas inter-layer interactions are represented
into off-diagonal blocks. Reproduced with permission from
[baggio2016multiplex].

In the next chapters we will use the tensorial formulation, when possible, for
the analysis of the topology of multilayer networks and to define dynamical
processes.

Figure 9: Supra-adjacency matrix representation of three distinct types of
multilayer networks [gomez2013diffusion, de2013mathematical], where layers are
encoded by colors: (a) edge-colored multigraph (no inter-layer connectivity);
(b) interconnected multiplex network (diagonal coupling in the off-diagonal
blocks), and (c) general interconnected case, where inter-layer connections
are not restricted to replicas (exogenous interactions). Latin letters denote
nodes, while Greek letters are used for layers. $D(i;\alpha\beta)$ indicates
the intensity of the connection between state nodes $(i,\alpha)$ and
$(i,\beta)$. Figure from [bertagnolli2020diffusion].

##  III Multilayer structural analysis

In this chapter, we will operationally define some of the most important
theoretical and computational tools for the analysis of a multilayer network.
We start with walks (Sec. III.1) to define distinct types of connected
components (Sec. III.2) and we will quickly move to describe several measures
of node importance within the system, i.e., multilayer versatility, which is
the generalization of the concept of node centrality (Sec. III.3).

It will follow an overview of methods to characterize the mesoscale
organization of a system (Sec. III.4), by identifying how nodes group together
into small clusters (Sec. III.4.1), such as triangles, and in larger
functional modules or communities (Sec. III.4.2). We will show how to use
these concepts to understand how units in multilayer systems are integrated or
segregated, in terms of information flow (Sec. III.4.3). Technically speaking,
connected components should be described here, but we opted to keep the
corresponding section after the description of walks for sake of simplicity.

We will conclude this chapter by discussing existing methods to quantify the
correlations between pairs of layers (Sec. III.5): this class of analytical
tools is fundamental to better understand results from versatility and
mesoscale analysis, since the existence of layer-layer correlation patterns is
reflected in structural measures based on walks and paths.

###  III.1 Basic definitions

In a generic network, a _walk_ is defined as a sequence of adjacent nodes and
edges visited by a hypothetical walker. In general, a given edge (but also
node) can be traversed more than once, but it is possible to apply some
restrictions to a multilayer walk to specify special walks. For instance, one
can define a _multilayer trail_ as a walk where links can be traversed only
one time. A further restriction can be also applied to the identity of origin
and destination nodes, and we can define a _multilayer path_ as a multilayer
trail with the restriction that repeated nodes are not allowed, whereas a
_multilayer cycle_ is a closed trail where only the origin and destination
nodes are repeated. Finally, a closed trail with more than one repeated node
is defined as a _multilayer circuit_. An illustration of different types of
walks on a multilayer network is shown in Fig. 10.

The _length of a walk_ is the number of edges traversed along the walks. It is
possible to calculate the number of walks of length $\ell$ from a node $i$ to
any other node $j$. For an unweighted network the element $A^{k}_{l}$ of its
adjacency matrix is 1 if there is an edge from $k$ to $l$, and is 0 otherwise.
Then, if there is a walk with $\ell=2$ from $i$ to $j$ via $k$, the product
$A^{i}_{k}A^{k}_{j}$ will be 1. By generalizing this argument, we can write
the rank–2 tensor encoding information about walk length between any pair of
nodes in the network as the $\ell$–th power of the rank–2 adjacency tensor
representing the network:

$\displaystyle\mathcal{W}^{i}_{j}(\ell)=(A^{i}_{j})^{\ell}=A^{i}_{j_{1}}A^{j_{1}}_{j_{2}}\dots
A^{j_{\ell-1}}_{j}$ (9)

The same formalism can be used in a weighted network, by defining the weight
of a walk as the product of the weights of the traversed links. The entries of
$\mathcal{W}^{i}_{j}(\ell)$ will give the sum of weights of the walks of
length $\ell$ connecting the corresponding pair of nodes.

An analogous approach is used to calculate the walk length for multilayer
networks. If $M^{i\alpha}_{j\beta}$ is the rank–4 adjacency tensor
representing the system, then the entries of the $\ell$–th power of this
tensor provide the number of multilayer walks of length $\ell$ between a node
$i$ in layer $\alpha$ and a node $j$ in layer $\beta$:

$\mathcal{W}^{i\alpha}_{j\beta}(\ell)=M^{i\alpha}_{j_{1}\beta_{1}}M^{j_{1}\beta_{1}}_{j_{2}\beta_{2}}\dots
M^{j_{\ell-1}\beta_{\ell-1}}_{j\beta}.$ (10)

The aggregate representation of multilayer networks allows us to aggregate
layers to obtain a network where the number of edges or the weight of an edge
is the number of different types of edges between a pair of nodes
[de2013mathematical, cozzo2015structure]. It is possible to highlight the
topological difference between multilayer networks and their aggregated
representations by using the aforementioned formalism. Let
$\bar{G}^{i}_{j}=M^{i\alpha}_{j\beta}U^{\beta}_{\alpha}$ be the aggregate
network which accounts for inter-layer links: the corresponding rank–2 walk
tensor is then given by

$\displaystyle\bar{\mathcal{W}}^{i}_{j}(\ell)$ $\displaystyle=$
$\displaystyle(\bar{G}^{i}_{j})^{\ell}=\bar{G}^{i}_{j_{1}}\bar{G}^{j_{1}}_{j_{2}}\dots\bar{G}^{j_{\ell-1}}_{j}$
(11) $\displaystyle=$ $\displaystyle
M^{i\alpha}_{j_{1}\beta_{1}}U^{\beta_{1}}_{\alpha}M^{j_{1}\beta_{1}}_{j_{2}\beta_{2}}U^{\beta_{2}}_{\beta_{1}}\dots
M^{j_{\ell-1}\beta_{\ell-1}}_{j\beta}U^{\beta}_{\beta_{\ell-1}}$
$\displaystyle=$
$\displaystyle\underbrace{\left(M^{i\alpha}_{j_{1}\beta_{1}}M^{j_{1}\beta_{1}}_{j_{2}\beta_{2}}\dots
M^{j_{\ell-1}\beta_{\ell-1}}_{j\beta}\right)}_{\mathcal{W}^{i\alpha}_{j\beta}(\ell)}\underbrace{\left(U^{\beta_{1}}_{\alpha}U^{\beta_{2}}_{\beta_{1}}\dots
U^{\beta}_{\beta_{\ell-1}}\right)}_{U^{\alpha}_{\beta}L^{\ell-1}},$

showing that the number of walks with length $\ell$ between two nodes in the
aggregate network is not a linear function of the number of walks of length
$\ell$ between the same pair of nodes in the multilayer network.

Figure 10:  Illustration of distinct walks possible on a multilayer network.
By definition, a walk in a multilayer system is the most general way to move
between its layers, nodes and links. To define walks characterized by distinct
features, it is sufficient to restrict the number of times nodes or links can
be traversed or to put constraints on origin and destination nodes. The
illustration shows a sequence of nodes and edges corresponding to each type of
walk (multilayer walk, trail, path, cycle and circuit) with a dashed line:
note that the sequence of nodes visited by each type of walk is explicitly
reported below the corresponding multilayer network. Reproduced with
permission from [de2020illustrations].

In a network, a walk that does not intersect itself is named _path_ and the
_shortest path_ is the shortest walk between a given pair of nodes. Shortest
paths have an important role in several network phenomena: they allow us to
model, for instance, how information is exchanged between two nodes by using
the least number of traversed nodes and links. In undirected networks the
length of the shortest path is often used to define a distance between nodes
whereas some properties of a _geodesic distance_ are, in general, no more
satisfied in presence of directed links.

###  III.2 Connected components

To identify clusters of nodes555Note that this notion is different from the
one of groups or communities or modules, although the terminology might be
sometimes misleading. that can exchange information in a network, it is useful
to analyze the connected components [newman2018networks, barabasi2016network].
Components are defined as separate parts of the network, i.e., subset of nodes
with at least one path between any origin/destination pair belonging to the
subset. A network with a single component is _connected_ whereas, if there is
more than one component, a network is _disconnected_. For instance, isolated
nodes count as disconnected components of the system. For a directed network,
a component is defined as _weakly connected_ if two nodes of the undirected
representation of the component are connected by one or more paths. If there
is a directed path in both directions between every pair of nodes, the
component is defined as _strongly connected_.

For networks with finite size, we define the _largest connected component_
(LCC) as the cluster with the maximal subset of nodes. If the size of the
network is infinite, such as in the thermodynamic limit, the LCC is usually
named _giant connected component_ (GCC).

In interdependent networks [gao2012networks] two systems $A$ and $B$ are
interconnected with links and the potentially functional clusters are
identified by _mutually connected components_. If we indicate by $\mathcal{A}$
the set of nodes in network $G(A)$ and by $\mathcal{B}$ the corresponding set
of nodes in network $G(B)$, they form a mutually connected component if:

  * •

each pair of nodes in $\mathcal{A}$ is connected by a path consisting of nodes
belonging to $\mathcal{A}$ and links of network $G(A)$;

  * •

each pair of nodes in $\mathcal{B}$ is connected by a path consisting of nodes
belonging to $\mathcal{B}$ and links of network $G(B)$
[buldyrev2010catastrophic].

In multilayer networks we can use the definition of path described in Sec.
III.1 to define a _multilayer connected component_ as the subset of nodes
connected by a multilayer path [de2014navigability]. This definition can be
used to identify connected components from the aggregate representation of the
multilayer network, because it is sufficient that a pair of nodes is connected
by a path to be part of the same component.

Figure 11:  In many cases of interest, multilayer networks can be represented
in different ways, e.g., (top-left) in terms of interconnected layers – i.e.,
a multiplex – or (bottom-left) by collapsing state nodes and explicitly
encoding layers as colors of an edge-colored multigraph. To define components,
it is important to identify origin and destination nodes and to consider paths
connecting them. A multilayer connected component is identified by paths
between state nodes. On the right-hand side of the figure, three distinct
definitions of multilayer connected components are illustrated: largest
connected component (top), largest intersection component (middle) and largest
viable component (bottom). See the text for details. Reproduced with
permission from [de2020illustrations].

However, different and complementary information about the structure of a
multilayer network can be provided by more restrictive definitions. For
instance, the _largest intersection component_ (LIC) is defined as the largest
cluster in which nodes are connected across all layers simultaneously
[baxter2012avalanche] and can be identified by intersecting the LCC of each
layer separately. Another alternative is to aggregate the multilayer system
with respect to the intersection of edges and then identify the largest
connected component of the resulting network [radicchi2015percolation]. That
definition has been recently used to better understand the emergence of
continuous or abrupt percolation phase transitions even in systems of finite
size.

However, one more definition is possible and useful for applications. The
_largest viable component_ (LVC) is defined as the maximal subset of _viable_
nodes and it consists of nodes that are connected by a path in each layer
simultaneously [baxter2012avalanche]. As a result, all nodes in the LVC are
essential to the function of the system and to define its structural core
[azimi2014k, grassberger2015percolation, stella2018multiplex]. The more
restrictive condition imposed by the LVC is responsible for a hybrid phase
transition which leads to the discontinuous emergence of the giant viable
cluster, at variance with ordinary percolation where a continuous phase
transition is observed [baxter2012avalanche].

Fig. 11 shows an illustration of multilayer connected components for both
interconnected multiplex and edge-colored multigraph representations. The more
restrictive the condition imposed, the smaller the size of the largest
connected cluster: in particular, the size of the LVC is equal to or smaller
than the size of the LIC, which in turn is equal or smaller than the size of
the LCC.

###  III.3 Measuring influence: versatility

Information exchange with a multilayer network is strongly dependent on the
organization of the underlying system into connected components. However,
given a (possibly connected) network, it is plausible to wonder _which
node(s)_ are more important than others for information flow. For practical
applications, several additional questions are plausible, and each question
encodes an operational definition of _node centrality_.

Revealing the most central nodes in complex networks is a key issue in a
variety of real-world scenarios [barabasi2016network]. In epidemiology, for
example, finding the most central multilayer nodes helps identifying the
pivotal disease spreaders [de2016physics], while in cascading failures, by
detecting the most central nodes we can recognize the most fragile actors able
to trigger the failure of other parts of the network
[buldyrev2010catastrophic].

Over the years, a wide spectrum of network centrality measures has been
proposed for single-layer networks (for interested readers, we refer to a
review from the computational social sciences [scott2017popularity]). The
generalization of such descriptors to the realm of multilayer systems is not
always immediate, since nodes peripheral in one layer might be extremely
central in another one [de2013mathematical, kivela2014multilayer,
boccaletti2014structure], and crucial nodes for a given dynamics might not be
central across all layers. If, for example, two distinct layers have only one
node in common, any exchange of information between those layers will involve
the passage through that common node, independently from its centrality in
each layer separately: hence, such a node will be highly central for the
considered process. _Multilayer versatility_ [de2013mathematical,
de2015ranking] is a measure able to capture these features of a node,
quantifying how important nodes are for the processes which characterize the
definition of specific centrality descriptors – e.g., in terms of information
diffusion or flow through shortest paths – in a multilayer network.

In 2006, Borgatti and Everett [borgatti2006graph] approached centrality from a
graph-theoretic perspective, claiming that all centralities consider the
involvement of a node in the walk structure of the network according to four
features: Walk Type, Walk Property, Walk Position and Summary Type. The Walk
Type dimension concerns the kind of walks considered as well as the kind of
constraints on such walks, for example considering only shortest-paths. The
Walk Property dimension distinguishes between measures that consider the
volume of walks a node is involved in (e.g., as in betweenness centrality) and
measures that regard the length of those walks (e.g., as in closeness
centrality). The Walk Position can be radial, where walks start in a given
node (e.g., closeness) or medial, where walks pass through that given node
(e.g., betweenness). Finally, the Summary Type dimension distinguishes
measures according to the chosen summary statistic used to obtain a centrality
score vector. Based on the possible combinations of these four features –
i.e., Walk Type, Walk Property, Walk Position and Summary Type – different
centrality measures can be obtained.

Besides identifying the most central nodes, the power of versatility lies also
in predicting their role in some emblematic dynamical processes, such as
diffusion and congestion [de2015ranking]. The importance of this aspect is
highlighted, for example, when nodes of a multilayer network are ranked by
their coverage [de2014navigability] (see also Sec. IV.1) using the PageRank –
see later in this section for its definition – obtained from the multilayer
model and its aggregated network. PageRank versatility, obtained from the
multilayer network, is a better estimator of the evolution of such dynamics,
outperforming its single-layer counterpart.

In the following, we present the most important centrality measures in the
case of multilayer networks, i.e. Multidegree, K-coreness, Eigenvector and
Katz versatility, HITS versatility, Random walk occupation centrality,
PageRank versatility, Random walk betweenness and closeness versatility,
Betweenness versatility and Interdependence centrality, while we refer the
interested reader to other interesting measures which exploiting multiplex
features to identify a node’s importance [posfai2019consensus], recently
validated to predict multiplex centrality in the rhesus macaque
[beisner2020multiplex].

####  III.3.1 Multidegree

Multidegree centrality ($k_{i}$) is the simplest indicator of node importance
at the local level, it is obtained by summing up all the links connected to
node $i$ across all layers [de2013mathematical]:

$\displaystyle
k_{i}=\sum_{\alpha,\beta=1}^{L}\sum_{j=1}^{N}M_{j\beta}^{i\alpha}=M_{j\beta}^{i\alpha}u_{\alpha}^{\beta}u^{j}.$
(12)

In some cases, e.g., when interlayer links are not explicitly considered, it
can be more suitable to evaluate the degrees coming from individual layers:

$\displaystyle k_{i}^{\alpha}=\sum_{j=1}^{N}A_{j}^{i}(\alpha),$ (13)

where $A_{j}^{i}(\alpha)$ represents the adjacency tensor of layer $\alpha$.
Other definitions partially related to this one, for the case of non-
interconnected multiplex networks, can be found in [bianconi2013statistical,
battiston2014structural].

####  III.3.2 K-coreness

In single layer network, the $k$-core of a graph is defined as a maximal
connected subgraph in which each vertex has at least degree $k$ within that
subgraph. The ensemble of all $k$-cores of a graph represents the core
decomposition of that specific graph [seidman1983network]. To extend this
concept to multilayer networks, we have to take into account the different
types of edges, encoding different types of interactions encoded into layers.
In this case, the $k$-core is defined as the largest subgraph in which each
vertex has at least $k_{\alpha}$ edges of each type, $\alpha=1,2,...,M$ where
$M$ is the total number of layers [azimi2014k] (see Fig. 12). For more detail
on how to efficiently compute the complete core decomposition of a multilayer
network, we refer to [galimberti2020core].

Figure 12: The ($k_{1},k_{2}$)-core decomposition of a multilayer network
characterized by two distinct types of edges. The first type is encoded by
solid black links, while the second type is encoded by dashed blue links. The
concentric circles help identify the ($k_{1},k_{2}$)-core decomposition: the
cores from the most external circle to the innermost are the (1,1)-core, the
(1,2)-core, the (2,2)-core, and the (1,3)-core, respectively. Figure from
[azimi2014k].

####  III.3.3 Eigenvector Versatility

Eigenvector centrality is a measure of influence of the single nodes in a
network. A particular node has high eigenvector centrality if it is high the
centrality of its neighbours [bonacich1987power]. In the monoplex case, the
recursive character of this definition is untangled by the eigenvalue problem:

$\displaystyle W{\bf x}=\lambda_{1}{\bf x},$ (14)

where $\lambda_{1}$ is the largest eigenvalue of $W$ and the element $x_{i}$
represents the centrality score of the node $i$ in the network described by
the weight matrix W [de2013mathematical, de2015ranking]. If $W$ is symmetric
with positive entries, the Perron-Frobenius theorem grants the existence and
uniqueness of this vector.

When complete knowledge of the intra-layer connectivity is available, a
natural extension of the definition of eigenvalue centrality for multilayer
networks can be easily obtained as [de2013mathematical]:

$\displaystyle\sum_{i,\alpha}M_{j\beta}^{i\alpha}\Theta_{i\alpha}=\lambda_{1}\Theta_{j\beta},$
(15)

where $\lambda_{1}$ is the largest eigenvalue of $M$ and $\Theta$ its
corresponding eigentensor, whose values represent the centrality of each node
in each layer. The problem of finding the eigenvector centrality consists in
computing
$\Theta_{j\beta}=\lambda_{1}^{-1}\sum_{i,\alpha}M_{j\beta}^{i\alpha}\Theta_{i\alpha}$,
which represents the multilayer generalization of Bonacich’s eigenvector
centrality per node per layer [de2013mathematical, de2015ranking]. By summing
up the scores of a node across all the layers, Eigenvector versatility can be
condensed across layers: $\theta_{i}=\sum_{\alpha}\Theta_{i\alpha}$.

Note that the summation across layers appears naturally through the tensorial
formalism [de2013mathematical, de2015ranking]. Nevertheless, we can opt for
other types of aggregation, based on heuristics specific to the nature of the
centrality vector, to summarize the centrality measures computed for all
layers into a unique descriptor. For other definitions, we refer the
interested reader to [sola2013eigenvector, halu2013multiplex] and
[taylor2021tunable], the latter proposing tunable eigenvector-based
centralities that can be applied to both temporal and multiplex networks.

####  III.3.4 Katz versatility

In the case of directed networks, nodes with only outgoing edges have by
definition eigenvalue centrality equal to zero. This may lead to meaningless
results for the eigenvalue centrality, such as in the limiting case of a
directed acyclic network, where all nodes have zero centrality
[newman2018networks]. This problem can be overcome by assigning to each node a
minimum value of centrality, so that all nodes are taken into account for
measuring the influence between neighbours [katz1953new]. This can be realized
[de2013mathematical] by redefining the eigenvalue problem as:

$\displaystyle\sum_{i,\alpha}(aM_{j\beta}^{i\alpha}\Phi_{i\alpha}+1)=\Phi_{j\beta}$
(16)

with $a<1/\lambda_{1}$, being $\lambda_{1}$ the largest eigenvalue of $M$. As
for the Eigenvector versatility, the overall Katz centrality of a node is the
sum of centrality scores $\Phi_{j\beta}$ across layers [de2013mathematical,
de2015ranking].

####  III.3.5 HITS versatility

In single-layer networks, the Hyperlink-Induced Topic Search or HITS
centrality (also known as hubs and authorities) was originally introduced for
web page rating according to their authority (e.g. their content) and their
hub value (e.g. the value of their links to other web-pages)
[kleinberg1999authoritative]. Again in the case of directed graphs, it can be
useful to recognise as important a node that is pointed to by important nodes
or, alternatively, a node that points to important nodes. These two behaviours
define two different roles a node can play in a directed network: hub and
authorities. In the case of multilayer networks, how extensively a node plays
the role of a hub or an authority is measured by the HITS versatility, that is
defined by two eigenvalue problems [de2013mathematical, de2015ranking]:

$\displaystyle\sum_{i,\alpha}(MM^{\dagger})_{j\beta}^{i\alpha}\Gamma_{i\alpha}=\lambda_{1}\Gamma_{j\beta}$
(17)
$\displaystyle\sum_{i,\alpha}(M^{\dagger}M)_{j\beta}^{i\alpha}\Upsilon_{i\alpha}=\lambda_{1}\Upsilon_{j\beta},$
(18)

where $\lambda_{1}$ is the largest eigenvalue (which is the same for the two
problems), $\Gamma$ represents the authority centrality, and $\Upsilon$, the
hub centrality.

####  III.3.6 Random walk occupation centrality

A random walk on a network is a stochastic process where a path is defined,
starting from a node of origin $X(0)$, with the node $X(t+1)$ chosen at random
from among the neighbours of node $X(t)$. For a multilayer network, the
probabilities of transition between pairs of nodes are represented as a tensor
$T_{j\beta}^{i\alpha}$ [de2013mathematical], that are often taken as
proportional to the edges’ weights. If we let $p_{i\alpha}(t)$ be a time
dependent tensor giving the occupation probability of a given node in a given
layer at time $t$, the random walk can be modelled as a Markov chain (note
that we are using the Einstein convention):

$\displaystyle p_{j\beta}(t+1)=T_{j\beta}^{i\alpha}p_{i\alpha}(t)$ (19)

The steady state for this equation $\Pi_{i\alpha}$ can be obtained as the
leading eigentensor in the eigenvalue problem:

$\displaystyle T_{j\beta}^{i\alpha}\Pi_{i\alpha}=\lambda_{1}\Pi_{j\beta}$ (20)

The probabilities $\Pi_{i\alpha}$ define the random walk occupation centrality
[sole2016random], a measure that highlights which nodes might experience
congestion due to insufficient outflow.

####  III.3.7 PageRank versatility

For single-layer networks, PageRank [page1999pagerank] is a centrality measure
originally developed for ranking webpages. It represents the occupation
probability of random walkers on the network subjected to teleportation: at
any time the walker can walk to a neighbour with a rate $r$ and be teleported
to any other node with rate $1-r$. Its extension to a multilayer network can
be based on different heuristics, in case of unknown nature of the layer’s
coupling [halu2013multiplex], or on the tensorial formulation
[de2013mathematical, de2015ranking] if the inter-layer connectivity is known.
In this second case, the dynamics of the random walker are regulated by a
transition tensor defined by:

$\displaystyle
R_{j\beta}^{i\alpha}=rT_{j\beta}^{i\alpha}+\frac{1-r}{NL}u_{j\beta}^{i\alpha},$
(21)

where $T_{j\beta}^{i\alpha}$ is the transition tensor of a classical random
walk in the absence of teleportation and $u_{j\beta}^{i\alpha}$ is the rank 4
tensor with all components equal to 1. If we denote by $\Omega_{i\alpha}$ the
eigentensor of $R_{j\beta}^{i\alpha}$, then the PageRank versatility
$\omega_{i}$ of a node is obtained by summing across layers
$\sum_{i,\alpha}\Omega_{i\alpha}$. In other words, the Page Rank centrality
for multilayer networks is the steady-state solution of the master equation
corresponding to the transition tensor $R_{j\beta}^{i\alpha}$.

It is important to remark that, for all nodes without out-going edges, the
transition tensor has to be redefined as
$R_{j\beta}^{i\alpha}=\frac{1}{NL}u_{j\beta}^{i\alpha}$ to grant the correct
normalisation of the transition probabilities. We refer to [halu2013multiplex,
battiston2016efficient, iacovacci2016functional] for existing variants on
PageRank definitions in the context of multiplex networks.

####  III.3.8 Random walk betweenness versatility

The random walk betweenness versatility [sole2016random] measures the
importance of nodes in terms of the number of times random walk paths in the
network pass by a given node. This quantity is suitable, for instance, for
identifying the critical nodes in the random spreading of pieces of
information, which are not necessarily following the shortest trajectories
[freeman1991centrality]. To analytically compute the random walk betweenness
versatility, it is convenient to take advantage of the concept of absorbing
random walks, i.e. walks that will end when a given node $d$ is reached. These
walks are defined by the absorbing transition tensor on a particular node $d$:

$\displaystyle(T_{[d]})^{i\sigma}_{j\beta}=\left\\{\begin{array}[]{l l}0&\quad
j=d\\\ T^{i\sigma}_{j\beta}&\quad j\neq d\end{array}\right.$ (24)

The average number of times a random walk (originating in node $o$ in layer
$\sigma$ and destination $d$, independently of the layer) will pass by a node
$j$ in layer $\beta$, regardless of the time step, is given by
[sole2016random]:

$\displaystyle(\tau_{[d]})^{o\sigma}_{j\beta}=\left[(\delta-
T_{[d]})^{-1}\right]^{o\alpha}_{j\beta},$ (25)

where $\delta^{i\alpha}_{j\beta}=\delta^{i}_{j}\delta^{\alpha}_{\beta}$ and
$\delta$ is the Kronecker delta. The average over all possible starting layers
$\sigma$ and the aggregation of the walks that pass through $j$ in the
different layers is obtained by

$\displaystyle(\tau_{[d]})^{o}_{j}=\frac{1}{L}(\tau_{[d]})^{o\sigma}_{j\beta}u^{\beta}u_{\sigma}.$
(26)

Finally, the overall betweenness versatility is given by averaging over all
possible origins and destinations:

$\displaystyle\tau_{j}=\frac{1}{N(N-1)}\sum_{d=1}^{N}(\tau_{[d]})^{o}_{j}u_{o}$
(27)

####  III.3.9 Random walk closeness versatility

Random walk closeness centrality [sole2016random] of a node $i$ is defined as
the inverse of the average number of steps required by a random walker,
starting from any other node in the multilayer network, to reach $i$ for the
first time. As for the case of betweenness, by using the concept of absorbing
random walks one can find that for walks originating in node $o$ in layer
$\sigma$ the probability of visiting node $j$ in layer $\beta$ after $t$ time
steps is given by $(T_{[d]})^{o\alpha}_{j\beta}$ to the power of $t$.
Considering the walk starting from node $o$ in layer $\sigma$, each tensor
encoding the mean first passage time to node $d$ is given by [sole2016random]:

$\displaystyle(H_{[d]})^{o\sigma}=\left[(\delta-
T_{[d]})^{-1}\right]^{o\sigma}_{j\beta}u^{j\beta}.$ (28)

The average mean first passage time $h_{[d]}$ to node $d$ is obtained by
averaging $(H_{[d]})^{o\sigma}$ over all possible starting nodes and layers
as:

$\displaystyle
h_{[d]}=\frac{1}{(N-1)L}u_{o\sigma}(H_{[d]})^{o\sigma}+\frac{1}{N}\Pi^{-1}_{[d]},$
(29)

where $\delta^{i\alpha}_{j\beta}=\delta^{i}_{j}\delta^{\alpha}_{\beta}$ and
$\Pi_{[d]}$ is the random walk occupation centrality. The last term
$\frac{1}{N}\Pi^{-1}_{[d]}$ is included to account for the average return
times, usually excluded when using absorbing random walks. Finally, the random
walk closeness centrality of node $d$ is simply obtained as the inverse of the
distance $\xi_{d}=1/h_{[d]}$.

####  III.3.10 Betweenness centrality

If a metric distance can be defined in a multilayer topology, it is possible
to extend the classical definitions of node betweenness centrality, edge
betweenness centrality and closeness centrality to obtain their multilayer
counterparts [morris2012transport, magnani2013combinatorial]. Given the
existence of layers, it is possible to select a subset of layers $\Omega$ and
consider only the paths belonging to that subset, defining the _cross
betweenness centrality_ [buccafurri2013measuring] of the node $i$:

$\displaystyle CBC(i,\alpha,\Omega)=\sum_{o\neq i,d\neq
i,\beta\in\Omega,\beta\neq\beta^{\prime}}\frac{\sigma_{o\beta}^{d\beta^{\prime}}(i,\alpha)}{\sigma_{o\beta}^{d\beta^{\prime}}},$
(30)

which counts the fraction of inter-layer shortest paths, having their
destination in $\Omega$, that pass through node $i$ of layer $\alpha$. Taking
advantage of this definition, the multiplex betweenness centrality (BC) can be
decomposed into the following contributions:

$\displaystyle
BC(i,\alpha)=CBC(i,\alpha,\Omega)+CBC(i.\alpha,\bar{\Omega})+IBC(i,\alpha),$
(31)

where $\bar{\Omega}$ indicates the layers not belonging to $\Omega$ and
$IBC(i,\alpha)$ is called _internal betweenness centrality_ and represents the
contribution at the BC of the paths that never leave the layer $\alpha$:

$\displaystyle IBC(i,\alpha)=\sum_{o\neq i,d\neq
i,\alpha=\beta=\beta^{\prime}}\frac{\sigma_{o\beta}^{d\beta^{\prime}}(i,\alpha)}{\sigma_{o\beta}^{d\beta^{\prime}}}.$
(32)

Other definitions are also available and define betweenness versatility as
another natural extension of the monoplex betweenness centrality
[sole2014centrality, SoleRibalta2016congestion].

####  III.3.11 Interdependence centrality

In general, the presence of more than one layer within a system increases the
number of available paths with respect to the case where only one single layer
is present. The richness of multilayer paths allows for the possibility that
single-layer shortest paths are longer than multilayer ones. To quantify the
value added by multiplexity to potential communication routes, it is possible
to define the reachability of each physical node $i$ in terms of
interdependence centrality [nicosia2013growing], defined as:

$\displaystyle\lambda_{i}=\frac{1}{N-1}\sum_{j\neq
i}\frac{\psi_{ij}}{\sigma_{ij}},$ (33)

where $\psi_{ij}$ is the number of shortest paths between nodes $i$ to $j$
that pass by more than one layer and $\sigma_{ij}$ is the total number of
shortest paths from $i$ to $j$. In both cases, the number of shortest paths
can be greater than 1 if there are multiple paths having the same, minimal,
length. Note that these quantities are not tensors or matrices: they can be
better understood in terms of matrix entries, since they encode pairwise
information about physical nodes. The global interdependence is computed by
summing up over the interdependence centrality of all nodes:
$\lambda=\langle\lambda_{i}\rangle$.

In some systems, such as human transportation, the flows between different
pairs of nodes might be strongly dis-homogeneous. In such cases, the global
interdependence can be redefined to take into account the real loads of the
network by weighting the aforementioned sums by $T_{ij}$, which is an origin-
destination matrix whose entries are normalised such that
$\sum\limits_{i,j}T_{ij}=1$. The weighted interdependence, also known as
_coupling_ [morris2012transport], is therefore given by

$\displaystyle\lambda^{\prime}=\sum_{i,j\neq
i}T_{ij}\frac{\psi_{ij}}{\sigma_{ij}}.$ (34)

###  III.4 Mesoscale Organization

####  III.4.1 Triadic Closures

Figure 13:  An elementary cycle on a multilayer network can include edges from
up to three different layers, identifying a triad of physical nodes. The
starting point of the cycle is highlighted in orange, intra-layer edges are
dotted curves and intra-layer edges are solid lines. To simplify the reading
of the cycle, the second intra-layer step is represented by a yellow line.
Above each cycle it is indicated the tensorial form of the cycle in terms of
the tensors $\mathcal{A}$ and $\mathcal{C}$ described in the text. Figure from
[cozzo2015structure].

The clustering coefficient, a measure of the cohesiveness of triads, is an
important topological property of a node and a whole network. Nevertheless, it
can also be defined as the number of 3-cycles that start and end at the same
node [cozzo2015structure], relating this network descriptor to how information
flows and localizes. As exemplified in Fig. 13 for an interconnected multiplex
network, a 3-cycle can go through different layers, but is characterized only
by 3 intra-layer steps. These type of walks, in which after or before each
intra-layer step the walker can choose between continuing on the same layer or
changing to some other adjacent layer, can be explicitly decomposed into
combinations of intra- and inter-layer steps by using adequately defined
tensors:

$\mathcal{A}_{i\alpha}^{j\beta}=\left\\{\begin{array}[]{l
l}M_{i\alpha}^{j\beta}&\quad\alpha=\beta\\\
0&\quad\alpha\neq\beta\end{array}\right.$ (35)

encoding intra-layer steps, and

$\mathcal{C}_{i\alpha}^{j\beta}=\left\\{\begin{array}[]{l
l}M_{i\alpha}^{j\beta}&\quad i=j\\\ 0&\quad i\neq j\end{array}\right.$ (36)

encoding inter-layer steps. Since changing layer is optional, we represent
layer’s choice as

$\mathcal{\hat{C}}_{i\alpha}^{j\beta}=bI_{\alpha}^{\beta}+g\mathcal{C}_{i\alpha}^{j\beta},$
(37)

where $I$ is the identity matrix, $b$ represents the weight of staying in the
layer $i$ and $g$ the weight of stepping from layer $\alpha$ to layer $\beta$.

With this definition, two types of cycles can be defined, and therefore two
different number of cycles starting in $i$ can be obtained. In the first case,
in which two consecutive inter-layer steps are forbidden, this number is
$t_{(w),i\alpha}=2\left((\mathcal{A}\mathcal{\hat{C}})^{3}\right)_{i\alpha}^{i\alpha}$,
corresponding to the examples shown in Fig. 13), while in the second case, if
those steps are permitted, this number is
$t_{(sw),i\alpha}=\left((\mathcal{\hat{C}}\mathcal{A}\mathcal{\hat{C}})^{3}\right)_{i\alpha}^{i\alpha}$.

To compute the multiplex clustering a normalisation factor $d_{(*),i\alpha}$
is also required, where with the label (∗) we indicate both types of cycles,
(w) and (sw). The corresponding definitions for both factors $d$ can be
obtained from the $t$ numbers by replacing the tensor $\mathcal{A}$ in the
second intra-layer step with the tensor:

$F_{i\alpha}^{j\beta}=\left\\{\begin{array}[]{l l}1&\quad i\neq j\ \text{and}\
\alpha=\beta\\\ 0&\quad\text{elsewhere}\end{array}\right.$ (38)

For example, for the (w) type of cycles we have
$d_{(w),i\alpha}=2\mathcal{A}\mathcal{\hat{C}}F\mathcal{\hat{C}}\mathcal{A}\mathcal{\hat{C}}$.
The local multiplex clustering coefficient for the state node $(i,\alpha)$ can
be then calculated as:

$C_{(*),i\alpha}=\frac{t_{(*),i\alpha}}{d_{(*),i\alpha}},$ (39)

while the global multiplex clustering coefficient, the scalar value
representative of triadic closure in the whole system, is given by

$C_{(*)}=\frac{\sum\limits_{i\alpha}t_{(*),i\alpha}}{\sum\limits_{i\alpha}d_{(*),i\alpha}}$
(40)

Figure 14:  Relation between local clustering coefficient, an indicator of
triadic closure, and node’s unweighted degrees for $m$-layers cycles
$C^{(m)}_{M,u}$ as defined in Eq. (41), and for all cycles but considering the
multilayer network with its unweighted $C_{M,u}$ or weighted $C_{Z,u}$
representation. Two empirical systems are considered: a) Kapferer tailor-shop
social network; b) the airport network obtained from openflights.org. The
curves interpolate the local clustering coefficient for each type of cycle.
Figure from [cozzo2015structure].

We can decompose Eq.( 40) into its contributions coming from cycles traversing
a different number $m$ of layers as

$C^{(m)}_{(*)}=\frac{\sum\limits_{i\alpha}t_{(*),i\alpha}^{(m)}}{\sum\limits_{i\alpha}d_{(*),i\alpha}^{(m)}},$
(41)

where $t_{(*),i\alpha}^{(m)}$ and $d_{(*),i\alpha}^{(m)}$ are restricted to
cycles involving exactly $m$ layers, with $m=1,2,3$. This perspective can be
useful to appreciate the differences observed for cycles of a different
nature, evident in Fig. 14, where local clustering coefficients are drawn as a
function of the nodes’ degrees. In several empirical multiplex systems
[cozzo2015structure], it has been observed that $C<C^{(1)}$ and
$C^{(1)}>C^{(2)}>C^{(3)}$ as a general pattern which provides a kind of
hierarchy of triadic closure across layers.

It is worth noticing that triadic closure can also be used to predict future
or missing links on a multilayer network, as recently shown in
[aleta2020link]. As can be seen from Fig. 15, depending on the location of the
links there are four possible triadic relations which allow us to to predict a
missing link by closing a cycle of two links, $u$ on layer $\alpha$ and $v$ on
layer $\beta$, with a third link on layer $x$. This observation allows for a
multilayer generalization of the Adamic-Adar method [adamic2003friends], which
predicts links on the basis of a score based on the number of common neighbors
weighted by their degree.

Figure 15: Four possible types of triadic closures on a multilayer network
that can be used for predicting the existence of missing links. (a)
$\mathcal{T}_{xx}$, if the two links belong to the same layer $x$ where the
cycle is closed. (b) $\mathcal{T}_{x\alpha}$, if one link belongs to $x$ and
another to a different layer $\alpha$. (c) $\mathcal{T}_{\alpha\alpha}$ if
both existing links belong to the same layer, different from $x$. (d)
$\mathcal{T}_{\alpha\beta}$ if the two existing links belong to two separate
layers, both different from $x$. Figure from [aleta2020link].

In fact, in multilayer networks, the neighbors of a node can belong to
different layers. An appropriate “Multilayer Adamic-Adar” (MAA) score counts
the common neighbors closing the triads of each of the aforementioned types,
weighting each contribution by the logarithm of the degree as usually done to
calculate the Adamic-Adar score:

$MAA(u,v)=\sum_{\alpha,\beta}\sum_{w\in\mathcal{T}_{\alpha\beta}}\frac{\eta_{x\alpha}\eta_{x\beta}}{\sqrt{\langle
k\rangle_{\alpha}\langle
k\rangle_{\beta}}}\frac{1}{\sqrt{\ln(k_{w}^{\alpha})\ln(k_{w}^{\beta})}}$ (42)

where $\langle k\rangle_{\alpha}$ indicates the average degree of nodes in
layer $\alpha$, $k_{w}^{\alpha}$ the degree of node $w$ in layer $\alpha$ and
$\eta_{x\alpha}$ are free coefficients allowing to control the relative weight
of each type of triadic closure in the link’s total score.

Fig. 16 shows the AUC (panel a) and the precision (panel b) of MAA for three
empirical systems, as a function of the coefficients $\eta_{x\alpha}$. The
results show that single-layer link prediction can be largely improved by
exploiting information encoded in other layers.

Figure 16: How areas under the curve (AUC) (a) and precision (b) change for
different values of the coefficients $\eta_{x\alpha}$, for three real-world
multiplex networks. For sake of simplicity, in each network only three layers
are considered. In both panels, each column corresponds to a different data
set, from left to right: a social network of students in Copenhagen, the _C.
elegans_ genetic network, and the CKM social network of physicians. Each row
corresponds to a prediction in a different layer. The cross marker indicates
the location of the maximum value for each plot, corresponding to the
combination of coefficients optimizing the indicator. Figure from
[aleta2020link].

####  III.4.2 Communities and modules

Complex systems are characterized by the mesoscale organization of their units
into groups, also known as _modules_ or _communities_ [fortunato2010community,
newman2012communities, fortunato2016community]. However, the definition of
what a group exactly corresponds to is an open problem. Here, we will briefly
describe the major advances in this direction, considering four methods widely
used in the literature, namely: multilayer modularity maximization, multilayer
tensor factorization, multilayer Bayesian inference and multilayer description
length minimization through the map equation. A variety of methods is
available, including the analysis of intermittent communities in time-varying
networks [aslak2018constrained], but they are beyond the scope of this work
and deserve a dedicated review [holme2012temporal, holme2015modern].

##### Multilayer modularity maximization

Figure 17:  Illustration of a _multislice_ network representing, for instance,
a time-varying network where each temporal snapshot is encoded by a layer and
inter-layer connectivity encodes temporal relationships among state nodes.
Figure from [mucha2010community].

The simplest – and also one of the most popular ones – definition concerns the
density of links within a group with respect to inter-group density of links:
for a module, we expect units to be mostly connected with other units inside
the module and poorly connected outside. This approach is based on the
calculation of a function named _modularity_ which, for a given partitioning
of the units into groups, quantifies the deviation of the number of links from
what is expected by chance according to the corresponding configuration model
[newman2006modularity]. This function is therefore calculated for all possible
partitions, and the partition with maximum modularity is chosen as the most
representative mesoscale organization of the underlying network, providing a
hint about its structural and functional organization. For this reason,
identifying such an organization is of paramount importance for applications
in many disciplines, from social sciences to biology [girvan2002community,
guimera2005functional].

Figure 18:  Evolution over time of political groups, identified by multilayer
modularity maximization, in the U.S. Senate across more than two centuries
(110 layers encoding the number of 2-year Congresses from 1789 to 2008).
Intra-layer links are given by roll call vote similarities and inter-layer
connectivity is weighted by $\omega=0.5$. (A) Results from senators’
assignments. (B) Results from state affiliations. Note that less sophisticated
analytical techniques, applied for instance on aggregated representation of
this system, would not identify the same patterns. Analysis and Figure from
[mucha2010community].

The definition of the multilayer generalization of modularity was given by
[mucha2010community] a decade ago, become a standard for applications
[carchiolo2011communities, amelio2019modularity]. It assumes that the system
can be represented by a three-dimensional tensor, encoding within-layer
connectivity as for edge-colored multigraphs and time-varying networks666Where
each snapshot is encoded by a layer. (see Fig. 17) and by another three-
dimensional tensor encoding inter-layer connectivity. To allow for a multi-
resolution analysis777We refer the interested reader to
[fortunato2007resolution, arenas2008analysis, Taylor2017superresolution] for
more information about the resolution and super-resolution problems in
community detection. Note that recent studies have shown that by aggregating
layers by summation, for instance, it is also possible to enhance the
detectability of community structures [Taylor2016detectability]., a resolution
parameter $\gamma$ is introduced, and inter-layer connectivity is weighted by
another tunable parameter $\omega$. Using the tensorial notation, considering
the more general systems that can be represented by the multilayer adjacency
tensor $M^{i\alpha}_{j\beta}$, and indicating by $P^{i\alpha}_{j\beta}$ the
tensor encoding the connectivity values obtained from a null model (e.g., the
configuration model), one can define the rank-4 tensor

$\displaystyle
B^{i\alpha}_{j\beta}=M^{i\alpha}_{j\beta}-P^{i\alpha}_{j\beta},$ (43)

which is then used to define the modularity function

$\displaystyle Q\propto S^{a}_{i\alpha}B^{i\alpha}_{j\beta}S^{j\beta}_{a}$
(44)

to be maximized while varying the membership tensor $S^{i\alpha}_{a}$, with
entries equal to 1 if node $i$ in layer $\alpha$ belongs to the community $a$,
and zero otherwise [de2013mathematical]. This approach can be suitably used to
identify groups in time-varying systems, like the U.S. Senate (Fig. 18). We
refer to [bazzi2016community] for an enhanced version of this method and to
the recent work demonstrating that, under certain conditions, modularity
maximization corresponds to maximizing the posterior probability of community
assignments under suitably chosen stochastic block models
[pamfil2019relating]. For further information about the latter, see also the
subsection on multilayer Bayesian inference later in this text.

Figure 19:  Illustrating the Kruskal decomposition of a rank-3 tensor
representing a time-varying network. The tensor is factorized and decomposed
into the sum of the outer product of rank-1 tensors. Figure from
[gauvin2014detecting].

##### Multilayer tensor factorization

At variance with modularity maximization, this approach is based on
decomposing the three-dimensional tensor $T_{ij\alpha}$, encoding the inter-
layer connectivity of a time-varying network where the order of tensors
follows the arrow of time, by means of non-negative factorization. This
procedure maps a tensor into the sum of outer products of rank-1 tensors as
schematically shown in Fig. 19.

Figure 20:  Non-negative tensor factorization of temporal snapshots (the
layers) describing close-range interactions (the links) of students and
teachers (the nodes), divided into classes. The panels show the distribution
of membership weights encoded by the first output of the procedure. Figure
from [gauvin2014detecting].

The result of this procedure is a coarse-grained representation of the system
in terms of two assignments: nodes to groups, giving the membership weight of
units, and time intervals to groups, giving the activity level of each group
at different temporal snapshots. When applied to a contact network, for
instance, the result of this approach is shown in Fig. 20.

##### Multilayer Bayesian inference

Figure 21:  Summarizing a non-parametric generative process in the case of a
monoplex network, according to degree-corrected stochastic block modeling
(DCSBM): (a) sampling the partition of units; (b) counting edges between and
within groups; (c) imposing the observed degree sequence; (d) sampling the
overall network by accounting for (a–c). Figure from
[peixoto2017nonparametric].

A different approach is to consider the problem of recovering the best
partition of the nodes, where one assumes that there is an underlying model
encoding some mechanisms at work to generate the observed system. This class
of methods, originally developed for analysis of social networks
[holland1983stochastic, snijders1997estimation, nowicki2001estimation],
assumes that groups can be encoded by blocks: probability to have links within
and between blocks, as well as the number of blocks, are the parameters of the
model to be calculated. This model is known as the Stochastic Block Model
(SBM), and it finds extensive applications in machine learning
[airoldi2008mixed, goldenberg2010survey, qin2013regularized,
anandkumar2014tensor]. More recently, a variation based on Bayesian inference
and statistical physics has been proposed for applications to multilayer
networks [peixoto2015inferring, valles2016multilayer]. A thorough discussion
of the mathematical framework behind this powerful approach is beyond the
scope of this section, therefore we refer the interested reader to an accurate
and recent work fully devoted to this purpose [peixoto2019chapter].

Figure 22:  Generative models of multilayer networks. Left: generate the
aggregated (collapsed) representation of the system and then create layers
whose aggregation leads to the same representation. Right: each layer is
generated independently from the others. Figure from [peixoto2015inferring].

Here, we outline the basic idea behind this method [peixoto2015inferring],
which is schematically summarized in Fig. 21 for the case of a monoplex. The
input is considered an edge-colored multigraph or a time-varying network that
can be represented by an array of matrices, similar to the ones we have
described in the previous sections. Let us indicate this object with
$M_{ij}^{\alpha}$, and its aggregation obtained by entry-wise summation across
layers as $A_{ij}=M_{ij}^{\alpha}u_{\alpha}$. Let $\boldsymbol{\Theta}$
indicate the set of all parameters of the model to fit. The Bayes theorem in
this case reads

$\displaystyle
P(\boldsymbol{\Theta}|M_{ij}^{\alpha})=\frac{P(M_{ij}^{\alpha}|\boldsymbol{\Theta})P(\boldsymbol{\Theta})}{P(M_{ij}^{\alpha})},$
(45)

where $P(\boldsymbol{\Theta})$ is the prior probability on the parameters,
$P(M_{ij}^{\alpha})$ is a normalization factor and
$P(M_{ij}^{\alpha}|\boldsymbol{\Theta})$ is the likelihood of observing the
multiplex system encoded by $M_{ij}^{\alpha}u_{\alpha}$ given the parameters.
Here, $P(\boldsymbol{\Theta}|M_{ij}^{\alpha})$ is the posterior likelihood:
the higher its value the more likely our model describes the data. By noticing
that:

  * •

$P(\boldsymbol{\Theta})=e^{-\mathcal{L}(\boldsymbol{\Theta})}$, where
$\mathcal{L}(\boldsymbol{\Theta})$ is the microcanonical entropy of the
parameter ensemble;

  * •

$P(M_{ij}^{\alpha}|\boldsymbol{\Theta})=e^{-\mathcal{S}(M_{ij}^{\alpha})}$,
where $\mathcal{S}(M_{ij}^{\alpha})$ is the microcanonical entropy of the
model ensemble;

  * •

the term $P(M_{ij}^{\alpha})$ can be neglected since it does not depend on
parameters and only acts as a constant;

it is straightforward to show888Hint: take the natural logarithm of both sides
of the Bayes formula, discard the constant term and solve by $\Sigma$. that
minimizing the function

$\displaystyle\Sigma=\mathcal{S}(M_{ij}^{\alpha})+\mathcal{L}(\boldsymbol{\Theta}),$
(46)

known as the _description length_ , is equivalent to maximizing the posterior
probability. Note that $\Sigma$ has a very nice interpretation in terms of
information theory: it is the number of bits required to describe the data
summed to the number of bits required to describe the model.

Figure 23:  A non-interconnected multiplex social network of physicians is
reproduced by means of two distinct generative models. The layout is specially
designed to highlight the presence of hierarchies, indicated by red markers
and links, of nodes – the circular dots arranged on a circle, linked by
observed edges drawn by means of bundling to better put in evidence inter-
group connections. (a) the inferred degree-corrected SBM is directly applied
to the aggregated network; (b) inference is performed while accounting for
edge covariates, by assigning each layer to a type of social interaction.
Below each panel, the value of the posterior odds ratio is shown. Figure from
[peixoto2015inferring].

The explicit definition of the prior probability $P(\boldsymbol{\Theta})$ is
often a subtle issue involving the choice of a generative process for the
parameters. This is a problem-specific task which depends on the a priori
knowledge and assumptions about the data. In a more general setting, and in
those cases in which the prior knowledge is missing, it is desirable to choose
uninformative priors. In [peixoto2015inferring] the author proposes a nested
uninformative prior which minimizes the influence on the posterior.

Figure 24:  A proximity network of high-school students, varying over time, is
considered. The probability of existence of an edge is shown, across time, in
the top panel: the colored areas, labeled from (a) to (j), are the ones
obtained from the inference process as the best partition into layers. For
each area, the corresponding partitioning is shown in the bottom panels by
using the hierarchical layout. Figure from [peixoto2015inferring].

One can choose different models to generate the observed network and,
specifically, plausible options are to either generate layers conditional to
the aggregate representation of the system or to generate each layer
independently from the others (Fig. 22), as well as other generative models –
such as the M-DCSBM – specifically designed for multilayer systems
[bazzi2020framework]. The reader could be concerned about how to choose the
best model: in fact, the framework also allows for model selection and one
should take the model which minimizes the description length. When one needs a
more refined approach, where a model or its alternative should be chosen, non-
parametrically, with some degree of confidence, it is possible to calculate
posterior odds ratio as

$\displaystyle\Lambda=\frac{P(\boldsymbol{\Theta}_{a}|M_{ij}^{\alpha},\mathcal{H}_{a})P(\mathcal{H}_{a})}{P(\boldsymbol{\Theta}_{b}|M_{ij}^{\alpha},\mathcal{H}_{b})P(\mathcal{H}_{b})}=e^{-\Delta\Sigma}\frac{P(\mathcal{H}_{a})}{P(\mathcal{H}_{b})},$
(47)

where each model class is encoded into the hypothesis $\mathcal{H}_{p}$
($p=a,b$), $P(\mathcal{H}_{p})$ is the prior belief for that hypothesis and
$\Delta\Sigma=\Sigma_{a}-\Sigma_{b}$ is the difference between the description
length of each hypothesis. The operational prescription is that if $\Lambda<1$
then $\mathcal{H}_{b}$ should be preferred to $\mathcal{H}_{a}$, whereas when
$\Lambda\approx 1$ both models are equally plausible. More generally, if
$\Lambda\in[1/3,1]$ the difference between the two hypotheses can be
considered negligible. An application of this procedure is shown in Fig. 23 in
the case of an empirical social network.

Finally, this procedure can be used to identify the best division of time-
varying edges into layers, as shown in the case of a real-world contact
network in Fig. 24.

We refer the interested reader to a more recent work [pamfilinference], where
an SBM-like model – not relying on sampling of edges in different layers
independently – is proposed to account for edge correlations, and to define a
new measure of layer-layer correlation which incorporates similarity between
connectivity patterns in different layers. It is used for link prediction in
multilayer networks.

##### Multilayer description length minimization through the map equation

Figure 25:  Detection of modules based on the Infomap algorithm. (A) The solid
line represents a trajectory of random walks on the top of a classical
network. (B) By encoding each node with a codeword, e.g. using a Huffman
scheme, one can efficiently describe the trajectories in terms of strings by
appending the codeword of each node visited by the walker. The trajectory
shown in (A) is encoded into a string of 314 bits. (C) Given a partitioning of
the system, the persistence of the trajectory in a certain module is encoded
by a certain number of bits: the lower this number the more persistent the
trajectory in that module. This regularity can be used to coarse-grain the
network into functional modules as in (D). Figure from [rosvall2008maps].

Instead of relying only on topological information – such as the abundance of
edges or the presence of blocks – to identify groups, one can observe how
information flow is trapped within modules, the rationale being that
information is more likely to flow within a group than between groups
[pons2006computing, lambiotte2014random]. Another approach might be to map the
potential regularities of the flow into strings and to look for the
partitioning of the system that minimizes the corresponding description length
[rosvall2007information, rosvall2008maps, esquivel2011compression], as
schematically described in Fig. 25.

Figure 26:  A complex system is analyzed using a conventional approach (a) and
the higher-order framework (b) provided by the map equation. The latter allows
us to account for the richness of the data, such as weights and
directionality, without being affected by limited detectability (e.g., due to
the resolution problem mentioned for modularity maximization) and, remarkably,
is suitable to detect multilevel and overlapping modular information within
the same elegant mathematical framework. Figure from [edler2017mapping].

The generalization to multilayer networks of this method, known as Multiplex
InfoMap, was introduced in [de2015identifying] and, more recently, it was
better formalized under the perspective of higher-order modeling
[rosvall2014memory, salnikov2016using, edler2017mapping,
lambiotte2019networks], as schematically illustrated in Fig. 26.

Figure 27:  Toy multilayer system with 4 physical nodes and 3 layers,
represented in terms of structure (a-c) and dynamics (b-d) with a layered
perspective (a-b) and with respect to state-node representation (c-d), the
latter highlighting the existence of memory (encoded by small colored nodes)
which alters the information flow (solid, colored, lines). Figure from
[de2015identifying].

Multiplex InfoMap is based on the same principles of its single-layer
counterpart, with some important differences. First, the monoplex InfoMap is
based on Markovian dynamics, whereas this feature is kept only for within-
layer dynamics. Second, coding only captures the visits to physical nodes, not
the ones to state nodes, an effective non-Markovian dynamics across layers.
The remainder of the method is the same: compressing persistent multilayer
trajectories999The curious reader might wonder about the governing equations
of the corresponding random walk: this will be discussed in the chapter
devoted to dynamics and, specifically, in Sec. IV.1. to identify modules
according to the map equation

$\displaystyle\Sigma(\mathcal{M})=q_{in}\mathcal{S}(\mathcal{Q})+\sum\limits_{\ell=1}^{m}p_{\ell}\mathcal{S}(\mathcal{P}_{\ell})$
(48)

which encodes information flows within and across layers, as shown in Fig. 27.
This equation deserves a careful explanation of its terms.

Here, $\mathcal{M}$ is a partitioning of the system – i.e., a model for the
observed modules – and $\Sigma$ is its description length, in bits. Indicating
with $q_{\ell,in}$ and $q_{\ell,out}$ the transition rates at which a random
walker enters and exits, respectively, a module $\ell$, with
$q_{in}=\sum_{\ell}q_{\ell,in}$ the sum of those rates and with
$\mathcal{S}(\mathcal{Q})$ the information entropy of the normalized
probability distribution of the transition rates
($Q=\\{q_{\ell,in}/q_{in}\\}$), the first term in the right-hand side of the
map equation captures the number of bits required to describe the dynamics.
Indicating with $p_{l\in\ell}$ the visit rates of the physical nodes for the
module codebook $\ell$, with $p_{\ell}$ the sum of those rates and with
$\mathcal{S}(\mathcal{P}_{\ell})$ the information entropy of the corresponding
normalized probability distribution
($\mathcal{P}_{\ell}=\\{p_{l\in\ell}/p_{\ell}\\}$), the second term of the map
equation captures the number of bits required for coding. An application of
Multiplex InfoMap to a toy system is schematically summarized in Fig. 27.

The interested reader can play interactively with an online version of this
method, which is publicly
available101010http://www.mapequation.org/apps/sparse-memory-
network/index.html..

####  III.4.3 Integrated and segregated systems

To understand how a network operates, it is essential to understand if
information flows in such a way that, from a functional perspective, the
system is either integrated, i.e., operating like a whole, or segregated,
i.e., operating like independent modules or groups of units (see Fig. 28).

Figure 28: Segregated communities and integrating hubs. (a) Network with four
communities interconnected by four central hubs. (b) Functional segregation
due to community structure is estimated by within-community coupling. (c)
Functional integration favors information flow across the whole network
through a subset of interconnecting hubs and is measured by the hubs or the
paths crossing the network. Figure from [sporns2013network].

How to correctly study integration and segregation of a complex network is a
still debated topic, and several proxies have been proposed across a broad
spectrum of disciplines. The analysis stem almost independently from sociology
and neuroscience, in both cases within the frameworks of classical single-
layer networks [latora2001, latora2003economic, achard2007efficiency,
rubinov2010complex, centola2015social, louf2016patterns, yamamoto2018impact,
gallotti2019disentangling, bertagnolli2021quantifying].

In multilayer systems, where multiple relationships co-exist simultaneously,
the evaluation of integration must be accounted for by more complex
topological models. In fact, it is critical to choose which kind of paths are
to be evaluated, e.g., by adding a cost to inter-layer links used in paths
that allow us to switch between two or more layers. The literature on this
topic is rather poor and an agreement on the most suitable methodology to
adopt is far from being reached. In the following, we describe two recent
attempts in this direction.

Figure 29: Multilayer modeling of structure and function of urban systems. (A)
Urban backbone of 10 megacities described by their street networks. (B) Urban
functional networks described by the flows between cells of size
500m$\times$500m. The Foursquare data used for this analysis allow us to
disentangle mobility flows into different layers where subsequent check-ins
have been made between activities of the same type (e.g., leisure, eating,
shopping, sport, etc.). (C) Measured flows for New York City represented as
the weighted edges, which may be very different across distinct activity
layers. Figure from [gallotti2019disentangling].

The first application comes from the analysis of socio-technical systems, and
in particular movement flows. One important approach in this direction is
communicability [akbarzadeh2018communicability], for which an explicit
multilayer definition is available [Estrada2014Communicability]. However, the
simultaneous analysis of integration and segregation can provide a more
comprehensive perspective of a city’s structural form and functional behavior.
In a city, distinct layers can represent different activities (e.g., leisure,
eating, shopping, sport, etc.; see Fig. 29) that can be performed in the same
area: two areas are connected within or across layers if there is a human flow
between them in either direction. The analysis of empirical flow networks of
this type measured from 10 megacities around the world would reveal, for
instance, which activity was most influential for the urban system, from an
integration-segregation perspective [gallotti2019disentangling]. Using
modularity as a proxy for integration and normalised global efficiency as a
proxy for integration, allowed scholars to unravel the functional organization
of each megacity and to discover the presence of distinct “cities within a
city”. By comparing the deviations measured after aggregating all layers into
a single-layer network and after aggregating while excluding one layer at a
time, it can be shown that “transportation” is the activity whose disruptions
would remarkably hinder the urban function: in fact, the integration granted
by long range connections would be sensibly reduced, while segregation due to
the interrupted connections characterised by large flows would increase.
Conversely, the same analysis shows that activities such as restoration,
leisure and going to shopping malls poorly contribute to a city’s integration
[gallotti2019disentangling].

The fact that this result aligns with our expectations of macroscopic and
easy-to-interpret systems confirms that the methodologies discussed in this
section are suitable for characterising a multilayer network in terms of
integration and segregation. In the specific case of layered urban systems,
this type of analysis is potentially relevant to support policy makers with
quantitative information on which activities can be temporary limited or
promoted to achieve a desired amount of human flows integrated across the
city.

The second application concerns the human brain. In fact, a possible way to
analyze its functional connectivity is to map relationships between distinct
regions of interest across distinct frequency bands [de2016mapping,
de2017multilayer, DeDomenico2018], each one encoding a distinct layer, and we
can wonder whether these layers are integrated or segregated
[tewarie2016integrating]. However, since the metric used to this aim is edge
overlap, according to the aforementioned prescriptions this application is not
entirely compatible with the study of integration and segregation and would be
better understood as an analysis of layer-layer correlations. It will be
interesting to see more studies on this topic in the future, for instance,
based on multilayer modularity and multilayer generalization of global
efficiency.

###  III.5 Layer-layer correlations

Single layer networks usually exhibit correlations between entities acting in
the system. Specifically, these correlations occur between some specific
nodes’ properties. In the case of assortativity, for instance, there is a
preference for network’s nodes to attach to others that are similar in some
way (e.g. in degree), while it is the opposite in the case of
disassortativity. Expanding the concept of _correlation_ to multilayer
networks has to take into account the increased degrees of freedom introduced
by the multilayer structure, where nodes’ involvement across layers exhibits
nontrivial and more variegated patterns than those observed in the single
layer [nicosia2015measuring]. In fact, one can study the degree-degree
correlations of each layer in the network but, in the case of a multilayer
network, it is far more intriguing, for example, to explore how a given
property of a node at a certain layer is correlated to the same or other
properties of the same node at another layer. In the following we present the
main correlation measures designed for multilayer networks, divided into three
main sub-sections: interlayer degree correlation, overlap and degree of
multiplexity and pairwise multiplexity, although other interesting measures –
e.g., based on SBM-like models used to infer edge correlations and for link
prediction [pamfilinference], on estimating a joint probability distribution
describing edge existence over all layers to quantify correlations through
conditional mutual information [wu2020correlated] or on a set-theoretic
approach to quantify if a layer correlates with a second layer directly or via
the indirect mediation with a third layer [lacasa2021beyond]– are available.

####  III.5.1 Inter-layer degree correlation

This kind of correlation points out if high-degree nodes in one layer maintain
this property in other layers. We can measure the inter-layer assortativity by
studying the conditional degree distribution
$P(k^{\alpha^{\prime}}|k^{\alpha})$[nicosia2013growing], where
$k_{i}^{\alpha}$ is the degree of node $i$ at layer $\alpha$, evaluating the
average degree $\bar{k}^{\alpha^{\prime}}$ at layer $\alpha^{\prime}$ of nodes
having degree $k^{\alpha}$ at layer $\alpha$:

$\displaystyle\bar{k}^{\alpha^{\prime}}(k^{\alpha})=\sum_{k_{\alpha^{\prime}}}k^{\alpha^{\prime}}P(k^{\alpha^{\prime}}|k^{\alpha});$
(49)

if there is no correlation between the layers $\alpha$ and $\alpha^{\prime}$
we expect $\bar{k}^{\alpha^{\prime}}(k^{\alpha})=\langle
k^{\alpha^{\prime}}\rangle$, i.e. the average degree
$\bar{k}^{\alpha^{\prime}}$ does not depend on $k^{\alpha}$. If, instead,
$\bar{k}^{\alpha^{\prime}}(k^{\alpha})$ is an increasing function in
$k^{\alpha}$ the degrees have an assortative correlation, while they have a
dissortative correlation if the function is decreasing. The inter-layer
assortativity could also be defined using the Pearson, Spearman or Kendall
correlation between the degrees of the same node in the different layers
[boccaletti2014structure].

####  III.5.2 Overlap and degree of multiplexity

In this case correlation is evaluated in terms of node connectivity patterns
in different layers. Specifically, we refer to the correlations among the
links of the different intra-layer networks. In fact, the internal
connectivity in different layers of the multiplex can be in certain cases
correlated. To clarify with an example, you can be a friend of the same person
in two different social networks, and thus the link between you and your
friend exists in both layers of an imaginary multilayer network representing
your social connections. Many different measures have been proposed to
quantify this topological similarity. For the unweighted case, a first,
associated to a couple of layers, is called global overlap
[bianconi2013statistical], and is defined as the total number of pairs of
nodes connected at the same time by a link in both layers, or, in other words,
the total number of links that are in common between layer $\alpha$ and layer
$\alpha^{\prime}$:

$\displaystyle
O^{\alpha,\alpha^{\prime}}=\sum_{i<j}M_{j\alpha}^{i\alpha}M_{j\alpha^{\prime}}^{i\alpha^{\prime}}$
(50)

(for other variants on overlapping please refer to [battiston2014structural]).
A second quantity, associated this time to the whole multiplex, is the degree
of multiplexity [kivela2014multilayer], defined as the fraction of node pairs
that have multiple edge types between them. This quantity is obtained by
dividing the number of node pairs that have multiple edge types between them
by the total number of adjacent node pairs [kivela2014multilayer]. For the
specific case of weighted multilayer networks we refer to
[menichetti2014weighted] and to [boccaletti2014structure].

In particular, here we present two main weighted measures of multiplex
networks, i.e. multistrength and the inverse multi-participation ratio
[menichetti2014weighted]. The multistrenght for a node $i$ in layer $\alpha$,
indicated by $s_{i,\alpha}^{\bar{m}}$, is obtained by summing the weights of a
certain type of multilink $\bar{m}$ – i.e., the set of links connecting a
given pair of nodes in the different layers of the multiplex – incident to a
single node:

$\displaystyle
s_{i,\alpha}^{\bar{m}}=\sum_{j=1}^{N}a_{ij}^{\alpha}A_{ij}^{\bar{m}}.$ (51)

Given that for every layer $\alpha$ and node $i$ the non trivial
multistrengths must include multilinks $\bar{M}$ with $\bar{M}_{\alpha}=1$,
the number of non trivial multistrengths $\alpha$ $s_{i,\alpha}^{\bar{m}}$ is
given by $2^{M-1}$. It follows that the number of multistrengths that can be
obtained for each node in a multiplex network of $M$ layers is $M2^{M-1}$. The
inverse multiparticipation ratio for a layer $\alpha$ $Y_{i,\alpha}^{\vec{m}}$
is used to measure the heterogeneity of the weights of multilinks incident
upon a single node:

$\displaystyle
Y_{i,\alpha}^{\vec{m}}=\sum_{j=1}^{N}\left(\frac{a_{ij}^{\alpha}A_{ij}^{\bar{m}}}{\sum_{r}a_{ir}^{\alpha}A_{ir}^{\bar{m}}}\right)^{2}.$
(52)

Within the context of a weighted multiplex network, the importance of
considering the interacting layers in the analysis of complex system has been
proven by [menichetti2014weighted] which evaluated the additional amount of
information provided by weighted properties of multilinks over the one
contained in single layers, further strengthening the evidence that the
analysis of multilayer networks cannot be confined to the partial analysis of
single layers.

Figure 30:  (a) Airline, (b) APS and IMDb networks. (c) Top 20 airlines in
Europe by number of covered airports. Figure from [nicosia2015measuring].

####  III.5.3 Pairwise multiplexity

In this last case, correlation is provided in terms of correlated activity
patterns in the multilayer. A node $i$ is said to be active at a layer
$\alpha$ if $k_{i}^{\alpha}>0$. For each node $i$, we can associate a node
activity vector $\mathbf{b_{i}}=\\{b_{i}^{[1]},b_{i}^{[2]},...,b_{i}^{[L]}\\}$
where $b_{i}^{[\alpha]}=1$ if node $i$ has at least one edge at layer $\alpha$
and is 0 otherwise. The node activity $B_{i}$ of the node $i$ corresponds to
the total number of layers where $i$ is active [nicosia2015measuring]:

$\displaystyle B_{i}=\sum_{\alpha}b_{i}^{[\alpha]}.$ (53)

By definition, $0\leq B_{i}\leq L$. Analogously, the layer activity of
$\alpha$ is given by the number of active nodes in layer $\alpha$
[nicosia2015measuring]:

$\displaystyle N^{[\alpha]}=\sum_{i}b_{i}^{[\alpha]}.$ (54)

By definition, $0\leq N^{[\alpha]}\leq N$. We can now define the layer
pairwise multiplexity $Q_{\alpha,\beta}$, which is a measure of correlation
between the layers, as [nicosia2015measuring]:

$\displaystyle
Q_{\alpha,\beta}=\frac{1}{N}\sum_{i}b_{i}^{[\alpha]}b_{i}^{[\beta]},$ (55)

that corresponds to the fraction of nodes that are active in both layer
$\alpha$ and $\beta$. Examples of the distribution of the pairwise
multiplexity for a continental airports networks, for the papers published in
the journals of the American Physical Society (APS) and for the movies in the
Internet Movie Database (IMDb) are reported in Fig. 30. Similarly, we can
define the node pairwise multiplexity, measuring the correlation of activities
between two nodes, as the fraction of layers in which both node $i$ and node
$j$ are active [criado2012mathematical, boccaletti2014structure].

Figure 31: Different types of duplex couplings: uncorrelated, maximally
positive (MP) and maximally negative (MN) correlated multiplexity. Figure from
[lee2012correlated].

We conclude this section showing how the study of the correlation is
fundamental since it alters the critical properties of the network. In fact,
the pattern of correlated multiplexity (Fig. 31) is crucial since it affects
the multiplex system’s connectivity as shown in Fig. 32 where the different
choices of correlations have an impact on the largest connected component.

Figure 32: Size of the giant component of a ER duplex network for (a)
maximally positive correlated multiplexity, (b) uncorrelated and (c) maximally
negative correlated multiplexity. (d–f) Size of the giant component (red) for
$z_{2}=0.4$ and assortativity coefficient $r$ (blue). Figure from
[lee2012correlated].

In this chapter, we have described the most widely used measures for the
topological analysis of multilayer systems, although in some cases we have
used some dynamical process (e.g., random walks, continuous-time diffusion or
shortest-path routing) to define some network descriptors. We refer to
[brodka2018quantifying] for a taxonomy and an experimental evaluation of the
approaches to compare different layers in multiplex networks. The next chapter
will be entirely dedicated to introduce dynamical processes on (and of)
multilayer networks.

##  IV Multilayer dynamics

In the previous chapter we discussed many structural descriptors of a
multilayer network. However, some definitions adopted for the structural
analysis are, in fact, based on some notion of dynamical process on the top of
the network. In this chapter, the reader will have the opportunity to shed
light on a spectrum of important network-driven phenomena, such as diffusive
processes (Sec. IV.1), synchronization dynamics (Sec. IV.2), cooperation
dynamics (Sec. IV.3), intertwined/interdependent processes (Sec. IV.4),
percolation (Sec. IV.5), cascade failures (Sec. IV.6). Note that we refer to
[nicosia2013growing, Santoro2017pareto] and [kivela2014multilayer,
boccaletti2014structure] for further information about the dynamics of
multiplex networks, like growth processes.

Figure 33:  Two classes of dynamical processes on the top of a multilayer
network. Left: single dynamics, e.g. a random walk or the diffusion of some
quantity, is defined for all layers. Right: coupled dynamics, where a distinct
process is defined in each layer and interdependency is specified by means of
coupled governing equations. Reproduced with permission from
[de2020illustrations].

In the following we will distinguish two main classes of dynamics
[de2016physics] on the top of multilayer systems (see Fig. 33):

  * •

Type–I: _single dynamics_ is defined and it is the same for all layers. Here
we can find discrete [de2014navigability, sole2016random, Lacasa2018,
Valdeolivas2018] and continuous-time [gomez2013diffusion, sole2013spectral,
Requejo2016, Brechtel2018, Tejedor2018diffusion] diffusion, pattern formation
[asllani2014turing, asllani2015turing, kouvaris2015pattern,
busiello2018homogeneous], communicable information
[Estrada2014Communicability], synchronization [Singh2015, Skardal2015,
Zhang2015, Saa2018, Leyva2018], epidemic spreading of one disease [Wang2013,
sahneh2013generalized, buono2014epidemics, valdano2015analytical,
bianconi2017epidemic], adoption dynamics, diffusion of innovation and other
contagion processes [Yagan2012contagion, Hu2014viral, ramezanian2015diffusion,
traag2016complex, Wang2017] and opinion dynamics [Jang2015opinion,
Diakonova2016irreducibility, artime2017joint, amato2017opinion,
Antonopoulos2018]. In this class we find also system control [yuan2014exact,
Posfai2016], congestion dynamics [tan2014traffic, SoleRibalta2016congestion,
Chodrow2016demand], cooperation processes [gomez2012evolution,
matamalas2015strategical, battiston2017determinants], percolation
[leicht2009percolation, lee2015towards, kryven2019bond] and cascade failures
[buldyrev2010catastrophic, valdez2020cascading].

  * •

Type–II: _coupled dynamics_ , where multiple dynamical processes are defined
on each layer and their interdependency is operationally encoded into
adequately coupled equations, for instance by means of inter-layer
connectivity when present. In this class we can find processes which combine
epidemics spreading with human behavior [Funk2009, wu2012impact,
lima2013cellular, Granell2013, massaro2014epidemic, lima2015disease, funk2015,
wang2015coupled, azimi2016cooperative, VelasquezRojas2017interacting], simple
and complex contagion ([Czaplicka2016]), evolutionary game dynamics with
social influence [amato2017interplay], cooperative with competitive epidemics
spreading [Dickison2012epidemics, Cozzo2013social, Sanz2014, DeArruda2017,
danziger2019dynamic], transport with synchronization [nicosia2017collective],
interdependent human flows [SorianoPanos2018framework] and other
interdependent processes [GomezGardenes2015competition].

###  IV.1 Diffusive processes

Diffusive processes, either discrete or continuous in time, provide reasonable
models for a wide spectrum of phenomena, from random searches to population
dynamics, from system exploration and navigability to information diffusion.
In the following we will consider two emblematic classes of models: random
walks and continuous-time diffusion.

##### Random walks on single-layer networks

Random walks on networks [noh2004random, masuda2017random] are discrete
processes governed by a transition probability from $i$ to $j$ that, at each
time step, is given by $T_{ij}=\frac{W_{ij}}{s_{i}}$, where $s_{i}$ is the
outgoing strength of node $i$ or, equivalently, $s_{i}=\sum\limits_{j}W_{ij}$.
Let us introduce the matrix diagonal $\mathbf{D}$ whose entries are defined by
$D_{ii}=s_{i}$ and the vector $\mathbf{p}(t)\in\mathbb{R}^{N}$, whose $i$-th
entry gives the probability to find the random walker at time $t$ on node $i$:
the master equation governing the random walk dynamics can be written in
compact form as

$\displaystyle\mathbf{p}(t+1)=\mathbf{p}(t)\mathbf{T},$ (56)

where $\mathbf{T}=\mathbf{D}^{-1}\mathbf{W}$ is known as the transition
matrix. We have just used the standard notation but, for instance, the same
equation can be written using our tensorial formalism as

$\displaystyle p_{j}(t+1)=p_{i}(t)T^{i}_{j},$ (57)

where Einstein convention (sum over $i$) is used and $i,j$ indicate covariant
and contravariant indices, not the entries of the corresponding rank-1 and
rank-2 tensors.

Figure 34:  A random walker visits the nodes of a multilayer network: starting
from one state node in a specific layer, it is allowed to i) stay in the same
state node, ii) jump to a neighbor within the same layer, iii) switch to
another state node on another layer or iv) jump elsewhere while subjected to
the transition tensor $T^{i\alpha}_{j\beta}$. The walk is defined by the
ordered sequence of visited intra- and inter-layer connections. Reproduced
with permission from [de2020illustrations].

The reader familiar with the Markov chains formalism will immediately
recognize the continuous-time Markov chain on $N$ states (i.e., on the set of
nodes), with independent exponential holding times with parameter $\lambda=1$
and jumping probabilities given by the transition probabilities of the
discrete-time random walk, and which can be equivalently defined by the
forward equation

$\dot{\mathbf{p}}(t)=-\mathbf{p}(t)\tilde{\mathbf{L}}$ (58)

where $\tilde{\mathbf{L}}=\mathbf{I}-\mathbf{D}^{-1}\mathbf{W}$ is known as
the random walk normalized Laplacian and $-\tilde{\mathbf{L}}$ is the
generator of the continuous-time process with an initial condition
$\mathbf{p}(0)=\mathbf{p}_{0}$. Equivalently, in tensorial notation we can
rewrite the equation and initial condition as

$\begin{cases}&\dot{{p}}_{j}(t)=-\sum\limits_{i=1}^{N}\tilde{{L}}^{i}_{j}{p}_{i}(t)\\\
&{p}_{i}(0)=q_{i}\end{cases}$ (59)

with $\tilde{{L}}^{i}_{j}=\delta^{i}_{j}-T^{i}_{j}$ indicating the component
$(i,j)$ of the random walk normalized Laplacian tensor, $\delta^{i}_{j}$ the
Kronecker delta and $T^{i}_{j}$ the component $(i,j)$ of the transition
probability tensor. The solution of this equation is given by

$\displaystyle\mathbf{p}(t)=\mathbf{p}_{0}e^{-t\tilde{\mathbf{L}}},$ (60)

where the intervening exponential matrix is known as the _propagator of the
dynamics_. Note that the transition matrix can encode other rules governing
the jumps of a random walker from one node to another: by changing the
_flavor_ of the walk one can explore a broad spectrum of stochastic processes,
as we will see in the next sections, where we will generalize these dynamics
to multilayer networks, where a walker can jump between nodes but can also
switch between layers. In general, the transition matrix is replaced by a
rank-4 transition tensor $T^{i\alpha}_{j\beta}$ which governs the probability
that a walker on node $i$ in layer $\alpha$ will move to node $j$ in layer
$\beta$. An illustration of a random walk on the top of such a system is shown
in Fig. 34, while a schematic representation of the supra-adjacency matrices
corresponding to the multilayer systems considered in the reminder of this
section, is shown in Fig. 9, and is useful to get some visual insights about
the structure of the corresponding supra-transition matrices.

##### Random Walks on edge-colored multigraphs

On this class of systems (note that $M^{i\alpha}_{j\beta}=0$ for all
$\alpha\neq\beta$), it is possible to define random walk dynamics in two
distinct ways. On one hand, we can allow the walker to jump through edges with
different colors, which are considered as multiple edges, regardless of their
colors, where a node’s degree is the sum of its state nodes’ degree
[battiston2016efficient]. Operationally, this approach is equivalent to
performing a classical random walk on the top of the aggregate representation
of the system: in general this could not be desirable, because it is difficult
to quantify the effects of neglected information (i.e., connectivity patterns
per layer) on the dynamics.

On the other hand, we can allow walkers to explore each layers independently
and then integrate the transition matrices governing their dynamics
[de2020EnhanceTransport] as

$\langle
T^{i}_{j}\rangle=\sum\limits_{\alpha=1}^{L}m_{i,\alpha}T_{j\alpha}^{i\alpha},\qquad\sum\limits_{\alpha=1}^{L}m_{i,\alpha}=1\quad\forall
i\in V,$ (61)

where $m_{i,\alpha}\geq 0$ are weights enabling us to tune the importance of
each layer, relative to node $i$ in the final dynamics of the system. If
$m_{i,\alpha}=\frac{1}{L}$ for all $i$ and $\alpha$, then $\langle
T^{i}_{j}\rangle$ is simply the average transition matrix over all layers.
However, if we choose
$m_{i,\alpha}=\frac{1}{\mu_{i}}\mathbf{1}_{\\{s_{i}(\alpha)\neq 0\\}}$, with
$s_{i}(\alpha)$ being the out-strength of vertex $i$ in layer $\alpha$ and
$\mu_{i}=\sum\limits_{\alpha=1}^{L}\mathbf{1}_{\\{s_{i}(\alpha)\neq 0\\}}$ the
multiplicity of node $i$, i.e., the number of layers in which $i$ is not
isolated, then we can discard the effect of $i$ being isolated in one or more
layers. This approach preserves the diversity of information coming from
layers’ connectivity patterns and allows for more flexibility. In fact, in
Equation (61), we use a finite mixture of probability mass functions with
weights $m_{i,\alpha}\geq 0$ such that
$\sum\limits_{\alpha=1}^{L}m_{i,\alpha}=1$, but if additional information from
the data is available then the transition matrices in Eq. (61) can be combined
by means of weights encoding the relative importance given to each layer, thus
providing an operational way to prioritize the available information from
specific layers. Here, for the sake of simplicity, we are just considering
that all layers are equally important and adopting the weights
$m_{i,\alpha}=\frac{1}{\mu_{i}}\mathbf{1}_{\\{s_{i}(\alpha)\neq 0\\}}$.

##### Random walks on multilayer networks

In the following, we indicate by $D(i;\alpha,\beta)$ the inter-layer
connections between replicas across nodes, i.e., the entries of
$M_{j\beta}^{i\alpha}$ for $i=j$. In this framework, the outgoing strength of
node $i$ in layer $\alpha$ is given by
$s_{i}(\alpha)=\sum\limits_{j=1}^{N}M^{i\alpha}_{j\alpha}$, whereas the
multilayer strength – discarding the inter-layer edges – is given by
$s_{i}=\sum\limits_{\alpha}s_{i}(\alpha)=\sum\limits_{\alpha=1}^{L}\sum\limits_{j=1}^{N}M_{j\alpha}^{i\alpha}$.
The inter-layer strengths are similarly given by
$S_{i}(\alpha)=\sum\limits_{j}\sum\limits_{\beta\neq\alpha}M^{i\alpha}_{j\beta}$
and, consequently, $(s_{i}(\alpha)+S_{i}(\alpha))_{i=1,\alpha=1}^{N,L}$ is the
out-strength supra-vector with $NL$ components obtained as the row sums of the
supra-adjacency matrix. The reader will notice that the same results can be
simply obtained by using the tensorial formalism and Einstein convention as
$s^{i\alpha}=M_{j\beta}^{i\alpha}u^{j}u^{\beta}$.

Usually, for the sake of simplicity, isolated nodes and components, as well as
more complex patterns of inter-layer connectivity, are discarded from the
analysis, e.g., focusing on the largest connected component and considering
only nodes which exist in all layers, since a more general framework has a
non-negligible impact on the calculation of transition probabilities. In
general, one can overcome these limitations, typical of real-world systems, by
adding state nodes as isolated and adequately accounting for them in
calculations. In the following, we consider three cases:

  1. 1.

$S_{i}(\alpha)=0$, i.e., inter-layer connectivity is missing for node $i$;

  2. 2.

$s_{i}(\alpha)=0$, i.e., intra-layer connectivity is missing for node $i$ in
layer $\alpha$ (i.e., the node is isolated);

  3. 3.

$s_{i}=0$, i.e., the node is isolated everywhere.

In the last case one can simply exclude the node from the analysis, since it
is not part of the network, in practice. Therefore we can safely assume that
$s_{i}>0$ for all nodes. Concerning the cases 1 and 2, they require special
attention only if $S_{i}(\alpha)+s_{i}(\alpha)=0$ for some $\alpha$, since the
state node $(i,\alpha)$ is an absorbing state, i.e., the walker reaching that
node will stay there with probability equal to 1. To avoid absorbing states of
this type, one option is to add a _teleportation_ dynamics, i.e., a uniform
but small probability equal to $\frac{1}{NL}$ to jump elsewhere in the system
even if outgoing connections are missing. This approach has been successfully
used in Google’s PageRank to rank web pages [brin1998anatomy], as well as to
identify functional modules in monoplex [rosvall2008maps] and multilayer
networks [de2015identifying]. Here, we use a similar approach, with the effect
that the occupation probability of the state node $(i,\alpha)$ is effectively
reduced. Once the transition tensor is well defined, as we will see later in
this section for some specific processes, one can define the master equation
in discrete time as

$\displaystyle p_{j\beta}(t+1)=$
$\displaystyle\underbrace{T_{j\beta}^{j\beta}p_{j\beta}(t)}_{\text{stay}}+\underbrace{\sum_{\begin{subarray}{c}\alpha=1\\\
\alpha\neq\beta\end{subarray}}^{L}T_{j\beta}^{j\alpha}p_{j\alpha}(t)}_{\text{switch}}+\underbrace{\sum_{\begin{subarray}{c}i=1\\\
i\neq
j\end{subarray}}^{N}T_{j\beta}^{i\beta}p_{i\beta}(t)}_{\text{jump}}+\underbrace{\sum_{\begin{subarray}{c}\alpha=1\\\
\alpha\neq\beta\end{subarray}}^{L}\sum_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{N}T_{j\beta}^{i\alpha}p_{i\alpha}(t)}_{\text{switch and
jump}}$

where $p_{j\beta}(t)$ indicates the probability of finding a random walker in
node $j$ of layer $\beta$ at time $t$, and the contributions of jumps and
switches are made explicit. Note that using the tensorial formalism, this
expression would be compressed into the elegant master equation

$\displaystyle p_{j\beta}(t+1)=p^{i\alpha}(t)T_{i\alpha}^{j\beta},$ (62)

with the continuous-time version described by the forward equation

$\dot{p}_{j\beta}(t)=-\tilde{L}_{j\beta}^{i\alpha}p_{i\alpha}(t),$ (63)

where
$\tilde{L}_{j\beta}^{i\alpha}=\delta_{j\beta}^{i\alpha}-T_{j\beta}^{i\alpha}$
is the random walk normalized Laplacian tensor.

It is important to remark that a random walk dynamics is affected by both
structure and transition mechanisms: i.e., by fixing the topology one can
always define distinct transition rules encoding different stochastic
movements between nodes and exploration strategies, which we have previously
named flavors.

As mentioned before, in PageRank random walks (PRRW) a teleportation (or non-
local jumping) parameter tunes the ability of the random walker to escape from
absorbing states and to reach nodes which are not in the neighborhood of the
current node [brin1998anatomy, de2015ranking]. The random walk introduced in
the previous section for monoplex networks is known as a classical random walk
(CRW) [de2013mathematical, de2014navigability]. There are several walks that
can be defined on multilayer networks, which generalize their monoplex
counterparts or are specific to multilayer systems. For instance, one can
define multilayer diffusive random walks (DRW) – generalizing the one defined
in [de2014navigability] –, maximal-entropy random walks (MERW)
[de2014navigability] – generalizing monoplex walks which localize around
topological defects [burda2009localization], and a physical random walk with
relaxation (PrRW) [de2015identifying]. Transition probabilities for these
processes are reported in Tab. 1.

| CRW | PRRW | DRW | MERW | PrRW
---|---|---|---|---|---
$T_{j\beta}^{j\beta}$ | $\frac{M_{j\beta}^{j\beta}}{\sigma_{j}(\beta)}$ | $r\frac{M_{j\beta}^{j\beta}}{\sigma_{j}(\beta)}+\frac{1-r}{NL}$ | $\frac{s_{\max}+M_{j\beta}^{j\beta}-\sigma_{j}(\beta)}{s_{\max}}$ | $\frac{M_{j\beta}^{j\beta}}{\lambda_{\max}}$ | $(1-r)\frac{M_{j\beta}^{j\beta}}{s_{j}(\beta)}+r\frac{M_{j\beta}^{j\beta}}{s_{j}}$
$T_{j\beta}^{j\alpha}$ | $\frac{M_{j\beta}^{j\alpha}}{\sigma_{j}(\alpha)}$ | $r\frac{M_{j\beta}^{j\alpha}}{\sigma_{j}(\alpha)}+\frac{1-r}{NL}$ | $\frac{M_{j\beta}^{j\alpha}}{s_{\max}}$ | $\frac{M_{j\beta}^{j\alpha}}{\lambda_{\max}}\frac{V_{j\beta}}{V_{j\alpha}}$ | $r\frac{M_{j\beta}^{j\beta}}{s_{j}}$
$T_{j\beta}^{i\beta}$ | $\frac{M_{j\beta}^{i\beta}}{\sigma_{i}(\beta)}$ | $r\frac{M_{j\beta}^{i\beta}}{\sigma_{i}(\beta)}+\frac{1-r}{NL}$ | $\frac{M_{j\beta}^{i\beta}}{s_{\max}}$ | $\frac{M_{j\beta}^{i\beta}}{\lambda_{\max}}\frac{V_{j\beta}}{V_{i\beta}}$ | $(1-r)\frac{M_{j\beta}^{i\beta}}{s_{i}(\beta)}+r\frac{M_{j\beta}^{i\beta}}{s_{i}}$
$T_{j\beta}^{i\alpha}$ | $\frac{M_{j\beta}^{i\alpha}}{\sigma_{i}(\alpha)}$ | $r\frac{M_{j\beta}^{i\alpha}}{\sigma_{i}(\alpha)}+\frac{1-r}{NL}$ | $\frac{M_{j\beta}^{i\alpha}}{s_{\max}}$ | $\frac{M_{j\beta}^{i\alpha}}{\lambda_{\max}}\frac{V_{j\beta}}{V_{i\alpha}}$ | $r\frac{M_{j\beta}^{i\beta}}{s_{i}}$
Table 1: Entries of the transition tensor for distinct random walks mentioned
in the text, namely classical (CRW), PageRank (PRRW), diffusive (DRW),
maximal-entropy (MERW), and physical with relaxation (PrRW) random walks. To
keep the notation simple, we use $\sigma_{j}(\beta)=s_{j}(\beta)+S_{j}(\beta)$
and $s_{\max}=\max\limits_{i,\alpha}\left\\{\sigma_{i}(\alpha)\right\\}$ for
the outgoing strengths (see the text for details). Note that the teleportation
parameter is usually denoted by $\alpha$: to avoid confusion, here we indicate
it by $r$. Finally, $\lambda_{\max}$ is the largest eigenvalue of the
multilayer adjacency tensor (see Sec. III.3.3). It is worth remarking that for
CRW, DRW and MERW, these transition rules generalize to any multilayer system
the ones introduced in [de2014navigability] for multiplex networks; PrRW is
defined as in [de2015identifying]; PRRW generalizes the walk introduced in
[de2015ranking]. Reprinted table with permission from
[bertagnolli2020diffusion]. Copyright (2021) by the American Physical Society.

We conclude this section with a short discussion about the physical random
walk, introduced in [de2014navigability] to describe those dynamics where the
state nodes have a “common memory” such that information reaching one state
node is _instantaneously_ diffused to all state nodes of the corresponding
physical node. This mechanism encodes a variety of processes in the real
world. For instance, in a network of digital interactions among individuals,
such as social media platforms, we can get a rumor in a particular platform
through our intra-layer connections, but at the same time that rumor will be
known to our alter egos in the other social media platforms. In this scenario,
inter-layer connectivity does not carry physical meaning and, consequently, is
ignored. The physical random walk with relaxation (PrRW) [de2015identifying]
can be seen as a variant of PRW, where the knowledge of inter-layer
connectivity is not required and is therefore discarded: from Tab. 1 it is
easy to identify that its transition probabilities balance intra- and inter-
links by means of probability weights $1-r$ and $r$, respectively.

It is worth remarking here that several other types of walks can be defined,
including ones describing quantum processes. See [biamonte2019complex] for a
review of this topic. From a mathematical perspective, one can write a master
equation similar in shape to the one of a classical random walk:

$\displaystyle|\dot{\psi}\rangle=-i\mathbf{L}_{Q}|\psi\rangle,$ (64)

where $i$ is the imaginary unit and
$\mathbf{L}_{Q}=\mathbf{D}^{-1/2}\mathbf{L}\mathbf{D}^{-1/2}$, $\mathbf{L}$
being the combinatorial Laplacian matrix. The evolution equation for this type
of continuous-time quantum walk is time-reversibile and, unlike classical
random walks, quantum ones: i) do not admit a stationary distribution and ii)
are deterministic in their dynamic nature and stochastic with respect to
measurements [bottcher2020walks]. In the future, it will be of interest to
study the emergence of potentially unseen physical properties related to
quantum walks on the top of multilayer networks.

##### Continuous-time diffusion

Let us assume we have some quantity, e.g., water, free to flow in a network of
nodes through pipes, which can be represented by links between pairs of nodes.
Given an initial distribution, what is the level of this quantity in each node
at a given time $t$? This problem can be casted into continuous-time diffusion
on a network.

In the following, let us indicate by $x_{i}(t)$ the state vector, in tensorial
formalism, carrying information about the quantity in each node of a monoplex
network at time $t$. Let $x_{i}(0)$ be the initial state vector: its evolution
over time can be modeled by the diffusion equation

$\displaystyle\dot{x}_{j}(t)=\mathcal{D}\left[W^{i}_{j}x_{i}(t)-W^{i}_{k}u_{i}e^{k}(j)x_{j}(t)\right]\,,$
(65)

where $\mathcal{D}$ is a diffusion constant, $u_{i}$ is the vector of ones and