PHYSICAL REVIEW E 99, 062311 (2019)
Explosive transitions induced by interdependent contagion-consensus dynamics
in multiplex networks
D. Soriano-Paños,
1,2
Q. Guo,
3,4,*
V. L a t o r a ,
5,6,7,†
and J. Gómez-Gardeñes
1,2,‡
1
GOTHAM Laboratory, Institute for Biocomputation and Physics of Complex Systems (BIFI), University of Zaragoza, 50018 Zaragoza, Spain
2
Departamento de Física de la Materia Condensada, Universidad de Zaragoza, 50009 Zaragoza, Spain
3
School of Mathematics and Systems Science, Beihang University & Key Laboratory of Mathematics Informatics Behavioral Semantics
(LMIB), Beijing 100191, China
4
China Construction Bank, Beijing 100033, China
5
School of Mathematical Sciences, Queen Mary University of London, London E1 4NS, United Kingdom
6
Dipartimento di Fisica ed Astronomia, Università di Catania and INFN, Catania I-95123, Italy
7
Complexity Science Hub Vienna (CSHV), 1080 Vienna, Austria
(Received 20 December 2018; revised manuscript received 22 April 2019; published 25 June 2019)
We introduce a model to study the interplay between information spreading and opinion formation in social
systems. Our framework consists in a two-layer multiplex network where opinion dynamics takes place in
one layer, while information spreads on the other one. The two dynamical processes are mutually coupled
in such a way that the control parameters governing the dynamics of the node states at one layer depend on
the dynamical states at the other layer. In particular, we consider the case in which consensus is favored by
the common adoption of information, while information spreading is boosted between agents sharing similar
opinions. Numerical simulations of the model point out that, when the coupling between the dynamics of the two
layers is strong enough, a double explosive transition, i.e., a discontinuous transition both in consensus dynamics
and in information spreading appears. Such explosive transitions lead to bi-stability regions in which the
consensus-informed states and the disagreement-uninformed states are both stable solutions of the intertwined
dynamics.
DOI: 10.1103/PhysRevE.99.062311
I. INTRODUCTION
The functioning of a wide range of complex systems in
physics, biology, and social sciences is subject to collective
phenomena such as the onset of synchronization [1], the
emergence of norms and cooperation [2], or the diffusion
of epidemics [3], among others. In the past two decades,
a number of works have analyzed the role played by the
structure of the networks governing the interactions among
the constituents of a complex system in the emerging of its
collective dynamics [4–8]. Our understanding of the funda-
mental mechanisms driving these phenomena is of utmost im-
portance, as it provides a solid basis for modeling, predicting,
and controlling real dynamical systems [9–11].
Recently, complexity and network science have moved
one step forward in this direction by considering that, very
frequently, the elements of many real complex systems are
subject to different types of interactions at the same time.
In some cases, the dynamical processes that occur simulta-
neously in a system depend on each others. Examples of the
coexistence and non-trivial interdependence of two or more
dynamical processes are very common in social systems and
in the natural sciences. For instance, human prevention be-
*
quantongg@buaa.edu.cn
†
v.latora@qmul.ac.uk
‡
gardenes@unizar.es
haviors coexist and co-evolve with disease spreading [12–15],
while different dynamical processes interact in neural systems
and govern the structure-function relationships in the human
brain [16].
The study of such coupled dynamical processes has been
largely stimulated by the introduction of novel frameworks to
deal with networks with many layers or networks of networks
[17–20]. Multiplex networks are indeed the natural way to
model the existence of different dynamical interactions among
the same set of units [21–31].
In this work we introduce and study a model of two
coevolving socially inspired processes: formation of opinions
[32] and information spreading [33]. In our model the two
dynamics are mutually coupled in such a way that the trans-
mission of information from a spreader to a receiver is boosted
when the neighbors of the latter share similar opinions [34].
In addition to this, the alignment of the opinion of an agent
to those of her neighbors is fostered when such neighbors
spread the voice simultaneously. The model, and in particular
the adopted type of interdependence of the two processes,
captures everyday life examples in which the use of technol-
ogy or the adoption of new ideas by an individual happens
in virtue of the consensus found among her acquaintances
[35,36] and, in turn, the common adoption of these novelties
boosts the degree of homophily needed for the creation of
social consensus [37,38]. The study of the model reveals how
the interplay between opinion and spreading dynamics can
dramatically alter the critical properties of the two dynamical
2470-0045/2019/99(6)/062311(8) 062311-1 ©2019 American Physical Society
D. SORIANO-PAÑOS et al. PHYSICAL REVIEW E 99, 062311 (2019)
processes, leading to abrupt onsets of epidemics and con-
sensus. These explosive onsets are discontinuous transitions
leading to the appearance of bistable regions where the multi-
plex network can switch between active to inactive dynamical
phases triggered by small perturbations. It is remarkable that
the explosive transitions that we observe in our model are not
produced by any of the standard mechanisms at work in single
networks [39], but they are the result of the coupling induced
by the multiplex architecture.
The work is organized as follows. In Sec. II we introduce
the model and discuss the rationale under the adopted dynam-
ical coupling between layers. Then, in Sec. III, we illustrate,
through numerical simulations of the model, the emergence
of explosive transitions both in the consensus and in the
contagion dynamics. In Sec. IV, we numerically explore the
role that both topological and dynamical features of the model
play on the emergence of these abrupt transitions. Finally, in
Sec. V we summarize the main results of this work and we
discuss possible future directions.
II. MODEL OF INTERDEPENDENT DYNAMICS
To describe the delicate interplay between information
spreading and the formation of consensus in a social system,
we introduce here a model in which the two processes take
place at the two layers of a multiplex network with M = 2
layers and are mutually coupled. We deal with a multiplex
network following the assumption that there exists a one-
to-one correspondence between nodes (the social agents) in
different layers, so that each layer is composed by the same
set of N nodes. However, the topologies of the two layers can
in general be different and are described by the adjacency
matrices A
[1]
={a
[1]
ij
} and A
[2]
={a
[2]
ij
}, respectively. These
matrices are defined such that a
[1]
ij
= 1(a
[2]
ij
= 1) if a link
exists between nodes i and j in the first (second) layer,
while a
[1]
ij
= 0(a
[2]
ij
= 0) otherwise. We denote the degree of
node i in the first (second) layer as k
[1]
i
=
N
j=1
a
[1]
ij
(k
[2]
i
=
N
j=1
a
[2]
ij
).
Our model can be phrased in terms of a general
formalism for interdependent dynamical networks
proposed in Ref. [29]. If we denote, respectively,
as x(t) ={x
1
(t ), x
2
(t ),...,x
N
(t )}∈
N
and y(t) =
{y
1
(t ), y
2
(t ),...,y
N
(t )}∈
N
the states of the nodes at
the two layers, then the evolution of the system can be written
as
˙x
i
= F
ξ
i
(x, A
[1]
)
˙y
i
= G
η
i
(y, A
[2]
) i = 1, 2,...N, (1)
where the dynamics of state x
i
(y
i
) of node i in the first
(second) layer is governed by a function F
ξ
(G
η
)ofthe
dynamical state x (y) and of the structure A
[1]
(A
[2]
)ofthe
first (second) layer. Notice that, following Ref. [29], functions
F
ξ
and G
η
in our model are taken to be dependent on the
parameters ξ and η, and this is the key ingredient to connect
the two dynamical processes. Namely, we assume that the
parameter ξ
i
of function F
ξ
i
at the first layer is itself a function
of time depending on the dynamical states {y
j
(t )} at the
second layer of the neighbors j of node i at the first layer
FIG. 1. Left: Schematic representation of our model on a mul-
tiplex network with M = 2 layers and N = 5 nodes. The first (top)
layer accounts for the consensus dynamics, which is modeled by a
Kuramoto model as in Eq. (3), whereas the second (bottom) layer
describes the spreading of information according to the SIS model
as in Eq. (4). Right: The coupling strength λ between opinions (top)
as well as the contagion rate β (bottom) have been modified as in
Eqs. (7)and(6), respectively, to mutually couple the synchronization
process to the spreading of information.
(a
[1]
ij
= 1). Analogously, the evolution of the parameter η
i
at
the second layer depends on the states {x
j
(t )} at the first layer
of the neighbors of node i at the second layer (a
[2]
ij
= 1). In
this way, the system of Eq. (1) is completed by the following
system of equations:
ξ
i
(t ) = f
y
j
(t )|a
[1]
ij
= 1
η
i
(t ) = g
x
j
(t )|a
[2]
ij
= 1
i = 1, 2,...N, (2)
where f and g are two assigned functions.
As illustrated in Fig. 1, the first layer in our model accounts
for the dynamics underlying the formation of consensus in
a social system, while the second layer describes the conta-
gion processes mimicking the spread of ideas/products. The
dynamical state x
i
(t ) of node i at the first layer represents
the opinion of individual i, described as a phase variable,
i.e., x
i
(t ) = θ
i
(t ) ∈ [−π,π]. The time evolution of x
i
(t )is
modeled via the Kuramoto model of coupled phase-oscillators
[40–42], so that the first set of equations in the system of
Eq. (1) reads
˙
θ
i
(t ) = F
ξ
i
(θ,A
[1]
) = ω
i
+ λξ
i
(t )
N
j=1
a
[1]
ij
sin[θ
j
(t ) − θ
i
(t )],
(3)
where ω
i
is the natural frequency of node i. Notice that λ is
a global coupling strength, while the local coupling strength
associated to node i is modulated by the dynamical variable
ξ
i
(t ) that changes in time depending on the dynamics of the
second layer, as sketched in Eq. (2), in a way that will be
specified below.
The dynamical state y
i
(t ) of node i at the second layer
represents the probability of node i of being active as
user/spreader of an idea, namely, y
i
(t ) = p
i
(t ) ∈ [0, 1]. The
062311-2
EXPLOSIVE TRANSITIONS INDUCED BY … PHYSICAL REVIEW E 99, 062311 (2019)
time evolution of p
i
(t ) is modeled through a susceptible-
infected-susceptible (SIS) model. In this way we identify
susceptible (ignorant or uninformed) agents as those who
do not have and transmit information, whereas the infected
ones correspond to active users (informed and spreaders) who
disseminate the information to the rest of the population.
Under this framework, a susceptible node (a node in state
S) with an infected neighbor can be infected by it at time t
through the process S + I → 2I and becomes itself a spreader
(state I) with a probability βη
i
(t ). In addition, a spreader can
return to its ignorant state through the process I → S with a
probability μ. Such a dynamics can be cast in the form of
a Markov evolution for the probability p
i
(t ) that a node i is
spreader at time t as [43–45]
˙p
i
(t ) =−μp
i
(t ) + [1 − p
i
(t )]
×
⎧
⎨
⎩
1 −
N
j=1
1 − a
[2]
ij
βη
i
(t )p
j
(t )
⎫
⎬
⎭
. (4)
Notice that, at variance with the usual SIS model, here the
microscopic contagion probability βη
i
(t ) of node i may differ
from node to node, and also change in time due to the presence
of the time dependence in factor η
i
(t ), in close analogy with
the presence of the time-dependent factor ξ
i
(t ) in the effective
coupling of unit i at the consensus layer.
To completely define our model we finally need to assign
the time-dependent functions {ξ
i
(t )} and {η
i
(t )} that mutually
couple the consensus dynamics and the process of contagion
as sketched in Eq. (2). To define η
i
(t ), we need to capture
the influence that consensus at layer 1 has on the contagion
dynamics at layer 2. With this purpose we evaluate the local
degree of consensus r
i
(t ) around node i at time t by con-
sidering the values of θ
j
(t ) in the neighborhood of node i.
Notice, however, that the neighbors of node i are taken in
the second layer, where information spreading takes place.
This is because it is the consensus among potential spreaders
that facilitates the transmission of ideas. We therefore use the
adjacency matrix {a
[2]
ij
} to construct the neighborhood of node
i. The local degree of consensus of node i is defined as the
modulus of the complex function:
r
i
(t )e
iψ
i
(t )
=
1
k
[2]
i
N
j=1
a
[2]
ij
e
iθ
j
(t )
, (5)
so that we get r
i
0 in the absence of local consensus and
r
i
= 1, otherwise. Once evaluated r
i
(t ), we can write the
second of Eq. (2)as
η
i
(t ) =
1
1 + exp{−α[r
i
(t ) − r
∗
]}
. (6)
The use of the Fermi function with a tuning parameter α>0
implies that, for large enough values of α, when r
i
(t ) → 0,
i.e., when the local consensus around i is small, the conta-
gion probability toward i, βη
i
(t ), tends to 0. However, when
consensus among the neighbors of i increases, their influence
over i also grows, approaching β as r
i
(t ) → 1. In this way the
value r
∗
acts as a threshold, so that for r
i
(t ) > r
∗
(r
i
(t ) < r
∗
)
we have η
i
(t ) > 0.5(η
i
(t ) < 0.5). For the sake of simplicity,
in the following we fix r
∗
= 0.5.
Last, we model the influence that the contagion dynamics
of layer 2 has on the formation of consensus at layer 1. To
this aim, the node-depending coupling constant λξ
i
(t )ofthe
Kuramoto model at layer 1 is chosen to be dependent on the
number of spreaders around node i at layer 1. Specifically,
ξ
i
(t ) is defined as the fraction of spreaders among the neigh-
bors of node i in layer 1, so that the first of Eqs. (2) reads
ξ
i
(t ) =
N
j=1
a
[1]
ij
p
j
(t )
k
[1]
i
. (7)
Notice that this time we have made use of the adjacency
matrix {a
[1]
ij
} to construct the neighborhood of node i .
Summing up, in our model the state [θ
i
(t ), p
i
(t )] of each
node i, with i = 1, 2,...,N, evolves in time as in Eqs. (3)
and (4), where the two parameters ξ
i
and η
i
depend in turn
on the state [θ
i
(t ), p
i
(t )] as in Eqs. (7) and (6), mutually
coupling the two dynamical processes. Notice that in this
way both the infection probability βη
i
and the Kuramoto
coupling strength λξ
i
of a node i are obtained by taking
average over the neighbors in the layer that governs the
corresponding dynamics, i.e., layer 1 for ξ
i
and layer 2 for
η
i
. However, the averaged dynamical quantities correspond
to the node states at the other layer, i.e., the phases for η
i
and the probabilities of being infected for ξ
i
, thus closing the
feedback loop between spreading and consensus dynamics.
The way the interdependence between these two processes
has been modelled follows, as discussed above, the rationale
that the existence of consensus facilitates the adoption of ideas
and that it is the simultaneous spread of ideas that fosters the
alignment of opinions.
III. RESULTS
To characterize the effects of the interplay between spread-
ing and consensus dynamics, we explore the dynamical be-
havior of our coevolving model mainly focusing on the onset
of synchronization and on the appearance of an endemic state.
To this aim, we start by infecting a small fraction ρ of agents
and by initially setting the oscillator phases θ
i
at random
within a range θ
i
∈ (−π,π]. The natural frequencies of oscil-
lators {w
i
} are also randomly chosen within w
i
∈ [−0.5, 0.5].
The particular values of these individual properties of nodes
remain the same for all the numerical experiments. This way,
we avoid the stochastic noise inherent to a random assignment
of the initial conditions. We take λ and β as the natural control
parameters for Kuramoto and SIS dynamics, respectively. The
order parameters are also the usual ones for both dynamical
systems. Namely, the degree of global consensus is measured
by using the Kuramoto order parameter r defined by the
complex number:
r(t )e
iψ (t )
=
1
N
N
j=1
e
iθ
j
(t )
, (8)
which represents the centroid of all oscillators placed on the
complex unit circle. In its turn, for the SIS model we monitor
062311-3
D. SORIANO-PAÑOS et al. PHYSICAL REVIEW E 99, 062311 (2019)
FIG. 2. Average global consensus r (top) and fraction of infor-
mation spreaders I (bottom) as a function of the coupling constant λ
for the SF-ER multiplex configuration. These order parameters have
been computed adiabatically by increasing the value of the coupling
constant λ (Forward) from λ = 0 or by decreasing it (Backward)
from λ = 1.5. The contagion rate values used are for (a) and (b) β =
0.70 and for (c) and (d) β = 1. The rest of the model parameters are
set to α = 10 and μ = 1.0. All the layers are composed by N = 500
nodes with average degree k=4. The SF network follows a power
law distribution with exponent γ = 3.
the evolution of the fraction of infected individuals:
I(t ) =
1
N
N
j=1
p
j
(t ). (9)
As usual, the order parameters, r and I, are measured by
making a time average of r(t ) and I(t ), once the stationary
regime of the dynamics is reached. To reach this stationary
state, we integrate Eq. (3) by using the fourth-order Runge-
Kutta method and Eq. (4) using an Euler method, both with
time steps δt = 0.01.
The networks used to build the multiplex configurations are
random Erdös-Renyi (ER) and scale-free (SF) networks with
N = 500 nodes and average degree k4. In particular, the
SF networks are constructed according to the Barabási-Albert
method [46], so the degree distribution follows a power law
with exponent γ = 3. The use of these two topologies allows
us to study the role of degree heterogeneity in the evolution of
consensus and spreading dynamics. As anticipated above, we
denote the multiplex considering that the first layer contains
consensus dynamics whereas information spreading takes
place on top of the second one.
In Fig. 2, we have computed the diagrams for global
consensus and fraction of infected people using a SF-ER
multiplex by keeping fixed the contagion probability, β,in
the ER layer and varying the consensus coupling in the SF
one. To this aim, we have computed the forward (increasing
λ) and backward (decreasing λ) diagrams. Figures 2(a), 2(b)
FIG. 3. Diagrams of the average global consensus r (top) and
of the fraction of spreaders nodes I (bottom) as a function of
the infectivity β and the coupling parameter λ for SF-ER (left)
and SF-SF (right) multiplex networks. The color code encodes the
magnitudes of the two order parameters. The striped regions in
the panels highlight the parameter region (λ, β), where hysteresis
cycles appear due to the coexistence of two stable solutions: total
consensus-disagreement in the synchronization layer and active-
inactive spreaders at the information spreading layer. The black solid
lines denote the critical values of the coupling λ separating the two
solutions. All the layers consists of N = 500 nodes with average
degree k=4. The SF networks follow a power law distribution
with exponent γ = 3.
and 2(c), 2(d) show drastically different transitions. On one
hand, Figs. 2(a) and 2(b), which correspond to β = 0.70,
show an abrupt transition both for the degree of consensus and
the fraction of spreaders. These diagrams are characterized
by the existence of regions of bistability where the solutions
corresponding to absence of global consensus and infor-
mation spreaders coexist with those displaying macroscopic
coherence and spreading. On the other hand, Figs. 2(c) and
2(d), corresponding to β = 1, show a smooth and continuous
transition, i.e., the expected onset from the usual Kuramoto
and SIS models. As we show below the particular type,
smooth of explosive, of transition depends on the multiplex
configuration and on the β value.
To have a broader picture about the phenomenon described
above, in Fig. 3 we represent the diagrams for global consen-
sus r (top row) and fraction of spreaders I (bottom row) as a
function of both β and λ for SF-ER (left panels) and SF-SF
(right panels) multiplexes. At first sight, for both topolo-
gies, we observe that below a critical value β
c
, represented
with dashed points, the single stable solution is the absence
of global consensus and information spreaders. Above this
threshold, we can find different stable solutions depending on
the Kuramoto coupling constant λ. Namely, for small values
of λ, the stable solution is the absence of global consensus and
the presence of a small fraction of spreaders I which depends
062311-4
EXPLOSIVE TRANSITIONS INDUCED BY … PHYSICAL REVIEW E 99, 062311 (2019)
on the value of the contagion rate β. This constitutes a
surprising result, since one would expect that, for large values
of α, the absence of global consensus prevents the diffusion of
information. However, the finite-size effects associated to the
reduced set of neighbors of each node in the spreading layer
leads to the emergence of spreaders. In particular, despite the
fact that there not exist a global shared opinion in the network
for the non-synchronized regime, each agent observes some
degree of local consensus which fosters the spread of ideas in
this model. In the next section, we deeply study this region
of the phase diagram unveiling the role that dynamical or
topological features of the intertwined dynamics play on the
diffusion of information despite the absence of consensus.
For larger values of λ, the type of transition observed
depends on the value of β. In particular, for β values close
to the diffusion threshold, the number of spreaders in the
non-synchronized regime is very small, so an abrupt transition
takes place toward the state of global consensus and the
existence of a macroscopic set of spreaders. Interestingly,
these abrupt transitions incorporate bi-stability regions (see
striped areas in Fig. 3) where the coexistence of two solutions
(corresponding to large and small order parameters) explains
the hysteresis cycles shown in Fig. 2. However, for large
values of β, a macroscopic set of spreaders already exists in
the non-synchronized regime, thus giving rise to a continuous
transition in which both consensus and number of spreaders
continuously grow while increasing the coupling constant λ.
At this point, we can understand the role that each process
plays on the intertwined dynamics. It becomes clear that the
epidemics behaves as the limiting process, for the emergence
of consensus requires the existence of active spreaders but no
viceversa. In its turn, the synchronization dynamics, moni-
tored by the coupling constant λ, behaves as an external force
which drives the system from a practically inactive phase to
an active one.
Once described the diagrams in the (β,α) plane let us
identify the main differences between SF-ER and SF-SF
multiplexes. By comparing panels corresponding to the SF-
ER configuration and those corresponding to the SF-SF one
in Fig. 3, it is clear that the value of the critical coupling
λ
c
separating the regions corresponding to the absence and
presence of global consensus is lower for the SF-SF config-
uration than for the SF-ER one. To explain this, we must
take into account that the degree distributions of both layers
in the SF-SF configuration are positively correlated, so that
hubs promote the interplay between consensus dynamics and
information spreading, thus anticipating the explosive onsets.
Another difference between both multiplexes, is that the
bistable regime is hindered in the SF-SF configuration with
respect to the SF-ER one.
IV. NUMERICAL STUDY OF THE INTERPLAY BETWEEN
TOPOLOGY AND DYNAMICS
In contrast to the differences discussed in the previous
section regarding the phase diagrams of SF-ER and SF-SF
multiplexes, the value of the diffusion threshold β
c
is roughly
the same in both cases. This is an unexpected result, since the
presence of hubs in heterogeneous networks boosts spread-
ing phenomena in the vast majority of dynamical models.
FIG. 4. Fraction of spreaders in the nonsynchronized regime
(λ = 0) as a function of the transmissibility β for both SF-ER
(solid lines) and SF-SF configurations (dashed lines). The line color
denotes the value of the social pressure α. Regarding the underlying
topologies, all layers are composed by N = 500 nodes with average
degree k=4. The degree distribution of the SF layer follows a
power law distribution with exponent γ = 3.
However, in our model of intertwined dynamics, the proba-
bility that highly connected nodes diffuse information is also
affected by the degree of consensus among their acquain-
tances. In this sense, the more neighbors agents have, the
wider is the set of opinions to which they can have access
in the nonsynchronized regime. As a consequence, highly
connected nodes experience smaller values of local consensus,
thus hindering the spreading dynamics. Therefore, there is a
competition between two opposite effects (the existence of
more spreading routes and the lack of local consensus) whose
outcome is governed by the value of α.
To shed light on this phenomenon, Fig. 4 reports the
number of spreaders in the nonsynchronized regime (λ = 0)
as a function of (β,α) for both SF-ER (solid lines) and
SF-SF (dashed lines) configurations. Interestingly, for large
values of α, a small degree of consensus around hubs is
enough to hinder their spreading ability, thus yielding a higher
diffusion threshold than for the SF-ER multiplex. However,
for small values of α, the lack of local consensus around
hubs becomes less relevant and the presence of highly con-
nected nodes promotes ideas spreading, thus anticipating the
diffusion threshold. This way, α can be interpreted as a kind
of social pressure shaping the transmissibility of ideas as a
function of their acceptance in society.
To further understand the role that α plays on the inter-
twined dynamics, we analyze in Fig. 5 the phase diagrams
of the spreading dynamics varying this parameter as well as
the average degree of the underlying SF-ER configuration.
Note that the shift of the diffusion threshold as a consequence
of the aforementioned interplay between topology and social
pressure becomes more evident in this case. In addition, Fig. 5
highlights the role that social pressure plays on the emergence
of abrupt or smooth transitions. Namely, increasing the social
pressure over agents turns the emergence of spreaders and full
062311-5
