$e^{k}(j)$ is a canonical vector (see Sec. II.1). Since $s_{k}=W^{i}_{k}u_{i}$
is the outgoing strength and
$s_{k}e^{k}(j)x_{j}(t)=s_{k}e^{k}(j)\delta^{i}_{j}x_{i}(t)$, the diffusion
equation can be written in more compact form as:

$\displaystyle\frac{dx_{j}(t)}{dt}=-\mathcal{D}L^{i}_{j}x_{i}(t)\,,$ (66)

where $L^{i}_{j}=W^{l}_{k}u_{l}e^{k}(j)\delta^{i}_{j}-W^{i}_{j}$ is the
combinatorial Laplacian tensor [chung1997]. It is worth remarking here that
this tensor differs from the normalized random walk Laplacian introduced in
the previous section, although in some cases – as for a classical random walk
– they are related by the simple relationship
$\tilde{L}_{j}^{i}=(D^{-1})^{i}_{k}L^{k}_{j}$. The solution of Eq. (66) is
given by $x_{j}(t)=x_{i}(0)e^{-\mathcal{D}L^{i}_{j}t}$, similarly to what we
have seen for continuous-time random walks. However, at variance with random
walks, it can be shown that in the stationary regime, the solution takes the
form $x_{j}(\infty)\propto u_{j}$, i.e., one will find exactly the same
fraction of the quantity in each node, uniformly distributed.

We might wonder how fast diffusion happens. Since the Laplacian matrix is
semi-positive definite, it is possible to show that it can be decomposed as

$\displaystyle L^{i}_{j}=Q^{i}_{h}\Lambda^{h}_{k}(Q^{-1})^{k}_{j},$ (67)

where $Q^{i}_{h}$ is a matrix whose columns are the eigenvectors of the
Laplacian and $\Lambda^{h}_{k}$ is a diagonal matrix whose entries are the
Laplacian’s eigenvalues. This algebraic feature allows us to prove that
$e^{-\mathcal{D}L^{i}_{j}t}=Q^{i}_{h}(e^{-\mathcal{D}\Lambda
t})^{h}_{k}(Q^{-1})^{k}_{j}$, highlighting that the exponential decay of
$x_{j}(t)$ is dominated by the smallest positive eigenvalue, which usually is
indicated by $\Lambda_{2}$. Therefore, the diffusion temporal scale is given
by $\tau\approx 1/\Lambda_{2}$.

Figure 35:  Continuous-time diffusion on a duplex network, i.e., a multiplex
consisting of two layers, obtained from two Erdős–Rényi networks with
independently wiring probabilities $p_{1},p_{2}\in[0,1]$, used to connect two
pairs of nodes within the same layer. Diffusion speed is quantified by the
second smallest eigenvalue of the Laplacian tensor, $\Lambda_{2}$: depending
on the wiring probabilities, diffusion in the duplex can be faster (left-hand
side panel) or slower (right-hand side panel) than in each layer separately.
The condition for enhanced diffusion to happen is given by
$\Lambda_{2}^{\text{multiplex}}\geq\max\\{\Lambda_{2}^{\text{layer
1}},\Lambda_{2}^{\text{layer 2}}\\}$: a sharp change in the behavior of the
characteristic temporal scale can be observed for varying weight of the inter-
layer connections (left and right panels), with a clear transition between two
distinct regimes above a certain critical value of inter-layer coupling. The
middle panel shows when the enhanced diffusion condition holds (encoded by
colors), while varying the wiring probabilities. Figure from [de2016physics].

In the case of interconnected multilayer networks, the diffusion equation has
been generalized by means of the supra-adjacency matrix [gomez2013diffusion]
and the tensorial formulation [de2013mathematical]. In this new setup, a
quantity can diffuse through inter-layer connections as well. If we indicate
by $X_{i\alpha}(t)$ the rank-2 state tensor at time $t$, then the multilayer
diffusion equation can be written as

$\displaystyle\frac{dX_{j\beta}(t)}{dt}=M^{i\alpha}_{j\beta}X_{i\alpha}(t)-M^{i\alpha}_{k\gamma}U_{i\alpha}E^{k\gamma}(i\beta)X_{i\beta}(t)\,,$
(68)

where $U_{i\alpha}=u_{i}u_{\alpha}$ and
$E^{k\gamma}(i\beta)=e^{k}(i)e^{\gamma}(\beta)$. If we define the multilayer
combinatorial Laplacian tensor as

$L^{i\alpha}_{j\beta}=M^{l\epsilon}_{k\gamma}U_{l\epsilon}E^{k\gamma}(j\beta)\delta^{i\alpha}_{j\beta}-M^{i\alpha}_{j\beta},$
(69)

the diffusion equation can be written more compactly as

$\displaystyle\frac{dX_{j\beta}(t)}{dt}=-L^{i\alpha}_{j\beta}X_{i\alpha}(t)\,,$
(70)

whose solution is given by
$X_{j\beta}(t)=X_{i\alpha}(0)e^{-L^{i\alpha}_{j\beta}t}$, a clear
generalization of the result obtained in the case of single-layer networks.
Also in this case, the second smallest eigenvalue $\Lambda_{2}$ – calculated
from the supra-adjacency matrix representation – governs the speed of
diffusion [de2013mathematical, sole2013spectral, gomez2013diffusion], leading
to interesting phenomena (see Fig. 35 for details).

##### Topological transition with diffusive processes

Figure 36:  Algebraic phase transition. In (a) and (e), the intra-layer
connectivity of a duplex (i.e., a multiplex consisting of 2 layers) is shown,
where the color of each node refers to the value of the corresponding
component in the Fiedler vector, i.e., the eigenvector associated with the
second smallest eigenvalue $\lambda_{2}$. In (a) $p=p^{*}=0.602$ and in (e)
$p=0.603$, In the middle column is displayed the eigenvalue $\lambda_{2}$ (b),
the inner product $\langle v_{A}|v_{B}\rangle$ (c), and, in (d), inner
products of $|v_{A}\rangle$ and $|v_{B}\rangle$ with the unit vector, i.e.,
the sum of their elements (see the text for details). Figure reproduced from
[radicchi2013abrupt].

An important question concerns the emergence of such effects in spreading
processes, but not only. One can adapt the dynamical rules of models defined
in monolayers to encompass the fact that the topology is more complex and
this, as we will see, can lead to substantially different behaviors than those
observed in isolated networks. Another point of view, though, is to ask why
and when a certain model could effectively behave as in an isolated network
even though it runs in an a layered structure, or in the other way around, why
and when the multilayer dimension dominates, thus disregarding any role of the
individual networks of the layers. In some special cases, such as
interconnected multiplex networks with identical coupling, it is possible to
show the existence of two distinct regimes as a function of the inter-layer
coupling strength [radicchi2013abrupt], highlighting how the multilayer
structure can influence the outcome of several physical processes. By
considering a duplex network – i.e., a multiplex with 2 layers, $A$ and $B$ –
and by using spectral properties of multilayer systems, it was found an abrupt
transition between the two aforementioned regimes as a function of $p$, the
coupling strength. The two regimes are inferred by analyzing the behavior of
eigenvectors and eigenvalues of the supra-Laplacian matrix and are clearly
distinguishable, separated by a critical point $p^{*}$. For $p\leq p^{*}$, the
second smallest eigenvalue111111Note that here we are using $\lambda_{2}$
instead of $\Lambda_{2}$, as in the previous section, to keep the same
notation of the original paper by [radicchi2013abrupt]. $\lambda_{2}$, which
is associated with a plethora of network properties, is independent of the
structure of the layers and hence the dynamical processes can be studied
separately, while in the regime $p>p^{*}$, $\lambda_{2}$ tends to a value
independent of $p$, i.e., depending only on the details of the intra-layer
networks.

The eigenvector $|v\rangle$ associated with $\lambda_{2}(p)$ can be split into
$|v_{A}\rangle$ and $|v_{B}\rangle$, which correspond to the elements of
$|v\rangle$ associated with nodes of networks $A$ and $B$, respectively. It
can be proved [radicchi2013abrupt] that for, $p\leq p^{*}$,
$|v_{A}\rangle=-|v_{B}\rangle=0$ holds, where
$|v_{A}\rangle=\pm\frac{1}{\sqrt{2N}}|1\rangle$, while for $p>p^{*}$ we have
$\langle v_{A}|1\rangle=\langle v_{B}|1\rangle=0$. Physically, this means that
in the subcritical regime the layers are structurally independent whereas in
the supercritical regime the interlayer connection dominates, imposing the
same sign in the eigenvector for nodes across networks and alternating the
sign for nodes in the same layer. Therefore, the algebraic phase transition
can be visualized in several ways. Panels in Fig. 36 show the behavior of the
eigenvectors in the subcritical and supercritical regions. Here,
$\lambda_{2}(p)$ displays a singular point $p^{*}$ at which the first
derivative is not continuous, a sign of an abrupt transition. For $p\leq
p^{*}$, $\lambda_{2}$ grows as $2p$, while in the other regime it tends to the
value that would take for a weighted superposition of the two layers $A$ and
$B$, whose Laplacian is $(1/2)(\mathcal{L}_{A}+\mathcal{L}_{B})$. This
abruptness can be observed more directly in the middle and lower panels, which
show the behavior of $\langle v_{A}|v_{B}\rangle$, $\langle v_{A}|1\rangle$
and $\langle v_{B}|1\rangle$ as a function of $p$.

###  IV.2 Synchronization processes

Synchronization is an emergent phenomenon of a population of dynamically
interacting units that, usually with a second-order phase transition
[nelson1977recent, stanley1999scaling, dorogovtsev2008critical], start
operating in a collective, coherent way. Synchronization phenomena may be
found in biology, sociology, and ecology and include birds flocking, fireflies
flashing, people singing, and neurons spiking, just to mention a few examples.
Once a fully synchronized state is reached, the (linear) stability of such a
state is tested by studying the effects of a small perturbation of the system
state. This approach was introduced in [Pecora1998a] for simple network
configurations and has been widely used and extended to complex network
topologies. Here we briefly introduce the framework of Master Stability
Function for a network of oscillators and then extend the formalism to
multilayer networks. For a more exhaustive description see [Boccaletti2018,
Arenas2008].

It is worth remarking that, in the following, we will use the more traditional
vector notation where $\mathbf{x}(t)$ indicates the state of the system at
time $t$. In some cases, we will also use the Kronecker product operator
$\otimes$, as is usual in equations governing synchronization dynamics. Our
choice is to help the reader link these concepts to the original studies which
introduced them. Nevertheless, the whole section could consider the tensorial
formalism, where the system state is indicated by the rank-1 tensor
$x_{\ell}(t)$ and where products such as $\mathbf{A}\otimes\mathbf{B}$ are
indicated as $A^{\alpha}_{\beta}B^{\gamma}_{\delta}$.

Let us consider a network of $N$ identical oscillators in an $m$-dimensional
space, where, in the absence of any interaction, the dynamics of each node $i$
is described by:

$\dot{\textbf{x}_{i}}=\textbf{F}(\textbf{x}_{i}),\qquad
i=1,2,...,N;\textbf{x}_{i}\in\mathbb{R}^{m}.$ (71)

We introduce an interaction between oscillators due to the fact that they are
coupled in an unweighted network specified by the adjacency matrix
$\textbf{A}=\\{A_{ij}\\}$ and we define the output function
$\textbf{H}(\textbf{x})$ as the function that governs the interaction between
nodes. We also assume that the coupling between the oscillators is diffusive,
that is the effect that node $j$ has on node $i$ is proportional to the
difference between $\textbf{H}(\textbf{x}_{j})$ and
$\textbf{H}(\textbf{x}_{i})$. Then, the evolution of the state of node $i$ is
given by:

$\dot{\textbf{x}_{i}}=\textbf{F}(\textbf{x}_{i})+\sigma\sum_{i=1}^{N}A_{ij}[\textbf{H}(\textbf{x}_{j})-\textbf{H}(\textbf{x}_{i})]=\textbf{F}(\textbf{x}_{i})-\sigma\sum_{j=1}^{N}L_{ij}\textbf{H}(\textbf{x}_{j}),$
(72)

where L is the Laplacian matrix and $\sigma$ is the coupling strength .

##### Stability of synchronized states

The MSF approach to test the stability of a fully synchronized state,
described in Appendix A, was extended to multilayer complex systems
[DelGenio2016]. In this work, it is considered a network with $M$ different
layers, each layer representing a different kind of interaction between nodes.
Eq. (72) is thus extended to describe the dynamics of the whole system:

$\dot{\textbf{x}_{i}}=\textbf{F}(\textbf{x}_{i})-\sum_{\alpha=1}^{M}\sigma_{\alpha}\sum_{j=1}^{N}L^{(\alpha)}_{ij}\textbf{H}^{(\alpha)}(\textbf{x}_{j}),$
(73)

where $\alpha$ is the index accounting for layers. We obtain the $M$-parameter
equation describing the time evolution of perturbation error:

$\dot{\boldsymbol{\xi}_{i}}=\left[J\textbf{F}(\textbf{s})-\sum_{\alpha=1}^{M}\sigma_{\alpha}\lambda^{(\alpha)}_{i}J\textbf{H}_{\alpha}(\textbf{s})\right]\boldsymbol{\xi}_{i},$
(74)

where $\boldsymbol{\xi}_{i}$ is the eigenmode associated with the eigenvalue
$\lambda_{i}$ of L. As in the case of a dynamics evolving on top of a single
layer, in a multilayer system the stability of the synchronized state is
completely specified by the sign of the maximum conditional Lyapunov exponent
$\Lambda_{max}$. In particular, it is also found that stability of the
complete synchronization state may be reached even if each layer, taken
individually, is unstable: a very interesting feature for practical
application.

It is worth noting that, together with complete synchronization, a network may
exhibit other forms of synchronization where clusters of nodes have a
synchronized dynamics but different clusters evolve on distinct time
evolutions. This type of synchronization is called clustered synchronization
(CS) and it has been well studied in terms of cluster formation, stability and
role of network symmetries [Nicosia2013, Pecora2014, Sorrentino2016]. CS has
been also studied in multiplex [Jalan2016] and, more recently, in multilayer
[DellaRossa2020] networks, and cluster stability has been tested as a function
of intra- and inter-layer symmetries. To describe the evolution of the
perturbation error for complex synchronization patterns, such as in cluster
synchronization, it is useful to rewrite Eq. (105) with a more compact
formalism. To this end, we write the state vector as a vector of vectors,
$\textbf{X}=(\textbf{x}_{1};\textbf{x}_{2};...;\textbf{x}_{N})$, and the
variational equation assumes the form:

$\delta\dot{\textbf{X}}=\left[\textbf{I}_{N}\otimes
J\textbf{F}(\textbf{s})-\sigma\textbf{L}\otimes
J\textbf{H}(\textbf{s})\right]\delta\textbf{X},$ (75)

where $\textbf{I}_{N}$ is the identity matrix and $\otimes$ is the Kronecker
product. Eq. (75) can be decoupled into N independent equations by
diagonalizing L. However, to deal with cluster synchronization and multilayer
interactions we have to further generalize Eq. (75). The variational equations
for complex synchronization patterns on generalized networks has the form
[zhang2020]:

$\delta\dot{\textbf{X}}=\left[\sum_{l=1}^{L}\textbf{D}^{(l)}\otimes
J\textbf{F}(\textbf{s}^{l})-\sum_{l=1}^{L}\sum_{\alpha=1}^{M}\sigma_{\alpha}\textbf{L}^{(\alpha)}\textbf{D}^{(l)}\otimes
J\textbf{H}^{(\alpha)}(\textbf{s}^{l})\right]\delta\textbf{X}$ (76)

where the identity matrix has been replaced by the diagonal matrix
$\textbf{D}^{(l)}$, whose generic element $D_{ii}^{(l)}=1$ if node $i$ belongs
to $l$th dynamical cluster and $D_{ii}^{(l)}=0$ otherwise.

In a recent paper [zhang2020], it was established that to optimally decouple
Eq. (76), the matrices encoding the synchronization pattern and the
interaction pattern, that is $\\{\textbf{D}^{(l)}\\}$ and
$\\{\textbf{L}^{(\alpha)}\\}$, should be simultaneously diagonalized. In this
work, an algorithm was also developed to find the finest simultaneous block
diagonalization.

##### Synchronization in a network of phase oscillators

One of the first approaches to describe phase synchronization in an ensemble
of oscillators on multiplex networks was done in 2015 [Gambuzza2015], to
investigate the synchronization of indirectly coupled units using a system
composed by two layers, where the top layer was made of disconnected
oscillators and the bottom one, modeling the medium, consisted of oscillators
coupled according to a given topology and with a characteristic natural
frequency. Each node of the multiplex was modeled as a Stuart-Landau (SL)
oscillator, that is an oscillator with amplitude as well as phase dynamics.
Therefore, the Kuramoto model (KM, see Appendix B) can be retrieved as a
limiting case when the amplitude dynamics vanishes. The Kuramoto order
parameter can be generalized to the multiplex framework as:

$r^{\alpha\beta}_{ij}=\big{|}\big{\langle}e^{i[\theta^{\alpha}_{i}(t)-\theta^{\beta}_{j}(t)]}\big{\rangle}_{t}\big{|},$
(77)

while intra- and inter-layer coherence, respectively, can be also defined by

$r^{\alpha}=\frac{1}{N(N-1)}\sum^{N}_{i,j=1}r^{\alpha\alpha}_{ij},$ (78)

and

$r^{\alpha\beta}=\frac{1}{N}\sum^{N}_{j=1}r^{\alpha\beta}_{jj}.$ (79)

By studying a population of $N$ disconnected oscillators, indirectly coupled
through an inhomogeneous medium, authors have shown the onset of intra-layer
synchronization without inter-layer coherence, i.e. a state in which the nodes
of a layer are synchronized between them without being synchronized with those
of the other layer (see Fig. 37, left panel). Synchronization of units that
are not connected requires the presence of an amplitude dynamics as the regime
of intra-layer synchronization is not observed in purely phase oscillators,
such as those in the KM.

Figure 37: Left: Two-layer multiplex network with one-to-one coupling between
the layers. In the top layer ($\alpha$) the nodes only interact with those in
the bottom one whereas in the bottom layer ($\beta$) the nodes also interact
with other members of the same layer. Middle: Kuramoto order parameters
($r^{\alpha}$) and ($r^{\alpha\beta}$) vs. coupling coefficients
$\lambda=\lambda_{\alpha\beta}=\lambda_{\beta}/5$. Continuous lines refer to a
multilayer network of Stuart-Landau oscillators with $a=1$, whereas the dashed
ones to purely phase oscillators ($a\xrightarrow{}\infty$). Right:
Synchronization transitions in two-layer networks with a fraction of the nodes
adaptively controlled by a local order parameter $f=1$. Squares and circles
(triangles and stars) refer to the values of $R1$ ($R2$), and the insets show
the corresponding dependence of the width of hysteretic loop $d$ on $f$. Left
and middle panels readapted from [Gambuzza2015], right panel from [Zhang2015]

As previously mentioned, not only were continuous second-order transitions
observed in an ensemble of networked phase oscillators, but examples of an
abrupt first-order transition, named explosive synchronization (ES), were also
observed [Gomez-Gardenes2011, Leyva2013]. It was first observed in a network
of oscillators presenting a positive correlations between natural frequencies
and the degree of the nodes. Moreover, ES was also studied in systems where a
local order parameter for the $i$-th oscillator is defined [Zhang2015]:

$r_{i}(t)e^{i\Phi(t)}=\frac{1}{k_{i}}\sum_{j=1}^{k_{i}}\sin(\theta_{j}),$ (80)

and where the phase dynamics is expressed as:

$\dot{\theta_{i}}=\omega_{i}+\sigma\alpha_{i}\sum_{j=1}^{N}A_{ij}\sin(\theta_{j}-\theta_{i}).$
(81)

The overall amount of phase coherence in the network is measured by means of
the global order parameter $R$:

$Re^{i\Psi}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}},$ (82)

where $0\leq R\leq 1$ and $\Psi$ denotes the average phase.

Explosive synchronization onset has been reported in a system of two
interdependent networks (see Fig. 37, right panel), with the same size and
where nodes on the two layers are coupled in a one-to-one correspondence, so
that a group of oscillators in the first layer is controlled by the local
order parameters of the corresponding nodes on the second layer, and vice-
versa [Zhang2015]. Nevertheless, it was shown that ES is a property of a
generic multilayer network as long as some microscopic suppressive rule can
prevent the formation of the giant synchronization cluster which characterizes
second-order transitions.

###  IV.3 Game dynamics: cooperation processes

Human cooperation may be intended as a collective behavior that emerges as the
result of the interactions among individuals. In the past few years
cooperation has been studied in social sciences with methods of statistical
physics [Perc2017], in particular Monte Carlo methods and the theory of
collective behavior of interacting particles near phase-transition points.
That approach has proven very valuable for understanding cooperation and its
spatio-temporal dynamics. The mathematical framework used to study human
cooperation is usually evolutionary game theory, which quantitatively
describes social interactions using example games and formalizes the concept
of the social dilemma, intended as the conflictual choice that an individual
has between doing what is best for society or doing what is best for
themselves.

Figure 38: a) Mean final cooperation as function of the temptation to defect
$T$, for several multiplex networks with different numbers of layers $M$. b)
Multiplex with $M=2$ layers, where cooperation is studied as a function of
synergy factor $r$ and edge overlap $\omega$. Dashed lines separate regions
with either full defection $c=0$ (IA) or full cooperation $c=1$ (IB), from the
region of continuously evolving coexistence of cooperators and defectors
$0<c<1$ (II). Panel a) readapted from [gomez2012evolution] and panel b)
readapted from [battiston2017determinants]. c) Mean final cooperation (color
coded) for the Prisoner’s Dilemma as a function $M$ and the strength of degree
correlations $\nu$. d) Mean final cooperation as a function of the initial
density of cooperators $c_{0}$ and the strength of correlations $\nu$. Panels
c) and d) readapted from [Kleineberg2018].

Evolution of cooperation strongly depends on the population interaction
network: individuals who have the same strategy are more likely to interact
[Nowak1992, Nowak2010], effectively creating resilient cooperative clusters in
a structured population, a phenomenon named network reciprocity. When the
edges that determine the interaction among individuals were fixed on a time
scale, it was demonstrated [Santos2006] that in heterogeneous populations –
modeled by networks with degree distributions exhibiting a power-law behavior
– the sustainability of cooperation is simpler to achieve than in
homogeneously structured populations. Here we address the problem of how human
cooperation emerges in multiplex networks, with different interaction layers
that can account for different kinds of social ties an individual may be
involved in.

In particular, we consider a Prisoner’s Dilemma game implemented in a set of
$M$ interdependent networks, characterized by an average degree $\langle
k\rangle$, each of them containing the same number $N$ of nodes. Each
individual is represented by one node in each of the layers and we define a
set of adjacency matrices $\\{A^{\alpha}\\}$ so that $A_{ij}^{\alpha}=1$ when
two nodes are connected in layer $\alpha$ and $A_{ij}^{\alpha}=0$ otherwise.

At each time step, for each of the $k_{i}^{\alpha}$ games played, each
individual $i$ facing a cooperator collects a payoff $\pi_{ij}=R$ or
$\pi_{ij}=T$ when playing as a cooperator or as a defector, respectively.
Conversely, if $i$ faces a defector, $i$ collects a payoff $\pi_{ij}=S$ or
$\pi_{ij}=P$ playing as a cooperator or as a defector, respectively. For game
parameters it holds that $S<P<R<T$. The aggregated payoff of node $i$ is given
by the sum of the payoffs $\pi_{i}$ over all layers. Furthermore, after each
round of the game, individuals update their strategies with a rule that can
use the global knowledge about the benefits of the neighbors or be random
[Gomez-Gardenes2011b].

Finally, it worth noting that cooperation is also studied in public good
games, where the Prisoner’s Dilemma is played in overlapping groups of
individuals: cooperators contribute with a “cost” $d$ to the public good,
while defectors do not contribute and the total amount in the common pool of
each group is multiplied by a factor $r$ and distributed equally among all
members of the group [Santos2008a].

Previous research [gomez2012evolution] demonstrated that the resilience of
cooperative behavior can be enhanced by the multiplex structure through the
simultaneous formation of correlated clusters of cooperators across different
layers, i.e. through multiplex network reciprocity. However, it was also
proven [battiston2017determinants] that, to gain benefits from multiplex
structure, a high topological overlap is needed or, in other words,
individuals have to be similarly linked across different layers (see Fig. 38,
panels a) and b)). Other studies [Kleineberg2018] investigated the role of
degree correlation $\nu$ between nodes in different layers for the emergence
of cooperation. They found that in the absence of degree correlations,
increasing the number of layers only leads to mild changes. However, if degree
correlations are present we observe a mean final cooperation of $c=0.5$, and
this value is nearly independent of the game payoff parameters. This mechanism
is called _topological enslavement_ and can be understood by considering that,
if degree correlations are strong, hubs dominate the game dynamics, since they
have the potential to earn higher payoffs (because they play more games) and
they are more likely to be selected by other nodes as imitation candidates.
Furthermore, topological enslavement implies that the outcome of the evolution
of the system is determined by the initial conditions (see Fig. 38, panels c)
and d)).

Finally, as anticipated in Sec. III.5, layer-layer correlations might have a
deep impact on the dynamics on the top of multilayer systems. This is the case
for game dynamics, where cooperation might be hampered, rather than enhanced,
by specific correlation patterns combining assortative and disassortative
degree mixing across layers [wang_pre14a, duh_njp19] (see Fig. 39 for
details).

Figure 39: Fraction of cooperators as a function of the payoff of a cooperator
when playing with a defector ($S$) and the payoff of a defector when facing a
cooperator ($T$). The interaction layer is subject to disassortative mixing
(assortative coefficient $A_{I}<0$) while the updating network is subject to
assortative mixing (disassortative coefficient $A_{U}>0$). Results show the
disruptive effect of symmetry breaking for the evolution of cooperation.
Parameter values are $A_{I}=-0.1$ and $A_{U}=0.1$ (Left) and $A_{I}=-0.2$ and
$A_{U}=0.2$ (Right). Figure from [wang_pre14a].

###  IV.4 Interdependent processes

The second class of dynamical processes on multilayer networks is the one of
“interdependent” or “coupled” dynamics, consisting of systems characterized by
different processes on each layer. The inter-layer interactions couple such
processes and are responsible for emerging phenomena which could not be
detected in a single-layer network framework (see Tab. 3).

As in the case of the structure, it is possible to introduce a unifying
framework in terms of a dynamical SNXI decomposition, outlined in Eq. (8), to
describe dynamics on multilayer networks. Let $x^{[l]}_{i\alpha}$ (where
$l\in{1,2,...,C}$) denote the $l$-th component of a $C-dimensional$ vector
$\mathbf{x}_{i\alpha}$ that represents the state of node $i$ in layer
$\alpha$. The most general (and possibly nonlinear) dynamics governing the
evolution of each state is given by the systems of equations

$\displaystyle\dot{\mathbf{x}}_{i\alpha}(t)$
$\displaystyle=F_{i\alpha}(\mathbf{X}(t))=\sum_{\beta=1}^{L}\sum_{j=1}^{N}f_{i\alpha}^{j\beta}(\mathbf{X}(t))$
(83)
$\displaystyle=\underbrace{\sum_{\beta=1}^{L}\sum_{j=1}^{N}f_{i\alpha}^{j\beta}(\mathbf{X}(t))\delta_{\alpha}^{\beta}\delta_{i}^{j}+\sum_{\beta=1}^{L}\sum_{j=1}^{N}f_{i\alpha}^{j\beta}(\mathbf{X}(t))\delta_{\alpha}^{\beta}\left(1-\delta_{i}^{j}\right)}_{\text{intra-
layer dynamics }}$
$\displaystyle+\underbrace{\sum_{\beta=1}^{L}\sum_{j=1}^{N}f_{i\alpha}^{j\beta}(\mathbf{X}(t))\left(1-\delta_{\alpha}^{\beta}\right)\delta_{i}^{j}+\sum_{\beta=1}^{L}\sum_{j=1}^{N}f_{i\alpha}^{j\beta}(\mathbf{X}(t))\left(1-\delta_{\alpha}^{\beta}\right)\left(1-\delta_{i}^{j}\right)}_{\text{inter-
layer dynamics }}$
$\displaystyle=\underbrace{f_{i\alpha}^{i\alpha}(\mathbf{X}(t))}_{\text{self-
interaction }}+\underbrace{\sum_{j\neq
i}f_{i\alpha}^{j\alpha}(\mathbf{X}(t))}_{\text{endogenous interaction
}}+\underbrace{\sum_{\beta\neq\alpha}\sum_{j\neq
i}f_{i\alpha}^{j\beta}(\mathbf{X}(t))}_{\text{exogenous interaction
}}+\underbrace{\sum_{\beta\neq\alpha}f_{i\alpha}^{i\beta}(\mathbf{X}(t))}_{\text{intertwining
}}$
$\displaystyle=\mathbb{S}_{i\alpha}(\mathbf{X}(t))+\mathbb{N}_{i\alpha}(\mathbf{X}(t))+\mathbb{X}_{i\alpha}(\mathbf{X}(t))+\mathbb{I}_{i\alpha}(\mathbf{X}(t))$

where
$\mathbf{X}(t)\equiv\left(\mathbf{x}_{11},\mathbf{x}_{21},\ldots,\mathbf{x}_{N1},\mathbf{x}_{12},\mathbf{x}_{22},\ldots,\mathbf{x}_{N2},\ldots,\mathbf{x}_{1L},\mathbf{x}_{2L},\ldots,\mathbf{x}_{NL}\right)$
and we did not use the tensorial formalism to make explicit the contributions
of each term in the governing equation. In fact, this equation could also be
compactly written as

$\displaystyle\dot{x}_{j\beta}(t)=F[x_{j\beta}(t)]=\mathbb{S}[x_{j\beta}(t)]+\mathbb{N}[x_{j\beta}(t)]+\mathbb{X}[x_{j\beta}(t)]+\mathbb{I}[x_{j\beta}(t)].$
(84)

We call this equation “dynamical $\mathbb{SNXI}$ decomposition”. Similarly to
the structural decomposition in Eq. (8), we have decoupled the different
contributions of intra-layer and inter-layer dynamics, allowing us to classify
different dynamical processes in terms of the corresponding dynamical
$\mathbb{SNXI}$ components. The peculiar behavior of the interdependent
processes is ascribed to the exogenous and intertwining components. Although
the most-studied examples come from mixed spreading processes, which are
crucial for understanding phenomena such as the spreading dynamics of two
concurrent diseases in two-layer multiplex networks [Sanz2014,
salehi2015spreading, Dickison2012epidemics, Cozzo2013social, DeArruda2017],
and the spread of diseases coupled with the spread of information or behavior
[wang2015coupled, funk2015, Granell2013, lima2015disease, Funk2009,
granell2014competing], other types of dynamical interdependence are attracting
a growing interest.

##### Coupling diffusion with synchronization

An illustrative and pedagogical example has been proposed by
[nicosia2017collective]. In this work, the authors examine the interdependent
dynamics of the two processes presented in the two previous sections, namely,
diffusion and synchronization. They propose a model that mimics the interplay
between the neural activity and energy transport in brain regions, from which
a rich collection of behaviors emerges. These processes evolve in the two
layers of a multilayer network and are related to each other by the
correspondence between layers, which in this case is realized through the
functional relation between the parameters governing the two processes and the
state variables. In particular, the dynamics of the entire system is governed
by the following equations:

$\left\\{\begin{aligned}
\dot{x}_{i}&=F_{\omega_{i}}\left(\mathbf{x},A^{[1]}\right)\\\
\dot{y}_{i}&=G_{\chi_{i}}\left(\mathbf{y},A^{[2]}\right)\end{aligned}\quad
i=1,2,\ldots,N\right.$ (85)

where $\mathbf{x}=\left\\{x_{1},x_{2},\ldots,x_{N}\right\\}\in\mathbb{R}^{N}$
and $\mathbf{y}=\left\\{y_{1},y_{2},\ldots,y_{N}\right\\}\in$ $\mathbb{R}^{N}$
denote the states of the two dynamical processes, while the topologies of the
two layers are encoded in the adjacency matrices
$A^{[1]}=\left\\{a_{ij}^{[1]}\right\\}$ and
$A^{[2]}=\left\\{a_{ij}^{[2]}\right\\},$ respectively, such that
$a_{ij}^{[1]}=1\left(a_{ij}^{[2]}=1\right)$ if a link exists between nodes $i$
and $j$ in the first (second) layer, and $a_{ij}^{[1]}=0$
$\left(a_{ij}^{[2]}=0\right)$ otherwise. The evolution of the system depends
on a set of parameters $\omega$ and $\chi$, which in turn depend on the state
of the nodes in the adjacent layer:

$\begin{array}[]{l}\dot{\omega}_{i}=f\left(\omega_{i},y_{i}\right)\\\
\dot{\chi}_{i}=g\left(\chi_{i},x_{i}\right)\end{array}\quad i=1,2,\ldots N$
(86)

The authors assign to functions $F_{\omega_{i}}$ and $G_{\chi_{i}}$ a Kuramoto
dynamic and a continuous-time random walk, respectively. Subsequently, they
assign the functions $f$ and $g$ in Eqs. (86), respectively, relating the
frequency $\omega_{i}$ of the oscillator $i$ at layer 1 to the state $y_{i}$
at layer 2 , and the bias property $\chi_{i}$ of the random walkers at layer 2
to the oscillator phase $x_{i}$ at layer $1$. The natural frequency
$\omega_{i}$ of the oscillator $i$ evolves relaxing to values proportional to
the fraction of random walkers at the replica node $i$ in the other layer.
Analogously, the random walks are biased toward (away from) strongly
synchronized nodes.

Note that the coupling between the two layers is tunable through two
parameters $\lambda$ and $\alpha$ that represent, respectively, the intensity
of the influence of the random walk on the oscillators and vice versa. This
setup completely defines an interdependent process on a multilayer network:
depending on the coupling strengths $\lambda$ and $\alpha$ the collective
behavior exhibits special dynamics (see Fig. 40).

Figure 40:  Coupling diffusion with synchronization dynamics. a) Distribution
$P(y_{i})$ of steady-state random walker fractions $y_{i}$ at layer $2$ for
$\alpha=1.0$, when the oscillators at layer $1$ are incoherent ($\lambda=0.1$,
top, blue) and synchronized ($\lambda=0.4$ bottom, red). b) Synchronization
phase diagram showing the level of synchronization as a function of coupling
$\lambda$ and bias exponent $\alpha$. The bistable region is colored in white.
Figure adapted from [nicosia2017collective].

The random walkers are homogeneously distributed in the incoherent state,
while in the synchronized state the distribution is heterogeneous. Conversely,
at certain values of the tuning parameters the system encounters an explosive
synchronization, or a bistability region characterized by a hysteretic
behavior. Similar phase transitions can also be observed in single-layer
networks under certain conditions, as reported in [acebron2005kuramoto,
martens2009exact, buendia2021broad]. Importantly, the multilayer network model
offers a parsimonious explanation of the emergence of these collective
phenomena, considering explicitly the intertwined nature of the dynamics.

##### Coupling epidemics spreading with awareness diffusion

Many different phenomena in nature can be described, in their essence, by the
results of constructive or destructive relations between two or numerous
parts. For instance, the mutualistic or competitive relation between dynamical
systems gives rise to a wealth of fascinating behaviors.

The dynamics occurring in one layer can have positive or negative feedbacks on
another: for example, human behavior (e.g., information awareness) can inhibit
the spread of disease; social mixing between classes and mobility may produce
abrupt changes in the critical properties of the epidemic onset; cooperation
emerges where the classical expectation was defection.

A generalization of these results can be found in [danziger2019dynamic], in
which the authors describe some universal features of interdependent systems
with coupled dynamics. There is an interesting relation between the collective
behavior emerging from percolation processes and from the ones arising in
interdependent systems. In such processes, abrupt transitions and critical
dynamics may arise in certain circumstances and, in dynamically coupled
systems, other interesting phenomena may be also observed, such as hysteresis
and multistability, with functionality of nodes in one layer influencing the
functionality of their replicas in the other layers. In [danziger2019dynamic]
the functionality is quantified by an order parameter, which plays a role in
the coupling strength between the nodes in different layers, with the order of
a node dynamics affecting the order of its neighbors. This fact is the
dynamical counterpart of interdependent percolation, where the functionality
of a node is related to its belonging to the mutual giant connected component
(see Sec. IV.5 for details). It turns out that the dynamic interdependence
increases the vulnerability of the system [danziger2016vulnerability], as in
the case of percolation processes.

Figure 41:  Two (left) reciprocally enhanced and (right) reciprocally
inhibited disease-spreading processes of susceptible–infected–susceptible
type. The colors in the figure represent the prevalence levels of the diseases
at a steady state of Monte Carlo simulations. Note the emergence of a curve of
critical points (at a “metacritical point”) in which the spreading in one
layer depends on the spreading in the other. Figure from [de2016physics].

Epidemic models are another important example of systems in which
interdependencies play a crucial role (see [wang2015coupled, de2016physics]
for a review) and, also in this case, hysteretic behavior with abrupt
transitions may occur. This means that, for instance, explosive pandemics with
no early-warning can suddenly appear and, in a similar implosive way, can
disappear. Unfortunately, the value of the infection rate which can eradicate
the disease must be much lower than the value that triggered it
[danziger2019dynamic]. Other interesting behaviors arise in the case of
competitive diseases, where mutual exclusion or endemic coexistence may
spontaneously occur.

Figure 41 shows two phase diagrams of disease incidence of reciprocally
enhanced (right) and inhibited (left) disease-spreading processes. In the
figure is highlighted the existence of a curve of critical points that
separate endemic and non-endemic phases of the disease. Moreover, the curve
that separates the endemic and non-endemic phases, is in turn divided into two
parts: one in which the critical properties of one spreading process are
independent of the other (straight dashed line), and one in which the critical
properties of one spreading process do depend on those of the other layer
(solid curve). The point at which this crossover occurs is called a
metacritical point.

Figure 42:  SIS-SIS interacting diseases model. Left: multiplex representation
of the diseases spreading. Each individual is present in both the layers and
can be infected by one (or both) of the diseases, as indicated in the central
panel. Right: possible transitions between the different states of the SIS
model. The variables represent the densities of individuals of each type
having degree equal to $k$ in the first layer and degree $l$ in the second.
Figure from [Sanz2014].

In fact, traditional epidemic models can describe in details the spreading of
a single disease in different realistic situations, whereas the multilayer
representation can extend the possibilities to spreading processes with
interacting diseases. In [Sanz2014] the authors propose a framework to
describe the spreading dynamics of two concurrent diseases (see Fig.42). Using
SIS-SIS and SIR-SIR models with appropriate interlayer couplings describing
the influence of one disease on the other, they derive the epidemic threshold
for the two interacting diseases, showing that the onset of a disease’s
outbreak is conditioned to the prevalence levels of the other disease.

Epidemic thresholds can be influenced not only by the presence of a second
disease, but also by the individual’s change of awareness about an ongoing
epidemic. In [Granell2013], a multilayer network couples the dynamics of
disease spreading in a social network and the diffusion of awareness among
actors. The two layers are coupled in a multiplex as illustrated in Fig. 43.
The process of spreading of awareness (UAU) is akin to an SIS process, where
in place of susceptible (S) there are unaware (U) and in place of infected (I)
there are aware (A) actors. The probability of being infected is influenced by
the state of awareness of the individuals. The authors found the existence of
a metacritical point in which the state of awareness of the individuals can
control the onset of an epidemic.

Figure 43: Multiplex representation of awareness-disease spreading. The
spreading of awareness occurs in the upper layer, whereas the spreading of the
diseases takes place in the lower layer. Figure from [Granell2013].

##### Coupling game dynamics with opinion dynamics

Figure 44: Polarization of opinions and strategies in a multiplex with 5,000
nodes, in the presence of angular correlations between the layers. Top:
comparison between the two layers; the color is the time average of the state
of each node. Bottom: evolution of the density of cooperators in the angular
bins used to compute the interlayer correlations. Numeric labels indicate the
clusters of nodes that adopt the same strategy. Figure from
[amato2017interplay].

The multilayer representation of the dynamics of complex systems can also be
used to shed light on real world social dilemmas. The influence of factors
acting on different layers might explain the emergence of particular patterns
of cooperation between social agents. In [amato2017interplay] the authors
present a model coupling evolutionary game dynamics and opinion dynamics,
regarded as two processes evolving on distinct layers of a multiplex network.

To model the game dynamics on the first layer the authors adopt a replicator,
in which the individuals (nodes) copy the strategy of one of their neighbors
with a probability that is so much higher the higher the neighbors’ payoff
compared to themselves. The possible states are ”cooperate” and ”defeat”. The
opinion dynamics on the second layer is modeled using the voter model, where
individuals (nodes) adopt the opinion of a randomly selected neighbor with a
certain probability. The opinions of the nodes can be ”cooperate” or ”defeat”,
as in the game layer. The authors assume that imitating a cooperative opinion
is more likely than imitating a defection. This can be interpreted as the
influence of media campaigns or broadcasting agents. Depending on the specific
parametrization of both the game and the opinion dynamics, this model gives
rise to fascinating dynamical behaviors in which equilibria of different types
exist for the game dynamics. The impact of social influence on the decision of
individuals is conveyed by the interlayer coupling, which is encoded in the
parameter $\gamma$, representing the tendency of individuals to act in
agreement with their proclamations: the nodes in one layer copy their own
state from the other layer with probability $\gamma$.

The main result is that this model is sufficient to explain the emergence of
cooperation in scenarios where the pure game dynamics predicts defection. This
is due both to the intertwined dynamics and to the multilayer structure
itself. In fact, the authors proved that the geometric correlation between
layers has a significant impact on the stability of the system. Importantly,
under certain conditions of correlation between layers, the system can reach a
polarized metastable state (see Fig. 44). This result can explain the observed
polarization in real world social systems.

Ultimately, the emergence of cooperation in unexpected conditions is due to
the interplay between the coupled dynamics of strategies and opinions, the
complex topologies of the networks upon which the dynamics exert, and the
structural relations between the layers. Missing one of these elements could
hinder the right interpretation of the complex behavior of such systems.

##### Coupling epidemics spreading with social integration

A multilayer perspective on epidemic spreading can be a valuable support to
design more educated strategies to reduce disease risk. A recent example comes
from [bosetti2020heterogeneity], which integrates social dynamics, human
mobility and epidemic spreading to assess the risk of measles outbreak in
Turkey. During the past decade, Turkey has received more than 3.5 million
refugees coming from Syria. The levels of immunization of the two populations
are considerably different. The outbreak risk is analysed through a multilayer
transmission model, which takes into account the different levels of
immunization in the two populations, along with the mobility pattern and the
level of social integration. The main result of the study is that, in the case
of heterogeneous immunization, high levels of social interaction can
drastically reduce the spatial spread and incidence of a disease. This
apparently counter-intuitive result is due to the fact that the high
immunization coverage of one population (Turkish citizens) can shield the
other (Syrian refugees) from getting exposed to the infection, as an effect of
herd immunity. The network structure is defined by dividing Turkey into
patches corresponding to administrative areas that represent the nodes. The
layers of the network are encoded by Turkish $(T)$ and Syrian refugee $(R)$
populations. On this network, social dynamics, mobility and epidemic spreading
happen simultaneously.

Figure 45: (A) Model scheme. Prefectures of Turkey are the nodes of a network
of geographic patches. Turkish and Syrian populations are encoded by two
colors and move between patches following the inferred mobility pathways. The
two populations encode two layers, where social dynamics and epidemic
spreading happen simultaneously. (B) Mobility of Syrian refugees (upper) and
Turkish citizens (lower) inferred from mobile phone data. Figure from
[bosetti2020heterogeneity].

The authors denotes with $c_{ki}^{(p)}(p\in T,R)$ the elements of a matrix
$\mathbf{C}^{(p)}$ encoding the number of people belonging to population
$p\in\\{T,R\\}$ traveling from patch $k$ to patch $i$, and with $\alpha$ the
fraction of Syrian contacts with Turkish citizens. The force of infection for
each population in the $i$-th patch depends on the contribution of all patches
in the country:

$\displaystyle\begin{array}[]{l}\lambda_{i}^{(T)}\left(\alpha,\mathbf{C}^{(T)},\mathbf{C}^{(R)}\right)=\beta_{T}\sum\limits_{k=1}^{L}\left[\underbrace{c_{ki}^{(T)}\frac{I_{k}^{(T)}}{N_{k}^{(T)}}}_{\text{Endogenous
}}+\underbrace{\alpha
c_{ki}^{(R)}\frac{I_{k}^{(R)}}{N_{k}^{(R)}}}_{\text{Exogenous }}\right]\\\
\lambda_{i}^{(R)}\left(\alpha,\mathbf{C}^{(T)},\mathbf{C}^{(R)}\right)=\beta_{R}\sum\limits_{k=1}^{L}\left[\underbrace{\alpha
c_{ki}^{(T)}\frac{I_{k}^{(T)}}{N_{k}^{(T)}}}_{\text{Exogenous
}}+\underbrace{c_{ki}^{(R)}\frac{I_{k}^{(R)}}{N_{k}^{(R)}}}_{\text{Endogenous
}}\right]\end{array}$

where $\beta_{p}=\beta/P_{i}^{(p)}(\alpha,c)$ is the transmission rate for
population $p$ and $P_{i}^{(p)}(\alpha,c)$ is an appropriate normalization
factor. From this equation we can easily recognise the structure of Eq. (83),
where the contributions to the force of infection come from an endogenous
term, accounting for the infectivity due to individuals from the same
population, and an exogenous term, accounting for the infectivity due to the
other population. The parameter $\alpha$ is the level of social mixing and can
be changed according to different social integration scenarios. This is the
parameter which plays the role of the coupling between the layers. An
illustrative representation of the model is shown in Fig. 45.

The epidemic transmission dynamics is eventually regulated by the following
SIR dynamical model

$\displaystyle\left\\{\begin{aligned}
\dot{S}_{i}^{1}&=-\lambda_{i}^{1}(\alpha,c)S_{i}^{1}\\\
\dot{S}_{i}^{2}&=-\lambda_{i}^{2}(\alpha,c)S_{i}^{2}\\\
\dot{I}_{i}^{1}&=\lambda_{i}^{1}(\alpha,c)S_{i}^{1}-\gamma I_{i}^{1}\\\
\dot{I}_{i}^{2}&=\lambda_{i}^{2}(\alpha,c)S_{i}^{2}-\gamma I_{i}^{2}\\\
\dot{R}_{i}^{1}&=\gamma I_{i}^{1}\\\ \dot{R}_{i}^{2}&=\gamma
I_{i}^{2}\end{aligned}\right.,$

whose analysis suggests that the incidence of the measles can be reduced up to
90% in the case of very high levels of integration [bosetti2020heterogeneity].

Despite the different contexts, the multilayer dynamics can have positive or
negative feedbacks, leading to interdependence between the corresponding
critical points of the dynamics. As a consequence, two different regimes
exist: i) one in which the critical properties of one process depend on those
of the other, and (ii) in which the critical properties are independent of the
other. The two regimes are separated by a metacritical point, where a
crossover occurs (see a recent review from [de2016physics] and Fig. 41). This
regime-shift is analogous to the one occurring in percolation processes, which
is presented in the next section.

###  IV.5 Percolation

In this section we switch our attention to percolation
[stauffer2018introduction], where the goal is to analyse how different
properties of the network change as we remove some of its nodes or links. The
properties we are interested in are typically topological, such as the size of
the largest connected component or the distribution of small connected
components, and they are used as a first proxy to assess the functionality of
a system exposed to failures or attacks. Failures are modelled as random
removals, whereas attacks assume some a priori knowledge of the network, where
the elements are ranked given some criterion —frequently topological (degree,
betweenness, etc.) [albert2000error, cohen2001breakdown], albeit not strictly
necessary [artime2020effectiveness, artime2021percolation]— and removed
accordingly. We aim at presenting the basic mathematical framework to address
the problem of multilayer percolation and showcase some applications. Good
reviews to expand on what we present here can be found in
[bianconi2018multilayer, lee2015towards].

The first step towards a description of multilayer percolation processes is to
consider the generalization of the generating function methodology, frequent
in single-layer networks [newman2001random]. In [leicht2009percolation],
random percolation is studied in general multilayer networks described by the
set of degree distributions $\\{p_{k_{1}k_{2}\ldots k_{L}}^{\alpha}\\}$, where
$\alpha$ is the layer label and $k_{\beta}$ denotes the number of links toward
nodes in layer $\beta$. Expressions for the size of the giant component of
each layer $S^{\alpha}$ can be written as a function of the generating
functions. This framework allows us to consider percolation on correlated
multilayer networks, leading to non-trivial results. For example, in
[lee2012correlated] correlated duplexes of Erdős–Rényi networks are studied.
When degrees across layers are maximally anti-correlated – i.e., hubs in one
layer are low degree nodes in the other layer – the giant component appears at
a link density considerably higher than the value for which it appears in
uncorrelated duplexes. The giant component exists, though, for any nonzero
link density if interlayer degrees are maximally positively correlated. In
other words, the more correlated the degrees are, the fewer edges are needed
to make a macroscopic structure emerge. When it comes to attacks on the most
connected nodes, the latter statement is true up to a certain intra-layer mean
degree (assuming it is the same for both layers), after which the behavior is
reversed and maximally negative correlated networks become more robust
[min2014network].

This framework disregards any functional characteristics of the layers and,
from a phenomenological point of view, continuous phase transitions are always
found. This is no longer true if the nature of the layers is taken into
account, for example via the interdependence of the nodes or antagonistic
interactions. These more realistic scenarios define new conditions for a node
to remain in the network, and need alternative, more function-oriented metrics
to describe the state of the system. A widely accepted choice is the mutually
connected component [buldyrev2010catastrophic, son2012percolation] already
defined in Sec. III.2, which is the set of nodes that are connected within
each and every layer simultaneously. The most striking result is that when the
giant mutual component (GMC) is computed in interdependent networks, the
percolation phase transition changes its nature, becoming a discontinuous one
[son2012percolation]. This has serious implications for the robustness of the
system, since the disintegration occurs abruptly, i.e., it is hard to
anticipate. Mathematically, the idea is as follows.

Let us assume an edge-colored multigraph, let $p_{k_{1}\ldots k_{L}}$ be the
probability that a node has degree $k_{\alpha}$ to other nodes within the
layer $\alpha$, and let $q_{k_{1}\ldots k_{L}}$ be the corresponding excess
degree distribution. We indicate by $w_{\alpha}$ the probability that a node
does not belong to the GMC via a link in layer $\alpha$. Hence,
$w_{\alpha}^{k_{\alpha}}$ gives the probability that the node does not belong
to the GMC via any of its neighbors in layer $\alpha$. The condition to belong
to the GMC is that the node has to be connected to it in all the layers, i.e.,
the size of the GMC is proportional to
$\prod_{\alpha=1}^{L}(1-w_{\alpha}^{k_{\alpha}})$. We just need to rescale by
the occupation probability $\phi$ and average over the degree distribution,
yielding

$\displaystyle M$
$\displaystyle=\phi\sum_{k_{1}=0}^{\infty}\ldots\sum_{k_{L}=0}^{\infty}p_{k_{1}\ldots
k_{L}}\prod_{\alpha=1}^{L}(1-w_{\alpha}^{k_{\alpha}}).$ (88)

To compute $w_{\alpha}$, we first note that $1-w_{\alpha}$ is the probability
that a node at the end of a link in layer $\alpha$ belongs to the GMC. For
this to happen any of its remaining $k-1$ neighbors in layer $\alpha$ are in
the GMC as well. Moreover, due to the condition of mutual connectivity, in
every other layer, the node needs to belong to the GMC via any of its
neighbors. Rescaling by $\phi$ because the node needs to be present in the
network, and averaging, we obtain

$1-w_{\alpha}=\phi\sum_{k_{1}=0}^{\infty}\ldots\sum_{k_{L}=0}^{\infty}q_{k_{1}\ldots
k_{L}}\left(1-w_{\alpha}^{k_{\alpha}-1}\right)\prod_{\begin{subarray}{c}\beta=1\\\
\beta\neq\alpha\end{subarray}}^{L}(1-w_{\beta}^{k_{\beta}}).$ (89)

By inserting the solutions $\\{w_{\alpha}\\}$ of this system of equations into
Eq. (88) we readily obtain the size of the giant mutual component. See Fig. 46
to appreciate the emergence of the abrupt transition of the GMC for multiplex
systems with Erdős–Rényi networks in the layers. The discontinuity is also
present when considering multiplexes of scale-free networks but, at odds with
single-layer scale-free percolation, the occupation probability $\phi_{c}$ is
finite, making them more vulnerable to random failures than the single-layer
network [baxter2012avalanche]. Targeted interventions in the network can be
easily modeled as well by including a degree-dependent occupation probability
$\phi_{k_{1}\ldots k_{L}}$ inside the sums [min2014network,
zhao2016robustness].

The dependence of the mutual component $M$ shown in Fig. 46 for Erdős–Rényi
multiplexes is shared for other intra-layer topologies as well. That is, when
the number of layers increases, the value of the mean degree at which the
discontinuity appears becomes larger, thus broadening the parameter region
where the network is not functional. Moreover, the height of the discontinuity
jump becomes larger too, thus making the transition harder to anticipate when
going from the supercritical to the subcritical region. In light of these
results, the more layers, the more fragile the system is, which might seem
paradoxical from an evolutionary point of view: why would a system organize
itself in a layered structure if that reduces its robustness? In
[radicchi2017redundant], Radicchi and Bianconi provided a potential answer to
this conundrum by proposing a model of multilayer percolation that relaxes the
condition for functionality of the nodes. They argue that a node can be
functional, i.e., it is not removed, as far as it is functioning in at least a
pair of layers. This new condition for functionality allows them to conclude
that the addition of extra layers boosts the robustness of the system.

Figure 46: Percolation phase transition in a multiplex formed by intralayer
Erdős–Rényi networks, with the same mean degree $\langle k\rangle$. On the
left, we display $M$, which is the solution of $M=\phi(1-e^{-\langle k\rangle
M})^{L}$, as a function of the mean degree, for different number of layers,
with the occupation probability fixed to $\phi=1$. On the right, we study the
emergence of the giant mutual component for partial multiplexes, where $q$ is
the fraction of interdependent nodes. In this case $M$ is the solution of
$M=\phi(1-e^{-\langle k\rangle M})(1-qe^{-\langle k\rangle M})$. In the inset
we show the height of the jump of the order parameter at the transition point,
along with the theoretical value of the tricritical point (vertical dashed
line).

The abrupt nature of the transition induced by Eq. (89) holds as far as $L>1$,
see Fig. 46. When $L=1$ the mutual component coincides with the standard giant
component, so the transition is continuous. We cannot interpolate continuously
from $L=1$ to $L=2$ to understand how the nature of the phase transition
changes. However, there are other variables that we can tune to go from an
effective single-layer system to a multiplex.

The first of these variables is the multiplexity parameter. It might occur
that in real multiplex networks only a fraction $q$ of all nodes shares the
functional dependency across layers, hence a natural question is what is the
nature of the transition as a function of $q$. Does it suffice to have a non-
zero fraction of dependent nodes to observe the discontinuity, or, on the
contrary, is there a finite threshold, only above which we observe an abrupt
transition? To answer this question, let us focus on duplexes. If a node in
one of the layers has a dependency link, which occurs with probability $q$,
then the condition to belong to the GMC is the same as discussed earlier.
Instead, if the node does not have a dependency link, which occurs with
probability $1-q$, then it belongs to the GMC as far as any of its neighbors
of its very same layer belong to the component. Therefore, we can write

$\displaystyle M_{\alpha}$
$\displaystyle=q\phi\sum_{k_{1}=0}^{\infty}\ldots\sum_{k_{L}=0}^{\infty}p_{k_{1}\ldots
k_{L}}\prod_{\beta=1}^{L}(1-w_{\beta}^{k_{\beta}})$
$\displaystyle+(1-q)\phi\sum_{k_{1}=0}^{\infty}\ldots\sum_{k_{L}=0}^{\infty}p_{k_{1}\ldots
k_{L}}w_{\alpha}^{k_{\alpha}}.$ (90)

The set of probabilities $\\{w_{\alpha}\\}$ obeys the equations

$\displaystyle 1-w_{\alpha}$
$\displaystyle=q\phi\sum_{k_{1}=0}^{\infty}\ldots\sum_{k_{L}=0}^{\infty}q_{k_{1}\ldots
k_{L}}\left(1-w_{\alpha}^{k_{\alpha}-1}\right)\prod_{\begin{subarray}{c}\beta=1\\\
\beta\neq\alpha\end{subarray}}^{L}(1-w_{\beta}^{k_{\beta}})$
$\displaystyle+(1-q)\phi\sum_{k_{1}=0}^{\infty}\ldots\sum_{k_{L}=0}^{\infty}q_{k_{1}\ldots
k_{L}}(1-w_{\alpha}^{k_{\alpha}-1}).$ (91)

Focusing again on a duplex of Erdős–Rényi networks, we see in Fig. 46 that
there is a finite tricritical point $q_{c}$ at which the transition changes
its order, confirming that, in general, a multiplex needs a certain level of
interdependency between layers to experience the discontinuous transition.
Since many infrastructural networks are embedded in a two-dimensional space,
in a similar fashion one might ask what type of transition we encounter when
coupling low-dimensional networks, that alone show continuous transitions, in
a multiplex. It turns out that in this case the transition does not change its
order [son2011percolation, berezin2013comment].

The second variable that allows us to interpolate between multiplexes and
monoplexes is the link overlap across layers. For complete overlapping all
layers are equal and the problem is reduced to percolation of a single-layer,
showing a continuous transition. How much overlap do we need to observe the
abrupt transition? Different approaches and techniques have been used to
answer this question, and accordingly slightly different phenomenology has
been discovered. Although the details of the phase diagram depend on the
number of layers and the degree distributions, Hu and coathors have found that
the transition stays abrupt as far as we do not have complete overlapping
[hu2013percolation]. In [cellai2013percolation] it has been reported that a
critical value of the edge overlap exists that changes the nature of the phase
transition from a hybrid first-order to a continuous one. See
[cellai2016message] for the generalization of the theory to an arbitrary
number of layers. If the edge overlap is combined with other topological
correlations, the phase diagram displays multiple and recursive hybrid phase
transitions [baxter2016correlated]. In [min2015link], the problem of link
overlap has been formulated in a way that the discontinuous transitions
display hysteresis. The role of the edge overlap, together with other
topological correlations, has been also explored in the problem of multilayer
optimal percolation [osat2017optimal], that is, in the identification of the
smallest set of nodes that, when removed, cause the largest damage in the
network [santoro2020optimal]. This problem is known to be NP-hard in single-
layer networks [kempe2003maximizing], and although there is no equivalent
proof in multilayer architectures, the intuition indicates that it is so as
well. On this basis, heuristic approaches to identify the critical subset of
nodes predominate. In [santoro2020optimal], the efficiency of $20$ of these
strategies is evaluated. The authors find that when no structural correlations
exist, a family of Pareto-efficient strategies based on both structural
descriptors and multiobjective optimization is the best at dismantling the
network. However, when evaluated in real multilayer networks that present non-
trivial correlations, the variability in performance changes from one
dismantling strategy to another, suggesting that a fair assessment of
multilayer robustness requires a comparison between strategies.

Beyond node interdependency as a condition for functioning, other interesting
mechanisms have been considered in the literature. One of them is antagonistic
relations, where a node in one layer is functional — that is, it has not been
removed — only if its replicas are not [zhao2013percolation]. Examples of this
kind of competitive or non-cooperative relations might be relevant, for
example, in biological networks. Interestingly, it has been found that when
percolation is considered in this scenario, the abrupt transition shown in
what follows persists, but displays as well, at variance with the transition
of the mutual component introduced earlier, the typical hysteresis and bi-
stability behaviors of equilibrium thermal first-order phase transitions
[goldenfeld2018lectures]. Other mechanisms have followed to include these
behaviors too, e.g., in [min2014multiple, min2015link].

Note that all of these results derived with the machinery of generating
functions come with some underlying assumptions that are very rarely met when
studying real networks: (i) the network is tree-like, (ii) it has infinite
size and (iii) the quantities of interest are computed as an average over the
ensemble of networks with given degree distribution, although, in reality, we
only have access to one supra-adjacency matrix, among all the possible ones
that a graph model could generate. As a consequence, the analytical
predictions might deviate from the actual process of percolation, obtained for
example via simulations. Analytical approaches have been proposed to overcome
points (ii) and (iii) in [radicchi2015percolation, bianconi2016percolation],
where a percolation theory of multiplex and interdependent networks is
introduced that takes as input the adjacency matrix instead of the degree
distribution. Furthermore, in a real network we may need to assess the
robustness under perturbation scenarios related to node metadata. For example,
in [baggio2016multiplex] it is addressed, among other things, how a multiplex
network capturing the flow of subsistence-related goods and services among
households in several Alaskan Natives communities responds to perturbations
involving targeted removals of specific resources by category, e.g.,
terrestrial, marine or riverine, as a representation of natural disasters. An
analytical framework to take into account non-topological features in the
robustness of a network has been recently developed for single layers
[artime2021percolation], but, at present, there is not a generalization for
multilayer networks.

###  IV.6 Cascade failures

The percolation model has the limitation that failures and attacks are treated
statically. A more complete description of multilayer robustness and
resilience would include their time evolution, in order to better understand
under which conditions small perturbations can trigger global network-wide
effects, the so-called _cascades_. In this section we review some of the most
emblematic models for multilayer cascade failures.

The seminal work of Buldyrev and colleagues was one of the first to deal with
the propagation of failures across interdependent multiplexes
[buldyrev2010catastrophic]. The motivation for such a study was the 2003 Italy
blackout, for which it was hypothesized that the power grid and the Internet
network, the latter acting as a supervisory control and data acquisition
system, were interdependent, and that failures in the power stations hampered
the Internet communication and further propagated the malfunction across the
system, see top panels of Fig. 47).

However, the interest in the amplification of small perturbations throughout
the network is much more general, finding eventual applications in areas such
as biochemistry, where metabolic pathways interact in complex ways with key
tissues, or finance, where different banks might share the same asset in their
balance sheet [huang2013cascading], among many other examples. In
[buldyrev2010catastrophic] was introduced the concept of a mutually connected
component from a dynamical perspective, at odds with its static definition
presented in the previous section, see bottom panels of Fig. 47. It was shown
that these cascades evolving in multiplex networks yield a discontinuous phase
transition in the size of the giant mutual component
[buldyrev2010catastrophic], meaning that the failure of one single node from
an apparently healthy, functional infrastructure can generate the emergence of
a global cascade that collapses the entire system. It was later shown that
these cascading failures can be mapped to static percolation in multiplexes
[son2012percolation]. In fact, many of the results presented in Sec. IV.5 are
recovered at the final state of the cascade propagation of this type of model.

Figure 47: Top: evolution of a cascade of failures in a real interdependent
system composed by a power network (on the map) and an Internet network
(shifted above) and involved in the Italian blackout of September 2003. In a,
failure of a power station makes the Internet nodes to stop functioning,
indicated in red. In green is shown the nodes that disconnect from the largest
connected component in the next cascade step. These nodes will fail
immediately after ($b$), introducing a feedback of failures back in the power
network. Those nodes that have been isolated (green nodes in $b$) are the ones
that will fail in the next event $(c)$, inducing further interdependent
failures in the Internet network. Bottom: we sketch this process for further
clarification. In d, an initial attack to one node in network $\mathsf{A}$
occurs. In e, the attacked node is removed, along with its links, and the
corresponding interdependent node in network $\mathsf{B}$, with its links, is
also removed. In f, the actual cascade starts. All the $\mathsf{B}$-links
between $\mathsf{B}$-nodes connected interdependently to different
$\mathsf{A}$-clusters are removed. In d, the same rule applies but for links
in network $\mathsf{A}$: all $\mathsf{A}$-links between $\mathsf{A}$-nodes are
deleted if the interdependent nodes do not belong to the same
$\mathsf{B}$-cluster. Steps f and g are repeated, propagating the failures
back and forth, until the cascade cannot further evolve. The remaining
connected components at the end of the propagation of failures coincide with
the mutually connected components introduced in Sec IV.5. Figures from
[buldyrev2010catastrophic].

Despite the fact that the dynamical propagation of cascades needs a more
convoluted mathematical treatment than percolation, these techniques have been
very flexible at adapting to variations of the original work of Buldyrev, so
that the propagation of cascades is modeled in more realistic scenarios. We
briefly review some of them in the following, and refer the interested reader
to the more complete reviews [kenett2014network, shekhtman2016recent,
valdez2020cascading].

After [buldyrev2010catastrophic], the problem of dependency-based cascades has
been generalized to $L$ layers in [gao2012networks], allowing analytical,
closed solutions in certain interdependent setups. The problem of failure
propagation in multiplexes with partial interdependency, where not all nodes
have dependency connections, has been addressed in [Parshani2010reducing]. If
the fraction of interdependent nodes is decreased enough, the transition is no
longer abrupt but becomes second-order. In order to narrow the range of
parameters for which the abrupt transition occurs, different strategies to
choose the autonomous nodes, those without interdependencies, have been
proposed [schneider2013towards]. It turns out that selecting those with the
highest degree or highest betweenness, significantly reduces the likelihood of
an abrupt collapse [schneider2013towards]. Intra-layer correlations can be
included in the analysis as well. In [huang2013robustness], interdependent
multiplexes with a tunable average number of single links and an average
number of triangles per node are analyzed, finding that, for fixed average
degree, the higher clustered networks are less robust than those with lower
clustering. Intra-layer correlations might be coupled as well to inter-layer
ones: in [reis2014avoiding] a flexible model including both types of
correlations is studied, and the maximization of robustness is addressed as a
function of the tunable strength of the correlations and the failure
propagation dynamics, showing that the theoretical results match well with the
experimental results of coupled functional brain modules.

Beyond cascades driven by topological failures, there are other mechanisms via
which small perturbations can drive a multilayer network to collapse. In many
problems related to social sciences, such as the adoption of fads, the
diffusion of norms and innovations, and the changes in the collective
attention in a population, models of behavioral contagion are used to model
the decisions of the agents. Agents need to decide between two alternative
options/actions and the influence of their neighbors is crucial in the final
choice. A simple way to encompass this influence is by setting an activation
threshold. Initially all nodes are inactive but a controlled small fraction.
An inactive agent changes her state and becomes active only when the fraction
of her active neighbors is larger than a threshold. This process might require
several activation steps before reaching a frozen state. In [watts2002simple]
was given the range of parameters for which such global cascades, measured as
the number of active users when the process stops, emerge in monoplexes,
turning out that they are only possible when the network is neither too sparse
nor too dense. The generalization of this threshold model in multiplexes was
addressed in [brummitt2012multiplexity], where a propagation rule is proposed
such that an agent activates as far as the fraction of active neighbors in at
least one layer exceeds the threshold. Following this rule, multiplexity
widens the ranges of parameters in which global cascades can occur, therefore
increasing the network’s vulnerability. Interestingly, isolated layers in
which, owing to their topological properties, global cascades would not be
observed, when multiplex-coupled, cooperate in a way that facilitates agent
activation and therefore lead to cascades.

Figure 48: The power grid is a paradigmatic example of a network susceptible
to overload cascades of failures due to the energy transported along the
links. We see on the top an aggregated representation of the multilayer power
grid of the USA, where each layer (link color) corresponds to different
voltage ranges. In red is depicted the planned interconnections to transport
wind power. On the bottom, the Bak-Tang-Wiesenfeld sandpile model is simulated
in a duplex of random regular graphs. For a cascade starting in one of the
networks (network $a$), is shown the probability of observing a final cascade
size in $a$, $T_{aa}$, and in network $b$, $T_{ab}$, as a function of the
probability of interconnections $p$ between nodes in $a$ and in $b$. Both
$T_{aa}$ and $T_{a}=1/2(T_{aa}+T_{ab})$ display a minimum, indicating that an
isolated network can reduce the damage caused by a cascade by interconnecting
to another network, but only up to a certain level of interconnectedness
$p^{*}$. Figures readapted from [brummitt2012suppressing]. Figure 49:
Robustness of multilayer networks exposed to overload failures due to nonlocal
load redistribution. We show the behavior of $\langle S\rangle$, the size of
the largest connected component of the network at the end of the cascade, as a
function of the tolerance parameter $\alpha$ (see text). The faster $\langle
S\rangle$ grows, the more robust is the system because there the range of
tolerance values for which it disintegrates is smaller. Three mechanisms to
increase the robustness are presented. In $(a)$, the comparison between the
multilayer network and its aggregated counterpart. In $(b)$, dependence on the
number of layers. In $(c)$, $\langle S\rangle$ for different values of the
multiplexity parameter (fraction of nodes participating simultaneously in a
duplex), with the value indicated in the legend. Figure readapted from
[artime2020abrupt].

Another mechanism that can induce cascades is load redistribution due to
overloads. Arguably, the most famous stylized model accounting for this type
of dynamics is the Bak-Tang-Wiesenfeld sandpile model [bak1988self], a
paradigmatic example of self-organized criticality used in a variety of
contexts [jensen1998self]. Each node has an internal, discrete variable, the
load. At each unit of time the load of a node selected uniformly at random
increases in one unit. When the load reaches a threshold, assumed to be
degree, the node redistributes its load to its neighbors. The neighbors might,
in turn, exceed their capacity and reallocate their load to their own
neighbors, hence propagating the cascade. Once all nodes have their own load
below the threshold, the random addition of load again is restarted. To avoid
inundation of the system, a frequently used strategy is to dissipate load at a
certain rate when it is reallocated. In [brummitt2012suppressing], the BTW
model is used in the context of multilayer power grids, shedding light on the
benefits and dangers of the level of interconnectivity between the layers.
Surprisingly, it is found that there is an optimal value of the
interconnectivity for which the chance to observe large cascades is minimum,
see Fig. 48. This is because the addition of interlayer links between highly
isolated power grids helps mitigate the cascades by creating a reservoir that
absorbs excess load, a sort of alternative dissipative mechanism. This
behavior, though, is valid up to a certain point, from which the further
addition of interlayer links is no longer beneficial since it creates more
loaded systems due to larger thresholds and, therefore, chances of larger
cascades, as well as it creates a positive feedback due to the reentrant new
paths of load into a layer. Regarding the critical properties of the BTW
model, it has been found that multiplexity does not alter the mean-field
scaling behavior observed in monoplexes [lee2012sandpiles].

To close this section, we discuss the relevant, but still quite unexplored,
case of nonlocal cascade propagation. One might argue that nonlocality can be
realized in spatially extended systems by allowing interlayer dependencies
between different locations of the space [li2012cascading]. Yet, under this
description the failures propagate via first neighbors, in the interdependent
sense. In fact, all types of cascades discussed earlier evolve locally via
first-neighbors, something that need not to be true in real systems, as
happened in the 1996 disturbance of the Western Systems Coordinating Council
(WSCC) system [NERC2002system], in the 2003 blackout in the northeastern USA
[nerc2004technical] or in the air-traffic disruption due to the eruption of
the Icelandic volcano Eyjafjallajökull [eurocontrol2010ashcloud]. A plausible
description of this phenomenon is based on the load-capacity model of Motter
and Lai [motter2002cascade], where a load, defined as the number of shortest
paths crossing a node, and a constant capacity, a factor $1+\alpha$ larger
than the initial load, are assigned to every node. When the load of a node
exceeds its capacity, the node fails. An initial perturbation in the form of
node removal is applied to the network, that globally modifies the loads and
allows subsequent failures to occur not necessarily close to the prior
failure. If during this process nodes that overload, they also fail,
propagating the cascade. This model has been investigated recently in
[artime2020abrupt], with the finding that the size of the largest connected
component at the end of the cascade suffers an abrupt jump when the tolerance
parameter $\alpha$ is increased. Moreover, since the load redistribution
depends significantly on the topology of the network, the average path length
is identified as a metric that correlates well with robustness. Based on this,
the article proposes different strategies to increase the robustness of the
network, such as adding new layers or reducing the level of multiplexity (see
Fig. 49). Unlike the other cascades phenomena introduced in this section,
which can be analytically treated with generating functions or multi-type
branching processes, nonlocal cascades do not have a solid mathematical
machinery to be described, and this certainly represents a challenge for the
future.

##  V Frontiers

###  V.1 Kinematic geometry

The geometric approach [boguna2020network] to network analysis has garnered
growing interest in the past two decades. Here we focus on the geometry of
network-driven processes on multilayer networks, a family of kinematic
geometries generalizing the diffusion geometry introduced in
[de2017diffusion].

The structure of a wide variety of real-world complex systems is modular and
hierarchical [guimera2005functional] and the effect of these large scale
properties on the dynamics of such systems has been studied, during the past
decade. It has been shown that complex systems with such a mesoscale
organization [fortunato2010community] are characterized by topological scales
[arenas2006synchronization], exhibiting the emergence of functional clusters
which might be different from topological ones.

In [de2017diffusion] the authors investigate the multiscale functional
geometry of monoplexes to characterize functional clusters. This approach
defines the diffusion distance between any pair of units in a networked
system, based on random walk dynamics, shown to correspond, among others, to
the phase deviation of coupled oscillators close to metastable synchronization
state and consensus dynamics. The diffusion distance for single-layer networks
is a key concept to define a kinematic geometry based on network-driven
processes and it can be calculated as the $L^{2}-$norm of the difference
between rows of the propagator $e^{-t\tilde{\mathbf{L}}}$:

$D_{t}(i,j)=\left\lVert\mathbf{p}(t|i)-\mathbf{p}(t|j)\right\rVert_{2}.$ (92)

Exploiting the fact that diffusion geometry is based on random walk dynamics,
it is possible to extend it to the realm of multilayer networks. At variance
with walks on edge-colored networks presented in Sec. IV.1, where one can
obtain the transition rules for the multigraph as a weighted average of the
transition probabilities in each distinct layer, see Eq. (61), on a multilayer
network the walk type determines the probability of a random walker jumping
across and within layers (see Eq. 62 and Tab. 1). Consequently, in the edge-
colored case diffusion distances between nodes can be obtained directly from
Eq. (92), while for multilayers we have to introduce a diffusion distance
between state nodes.

Regardless of the random walk type, let us indicate the probability of finding
a random walker at a given node and layer at time $t$ by $p_{j\beta}(t)$.
Similarly to Eq. (92) we define the diffusion distance between state nodes
$(i,\alpha)$ and $(j,\beta)$ as

$D_{t}^{2}((i,\alpha),(j,\beta))=\sum\limits_{k,\gamma}(p_{k\gamma}(t|(i,\alpha))-p_{k\gamma}(t|(j,\beta)))^{2},$
(93)

where probabilities are conditional to using the corresponding state nodes as
the origins of random walkers at time $t=0$. One can summarize this supra-
distance matrix $\mathbf{D}_{t}=\left(D_{t}((i,\alpha),(j,\beta))\right)$ into
an $N\times N$ matrix by encoding the diffusion distance among the physical
nodes, across all the layers. Intuitively, resembling the parallel sum of
resistances in electrical circuits leading to the equivalent resistance, the
_equivalent diffusion distance_ can be written as

$D_{t}^{\text{eq}}(i,j)=\left(\sum_{\alpha=1}^{L}\frac{1}{D_{t}((i,\alpha),(j,\alpha))}\right)^{-1}.$
(94)

Usually, the functional distance (supra-)matrices are rescaled in $[0,1]$,
normalizing each by its maximum value to allow for comparisons. If,
additionally, one is interested in the most persistent patterns, the diffusion
distance is averaged over time and called the average diffusion distance. The
corresponding supra-distance matrix is indicated by $\bar{\mathbf{D}}_{t}$.

The diffusion distance between nodes is highly dependent on the type of random
walk dynamics, its propagation time, the topology of layers and the layer-
layer correlations. To better understand the effects of dynamics and
topological variations, let us consider three classes of two-layer networks,
namely:

  * •

Barabasi-Albert layers with preferential attachment of 4 links;

  * •

Watts-Strogatz layers, with rewiring probability 0.2;

  * •

Girvan-Newman model layers, where the inter-community connectivity probability
is 1, whereas cross-group connections exist with probability 0.05.

Additionally, we consider the five random walks of Tab. 1. For the three
cases, the link overlap across layers – i.e., the fraction of links present in
both layers among the same pairs of nodes [de2015structural]— is fixed to 10%.
As shown in Fig. 50, the diffusive walk on a scale-free topology leads to a
high level of mixed pathways across layers, while in small-world systems,
nodes tend to shape more distinct functional clusters.

Figure 50: Average diffusion distance supra-matrices $\bar{\mathbf{D}}_{t}$
for different combinations of multilayer topologies and random walk dynamics
(see the text for details). Figure reprinted with permission from
[bertagnolli2020diffusion]. Copyright (2021) by the American Physical Society.

The functional geometry framework has also been used to analyze empirical
systems in [bertagnolli2020diffusion], where the authors highlighted
dissimilarities in the diffusion spaces of the public transportation of London
and of the social network of Noordin terrorists, by looking at the projections
of the spaces and at the results of a Mantel’s test [mantel1967detection]
applied to supra-distance matrices.

Finally, it is worth mentioning that one can use other walks instead of random
ones: this is the case of walks that result from powers of the adjacency
matrix and provide the basis for _communicability_ , a measure recently used
to introduce a geometric framework with applications to single and multiplex
networks [Estrada2019].

###  V.2 Statistical physics of multilayer systems

####  V.2.1 Classical ensembles

To study the properties of an observed multilayer network, comparison with
null models is often essential. The maximum entropy approach is one of the
standard ways to obtain the required null models, in terms of ensembles of
networks exhibiting one (or more) specific property of the observed network,
while being maximally random with respect to the other properties. For
instance, in the case of single layer networks, one of the most famous null
models is the configuration model (CM) [molloy1995critical] which is an
ensemble of networks with the same degree sequence as the observed one. The
unbiased probability distribution of the members of this ensemble is indicated
as $P(G)$ and must maximize the Shannon information entropy

$S=-\sum\limits_{G\in\mathcal{G}}P(G)\log{P(G)}.$ (95)

It is worth mentioning that such a probability distribution does not
correspond to a physical process related to the second law of thermodynamics.
However, the discussed entropy maximization approach is mimicking the
mathematical machinery of statistical physics and, in the case of a fixed
degree sequence, it leads to a microcanonical ensemble where all members of
the CM have equal probability.

Depending on the type of study, the mentioned constraint of a fixed degree
sequence can be relaxed to obtain other ensembles [cimini2019]. For instance,
by fixing the average degree instead of the full degree sequence, we achieve
the grand-canonical ensemble of random graphs. Furthermore, if our knowledge
is limited to the degree distribution from which the degree sequence is
sampled, one obtains the hypersoft configuration model
[Krioukov2020MaxEntropy] which is a hypercanonical ensemble of random networks
all drawn from the fixed distribution (see Fig. 51).

Figure 51:  Schematic of sampling the network, represented by $G$, from
hypercanonical, grandcanonical, and microcanonical ensembles discussed in the
text. Figure from [Krioukov2020MaxEntropy].

Similarly, maximum entropy approaches have been used to study ensembles of
non-interconnected multiplex networks. Assume an edge-colored multigraph $M$
with $L$ layers, each forming a network $G^{(\ell)}$ ($\ell=1,2,...N$).
Remarkably, when there is no correlation between the layers, the probability
of observing $M$ can be obtained as the product of the probabilities of the
layers:

$\displaystyle P(M)=\prod\limits_{\ell=1}^{L}P(G^{(\ell)})$ (96)

and the corresponding Shannon entropy becomes the summation of Shannon
entropies of layers [bianconi2013statistical]:

$\displaystyle S=\sum\limits_{\ell=1}^{L}S^{(\ell)}.$ (97)

By imposing the soft constraints – e.g., fixing the average degree instead of
degree sequence – and using the Lagrangian multipliers method, one can obtain
the canonical multiplex ensemble with probability $P_{C}(M)$ maximizing the
Shannon entropy:

$\displaystyle
P_{C}(M)=\frac{e^{-\sum\limits_{\mu}\lambda_{\mu}F_{\mu}(M)}}{Z_{C}}=\prod\limits_{\ell=1}^{L}P_{C}(G^{(\ell)}),$
(98)

where $\lambda_{\mu}$ corresponds to the Lagrangian multipliers and
$F_{\mu}(M)$ determine the constraints on the network (e.g., average degree).
Note that, in presence of layer-layer correlations, we have
$P(M)\neq\prod\limits_{\ell=1}^{L}P(G^{(\ell)})$ leading to other probability
distributions extensively studied in [bianconi2013statistical].

####  V.2.2 Quantum-like ensembles

Complex systems include a wide variety of physical attributes and dynamics.
The interactions between their units vary, from the electrochemical signals
traveling among neurons in the human brain to the transport of goods between
different areas of a urban ecosystem and the spreading of an infectious
pathogen between individuals of a society. Regardless of the nature of these
interactions, they can be described as information exchange. In fact, complex
systems resemble one another in the way they handle information. Therefore, to
understand how these systems operate, from a physical point of view,
investigating their information dynamics is crucial.

It is important to note that the information flow between the units is
regulated by the underlying structure, depending on the local neighborhood in
certain classes of networks and on long-range communication between distant
components in other topologies. At the same time, it is essential to consider
the coupling between the structure and the dynamical processes governing the
flow of information, such as diffusion, random walks, synchronization and
consensus and, since it is often difficult to understand a system in terms of
microscopic features, a framework providing a statistical description of the
system might be relevant for applications.

Recently, a statistical field theory of complex information dynamics has been
introduced to unify a range of dynamical processes governing the evolution of
information on top of static or time-varying structures [de2020SFT]. This
framework describes the interactions among the units in terms of _information
streams_ , a set of operators that determines the direction of information
flow and provides a statistical ensemble to construct the statistical physics
of complex information dynamics. In fact, the formalism allows for defining
density matrices— i.e., statistical average of streams— from which a variety
of descriptors can be derived. Of course, density matrix formalism comes from
quantum mechanics, where vectors were found insufficient to represent the
mixed states and encode the pairwise coherence between quantum states.
Similarly, properties of networked systems cannot be fully described by
vectors or distribution functions, without information loss. For instance, a
system’s structure is often encoded in two dimensional data structures, being
the adjacency matrix, except for trivial symmetries like chains or paths.
Therefore, the counterpart of the quantum density matrix has been introduced
for complex networks, where off-diagonal elements provide a proxy for
interactions between the nodes, instead of the coherence between the states.
So far, this framework has only been used to analyze classical networks. In
the following, we provide a direct application of the theoretical framework to
multilayer systems.

Figure 52:  Schematic illustration of the spectral entropy for different
classes of networks has been presented. Remarkably, the entropic distances
provided by the framework can be used to compare network similarity and
characterize the network topology with high accuracy. Figure from
[de2016spectral].

##### Redundancy and reducibility of multilayer systems

As discussed previously, the constituents of complex systems must exchange
information efficiently in order to function properly. However, a deep
understanding of the structure’s role in enhancing or hindering the transport
properties – such as navigability [boguna2009navigability] – continues to
elude us, especially in the case of multilayer systems [gao2012networks,
de2013mathematical, gomez2013diffusion, radicchi2013abrupt,
kivela2014multilayer, boccaletti2014structure, de2016physics] where an
efficient flow might be hindered by the lack of synchronization between
different layers [gallotti2014anatomy] or the redundancy of pathways across
the layers. There are a number of methods to reduce multiplex networks by
structurally merging their most similar layers, based on information-theoretic
frameworks [de2015structural]. Despite their success, these methods are mostly
based on heuristics and have been proved inaccurate under specific
circumstances [de2020EnhanceTransport].

Enhancing the flow distribution in multilayer systems is challenging, since
adding links to the layers (e.g., highways, tubes, flights, synapses, etc.)
comes with a cost. Interestingly, for multilayer networks with interlinks,
changing the weights of interlinks can enhance the diffusion on top of the
networks [gomez2013diffusion]. When acting on the structure is not an option,
it has been shown that one can still enhance the transport properties (see
Fig. 53) using the statistical physics of complex information dynamics by
coupling layers dynamically, in a way that a dynamical process can not
distinguish the functionally coupled layers – e.g., airlines proposing shared
flights to their customers – while evolving. The functional reduction includes
coupling the layers with high similarity that are responsible for the
redundant diffusion pathways in the system, in order to identify the (sub)set
with maximally diverse layers [de2020EnhanceTransport].

Figure 53:  Schematic illustration comparing structural against functional
reduction of a multilayer network consisting of $L=4$ layers. The procedure is
similar for both approaches, but the relevant difference is that while
structural reducibility alters the topology of the system, function
reducibility allows us to functionally couple layers without altering their
structure. Figure from [de2020EnhanceTransport].

Diffusive processes, such as random walks, have been used extensively to model
the information transport within complex structures [masuda2017random]. Here,
we consider random walk dynamics governed by the normalized Laplacian
$\hat{\tilde{\mathbf{L}}}=\langle\hat{\tilde{\mathbf{L}}}^{(\ell)}\rangle$
given by
$\hat{\tilde{\mathbf{L}}}=\hat{\mathbf{I}}-\langle\hat{\mathbf{T}}^{(\ell)}\rangle$
playing the role of the quasi-Hamiltonian (see Eq. (61)). The density matrix
can be obtained for random walks on multiplex networks as
$\hat{\boldsymbol{\rho}}(t)=\frac{e^{-t\hat{\tilde{\mathbf{L}}}}}{Z(t)}$.

Interestingly, it has been shown that the partition function, which is
encoding the amount of the trapped field, is proportional to the average
return probability of random walk dynamics: $Z(t)=N\mathcal{R}(t)$, where
$\mathcal{R}(t)=N^{-1}\sum\limits_{i=1}^{N}e^{-t\lambda_{i}}$ is the average
return probability and $\lambda_{i}$ is the $i$–th eigenvalue of
$\hat{\tilde{\mathbf{L}}}$. Intuitively, the average return probability is
high when the structural symmetries and abundance of redundant diffusion
pathways slows down the information propagation between the units. Thus, it is
expected that breaking the structural regularities by adding long-range
interactions [watts1998collective] or increasing the diversity of diffusion
pathways across layers can lead to faster information flow and and lower
$Z(t)$.

Multiplexity of interactions among units generates non-trivial dynamical
correlations between layers that have no counterpart when layers are
considered in isolation. This important difference can be characterized in
terms of average entropy distance between the multiplex and its layers. This
measure, named intertwining, is defined by

$\displaystyle\mathcal{I}$ $\displaystyle=$
$\displaystyle\langle\mathcal{D}_{KL}(\hat{\boldsymbol{\rho}}||\hat{\boldsymbol{\rho}}^{(\ell)})\rangle=\frac{1}{L}\sum\limits_{\ell=1}^{L}\mathcal{D}_{KL}(\hat{\boldsymbol{\rho}}||\hat{\boldsymbol{\rho}}^{(\ell)})$
(99)

where
$\mathcal{D}_{KL}(\hat{\boldsymbol{\rho}}||\hat{\boldsymbol{\rho}}^{(\ell)})=Tr[\hat{\boldsymbol{\rho}}(\log_{2}\hat{\boldsymbol{\rho}}-\log_{2}\hat{\boldsymbol{\rho}}^{(\ell)})]$
is the quantum-like Kullback-Leibler (KL) divergence between layer $\ell$ and
a multiplex system as a whole.

Directly from intertwining, a fundamental inequality between the partition
function of a multiplex system as whole and the partition functions of its
layers can be derived. This inequality is important, as it relates the
transport phenomena of multiplex system and layers, through average dynamical
trapping (i.e., the partition function):

$\displaystyle Z(t)\leq\prod\limits_{\ell=1}^{L}Z^{(\ell)}(t)^{1/L},$ (100)

where equality holds if and only if all the layers are the same. Using
dynamical trapping as a measure of transport, Eq. (100) shows that a multiplex
network has better transport properties than the geometric average of layers,
adding an advantage to multilayer structures.

Furthermore, being reminiscent of statistical physics of particles, the
equality defines the non-interacting scenario where layer-layer correlations
do not alter the underlying dynamics: the entropy $S^{(\ell)}(t)$ of each
layer is calculated separately and the overall entropy is given by their
average $\mathcal{S}_{nint}(t)=\langle S^{(\ell)}(t)\rangle$. Conversely, any
topological alteration of the non-interacting scenario introduces a dynamical
correlation between layers, requiring the exploration of layers to gather more
information about the system: in this case the network consists of interacting
layers, where the diffusion dynamics on the whole multiplex network is
considered to measure the entropy $\mathcal{S}_{int}(t)$. To obtain another
form of Eq. (99), a mean-field approximation of the Von Neumann entropy is
given by

$\displaystyle\mathcal{S}^{MF}(t)=\frac{1}{\log
2}\left(t\frac{Z(t)-1}{Z(t)}+\log Z(t)\right),$ (101)

and can be used to prove that layer-layer interactions lower the system’s
entropy ($\mathcal{S}_{int}(t)\leq\mathcal{S}_{nint}(t)$) and partition
function $Z(t)$. Normalizing the intertwining by its upper bound, for values
of time $t$ sufficiently large and in absence of isolated state nodes, Eq.
(99) reduces to the relative intertwining:

$\displaystyle\mathcal{I}^{*}(t)=1-\frac{\mathcal{S}_{int}(t)}{\mathcal{S}_{nint}(t)},$
(102)

which is bounded between 0 (i.e., the layers are redundant) and 1 (i.e., the
layers are diverse and the system is irreducible).

We can show how intertwining is proportional to the functional diversity of
layers. The Laplacian matrix of the multiplex network is given by
$\hat{\tilde{\mathbf{L}}}=\langle\hat{\tilde{\mathbf{L}}}^{(\ell)}\rangle$.
Therefore, the Laplacian matrix of each layer
$\hat{\tilde{\mathbf{L}}}^{(\ell)}$ can be written as a perturbation of
multiplex Laplacian
$\hat{\tilde{\mathbf{L}}}^{(\ell)}=\hat{\tilde{\mathbf{L}}}+\Delta\hat{\tilde{\mathbf{L}}}^{(\ell)}$,
reflected in its eigenvalues as
$\lambda_{i}^{(\ell)}=\lambda_{i}+\Delta\lambda_{i}^{(\ell)}$ $(i=0,1,...N)$.
It is straightforward to show that
$\frac{1}{N}\sum\limits_{i=1}^{N}\Delta\lambda_{i}^{(\ell)}=\overline{\Delta\lambda^{(\ell)}}=0$
and that $\overline{\Delta\lambda^{(\ell)^{2}}}\geq 0$, the latter quantifying
the influence of the perturbation to each layer. The average of the variance
across all layers $\overline{\langle(\Delta\lambda^{(\ell)})^{2}\rangle}$
provides a measure of the overall spectral diversity of layers which,
interestingly, is demonstrated to be proportional to the relative
intertwining:

$\displaystyle\mathcal{I}^{\star}(t)$ $\displaystyle\approx$
$\displaystyle\frac{t^{2}}{2}\overline{\left\langle(\Delta\lambda^{(\ell)})^{2}\right\rangle}.$
(103)

This proves the sensitivity of intertwining to the spectral diversity of
layers. Furthermore, the partition functions of layers can be written in terms
of perturbations such as $Z^{(\ell)}(t)=Z(t)+\Delta Z^{(\ell)}(t)$, leading to

$\displaystyle\mathcal{I}^{\star}(t)$ $\displaystyle\approx$
$\displaystyle\frac{\overline{\Delta Z^{(\ell)}}(t)}{Z(t)-1};\
\overline{\Delta Z^{(\ell)}}(t)=\frac{1}{L}\sum\limits_{\ell=1}^{L}\Delta
Z^{(\ell)}(t).$ (104)

Equations (103) and (104) provide a fundamental result: they show that by
minimizing the partition function of the system one maximizes the relative
intertwining while favoring the maximum functional diversity of layers.
Additionally, it has been shown that an inverse proportionality holds between
the diffusion time ($1/\lambda_{2}$) and intertwining, as further evidence for
the role of intertwining in characterizing the transport properties.

These theoretical findings have been applied to a broad range of synthetic and
empirical systems, providing a transparent framework for coupling the most
similar layers of multiplex networks, in order to improve their transport
properties including dynamical trapping, diffusion time, and navigability
[de2020EnhanceTransport].

##  VI Conclusions

Network Science is one of the greatest achievements of the 21st century,
paving the way toward a mathematical approach for the analysis of disparate
complex systems, and allowing us to find regularities in apparently disordered
connectivity patterns. The past decade has seen the flourishing of analytical
techniques and models exploiting, or characterizing, the inherent
multidimensionality of empirical systems, from multiplexity – i.e., the
existence of distinct types of relationships among the same set of actors or
units – to interdependency – i.e., the existence of structural or functional
connections among sets of actors or units of a different nature. Such multiple
dimensions are nowadays easily encoded into layers of information.

Figure 54: Multilayer representation of a cell. Layers: i) regulatory
interactions involving RNA and protein expression, ii) protein-protein
interactions involved in signaling and responsible cell function, iii)
metabolic interactions with reactions and pathways crucial for cell function.
Figure from [vermeulen2020exposome].

Most of the advances in this direction are described in some detail or
referred to in this work. Starting from the mathematical representation of
multilayer networks (Chap. 2) we have introduced structural descriptors for
units, layers and the whole system (Chap. 3), providing an overview of the
micro- and meso-scale organization of such systems. We have discussed the rich
spectrum of phenomena, with no classical counterparts, related to dynamical
processes on the top of the networks and their intertwining (Chap. 4), guiding
the reader towards two promising research areas for the future, namely network
geometry and information dynamics (Chap. 5), although many other exciting sub-
field are emerging, e.g., higher-order modeling and analysis
[lambiotte2019networks, battiston2020networks].

At this point, the reader should be sufficiently familiar with multilayer
network science and, to conclude, we would like to make a quick journey
through the most recent applications of its paradigm, moving across different
spatial scales and ranging from cells to societies.

The first stop of this journey is exactly a cell, which can vary in diameter
between $10^{-6}$ m and $10^{-4}$ m (note that a DNA double helix is about
$10^{-8}$ m wide). The cell can be seen as a multilayer system consisting of
three interdependent layers (see Fig. 54). Here, applications are mostly
related to the emerging field of systems biology and network medicine,
promising to use biomolecular interactions to develop a deeper knowledge of
biology across scales with the ultimate goal to better understand diseases, to
prevent them, and to treat them with medical drugs that reduce side effects.

The second stop requires a big jump of about three orders of magnitude, to
explore one of the most famous multicellular organisms with a small-scale
neural system: the Caenorhabditis elegans, a nematode of about $10^{-3}$ m.
The nervous system of this small worm consists of synaptic and neuropeptide
interactions which can be seen as interconnected layers. Here, the multilayer
perspective provides the opportunity to better understand how the integration
between hard-wired synaptic or junctional circuits and extrasynaptic signals
can modulate the large-scale behavior of the worm [bentley2016multilayer].

At a larger scale, around $10^{-1}$ m, another emblematic neural system,
namely the human brain, is currently being characterized in terms of how its
structural and functional connectivity evolves over time, e.g., while
performing a specific task, or across groups, unraveling the existence of
modular and hierarchical structures which would remain hidden under the lens
of less sophisticated models and analytical techniques [bassett2017network].
In parallel, the functional connectivity of the human brain can be stratified
by frequency bands (see Fig. 55) where specific correlations or causal
relationships between regions of interests appear [de2017multilayer].
Multilayer analysis can be used to better characterize the relative importance
of all brain regions [williamson2021multilayer] and to enhance the accuracy in
discriminating between healthy and unhealthy patients starting from imaging
information [de2016mapping].

Figure 55: Multilayer representation of a human brain. Both the 3D and layered
visualization encode the functional brain of a schizophrenic subject (11 non-
overlapping frequency-band layers between 0.01 and 0.23 Hz). Figure from
[de2017multilayer], readapted from [de2016mapping].

At larger scales, on the order of hundreds of meters ($10^{2}$m), cooperative
systems like that of dolphins, can be characterized with respect to distinct
types of interactions allowing us to get unprecedented insights about the
underlying social organization and group dynamics. At similar spatial scales,
another emblematic example, accounting for the socio-spatial interdependence
typical of many other networks, concerns the organization of ecological
systems (Fig. 56): for instance, it has been recently shown that the
importance of dispersers for an ecosystem is better captured by multilayer
measures of importance rather than by standard metrics
[timoteo2018multilayer].

Figure 56:  Aggregate (top) and multilayer (bottom) representations of a real
seed-dispersal network. Reproduced from [timoteo2018multilayer] under Creative
Commons Attribution 4.0 International License.

At the human scale, relatively small-scale social systems
($10^{4}\leavevmode\nobreak\ $m) have been studied to better understand the
impact of external shocks on social structure and dynamics. For instance, it
has been shown that multilayer modeling can disentangle the different social
ties that conform to the proximity interaction networks of extant hunter-
gatherer societies, identifying social relationships with a key role in the
spread and accumulation of culture [migliano2017characterization]. In another
study, the analysis of mixed economies in three villages of Alaska unraveled
that factors related to climate change, such as global warming, have a non-
negligible effect on household structure, but the most important factors for
vulnerability are indeed due to social shift rather than resource depletion
[baggio2016multiplex].

Let us jump by three more orders of magnitude and discuss systems at the
planetary scale: the ones based on information exchange thanks to large-scale
communication infrastructures, such as the Internet. At this scale the
activity of a complex system might be very frenetic; think about an online
social media platform, where millions or billions of users worldwide
continuously produce content to be shared, which quickly travels and bounces
from one country to another. For instance, it is interesting to study how
rumors and memes spread on these systems, as recently proposed in
[d2019spreading] to capture the behavior of users who post information from
one social media platform to another and to provide a plausible explanation
for the heavy-tailed distributions of meme popularity that is usually observed
in empirical data.

Recent applications also include the analysis of trade networks and their
nested and modular structure [torreggiani2018identifying, a2018unfolding,
alves2019nested]. Additionally, as anticipated, multilayer models of financial
networks have been proposed very recently, to better explain how the financial
distress of one country can ignite a global financial crisis, like the one in
2008. By considering distinct asset types as layers where countries (nodes)
exchange financial assets (links), [del2020multiplex] have highlighted the
importance of both intra-layer and inter-layer connectivity contributions to
the propagation of contagions (Fig. 57).

Figure 57: Multilayer representation of a global financial network. Layers are
asset types, nodes are countries and links are cross-country financial
relationships. Figure from [del2020multiplex].

The long journey summarized in the last few pages of this work allowed us to
stress the importance of multilayer modeling across a broad spectrum of
disciplines, including cell biology, neuroscience, ecology and social
sciences, spanning 10 orders of magnitude from a spatial scale of about
$10^{-6}$ m to one of about $10^{4}$ m. The theoretical and computational
framework presented in this work provides researchers and practitioners with a
versatile and unified tool kit to shed light and gain new insights on the
complexity of natural and artificial systems.

##  Appendix A Master Stability Function (MSF) formalism

To test the stability of the synchronized state s we study how the
perturbation error $\delta\textbf{x}_{i}(t)=\textbf{x}_{i}(t)-\textbf{s}(t)$
evolves. By assuming small perturbations, we can write the variational
equation for $\delta\textbf{x}_{i}$:

$\delta\textbf{x}_{i}=J\textbf{F}(\textbf{s})\delta\textbf{x}_{i}-\sigma
J\textbf{H}(\textbf{s})\sum_{j=1}^{N}L_{ij}\delta\textbf{x}_{j}$ (105)

where $J\textbf{F}$ and $J\textbf{H}$ are the Jacobian of F and H,
respectively.

To find the solution of Eq. (105) we have to project $\delta\textbf{x}$ into
the eigenspace formed by the eigenvectors of the Laplacian matrix L, obtaining
a decomposition of the time evolution of perturbation error into $N$ decoupled
eigenmodes:

$\dot{\boldsymbol{\xi}_{i}}=[J\textbf{F}(\textbf{s})-\sigma\lambda_{i}J\textbf{H}(\textbf{s})]\boldsymbol{\xi}_{i},\quad
i=1,...,N$ (106)

where $\boldsymbol{\xi}_{i}$ is the eigenmode associated with the eigenvalue
$\lambda_{i}$ of L. By indicating as $\Lambda_{max}$ the maximum Lyapunov
exponent associated with the system of equations (106), then we can write the
time evolution of $\boldsymbol{\xi}$ as $|\boldsymbol{\xi}|\sim
e^{\Lambda_{max}t}$ and, finally, we found that a necessary condition for the
stability of a synchronized state is that $\Lambda_{max}<0$. The expression of
$\Lambda_{max}$ as a function of a generic parameter
$\alpha_{i}=\sigma\lambda_{i}$ is named the Master Stability Function (MSF),
and it usually assumes negative values in an interval
$\alpha\in(\alpha_{1},\alpha_{2})$. It means that, for a fixed coupling
strength $\sigma$, a network can reach and maintain complete synchronization
only if its structure, defined by its Laplacian matrix, is such that:

$\alpha_{1}<\sigma\lambda_{2}\leq\sigma\lambda_{\leq}...\leq\sigma\lambda_{N}<\alpha_{2}.$
(107)

or, equivalently,

$R\equiv\frac{\lambda_{N}}{\lambda_{2}}<\frac{\alpha_{2}}{\alpha_{1}},$ (108)

used, for instance, in [sole2013spectral] to unravel the existence of an
optimal value for the synchronizability of a multilayer system.

##  Appendix B Kuramoto model on networks

The first approach to model collective synchronization considers a population
of coupled limit-cycle oscillators whose natural frequencies are drawn from
some prescribed distribution and exert a phase-dependent influence on each
others. In formulae, we can write the Kuramoto model (KM) [Strogatz2000,
Arenas2008] as a system of $N$ oscillators, whose instantaneous phases
$\theta_{i}$ are described by the equation:

$\dot{\theta_{i}}=\omega_{i}+\frac{\sigma}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i}),$
(109)

where $\sigma$ is the coupling constant and $\omega_{i}$ is the natural
frequency of oscillator $i$, chosen from an unimodal distribution $g(\omega)$.

We can use a order parameter that describes the transition: it is defined as a
macroscopic measure that quantifies the collective rhythm produced by the
whole population:

$r(t)e^{i\Phi(t)}=\frac{1}{N}\sum_{j=1}^{N}\sin(\theta_{j}),$ (110)

where $\Phi(t)$ is the average phase, and $0\leq r(t)\leq 1$ measures the
phase coherence, where the two extreme values correspond to phase locked
($r=1)$ or incoherent oscillators. Equation (109) is then rewritten in terms
of the order parameter, and the equation for the instantaneous phase reduces
to

$\dot{\theta_{i}}=\omega_{i}+\sigma r\sin(\theta_{i}),$ (111)

which has two types of long-term behaviours. Oscillators for which
$|\omega_{i}|\leq\sigma r$ are phase-locked and form a mutually synchronized
cluster. Oscillators with frequencies in the tails of $g(\omega)$
distribution, where $|\omega_{i}|>\sigma r$ holds, are drifting with respect
to the synchronized cluster.

The Kuramoto model is then generalized on networks by including in Eq. (109)
information about network connectivity:

$\dot{\theta_{i}}=\omega_{i}+\sum_{j=1}^{N}\sigma_{ij}A_{ij}\sin(\theta_{j}-\theta_{i}).$
(112)

##  Appendix C Transitions in multilayer systems

Type of phase transition |  Mode |  Reference
---|---|---
Enhanced diffusion |  C |  [gomez2013diffusion]
Emergence of multiplexity |  D |  [radicchi2013abrupt]

Table 2:  Phase transitions (C: continuous, D: discontinuous, H: hybrid) in
multilayer networks due to the algebraic, topological and simple dynamics
described in Sec. IV.1.

Type of dynamics |  Mode |  Reference
---|---|---
Synchronization |  C |  [Gambuzza2015]
Explosive synchronization with local order parameter |  D |  [Zhang2015]
Synchronization of oscillators with random walks |  C & D |  [nicosia2017collective]
Interacting diseases |  C |  [Sanz2014]
Epidemic onset control by raising of awareness |  C |  [Granell2013]
Explosive pandemics with no early-warning |  D |  [danziger2019dynamic]
Cooperative behaviour in coupled competitive games – social influence |  C |  [amato2017interplay]
Cooperative behaviour in Prisoner’s Dilemma game |  C |  [gomez2012evolution]
Cooperative behaviour in Public Good game |  C |  [battiston2017determinants]

Table 3:  Phase transitions in multilayer networks due to the simple and
couple dynamics described in Sec. IV.2, Sec. IV.3 and Sec. IV.4.

Type of network |  Mode |  Reference
---|---|---
Multilayer network |  C |  [leicht2009percolation]
Multiplex |  D |  [son2012percolation]
Partial Multiplex |  C & D |  [son2012percolation]
Multiplex with overlap |  C, D & H |  [hu2013percolation, baxter2016correlated, cellai2016message]
Multiplex with spatially embedded networks |  C & D |  [son2011percolation, bashan2013extreme, grassberger2015percolation]
Bond Perc. on Multiplex |  C |  [hackett2016bond]
Multiplex & Local cascades |  D |  [buldyrev2010catastrophic]
Multiplex & Non-local cascades |  D |  [artime2020abrupt]

Table 4:  Type of transitions for different scenarios in percolation and
cascades propagation.