{'fraction_non_alphanumeric': 0.03562937062937063, 'fraction_numerical': 0.011975524475524475, 'mean_word_length': 4.284157043879907, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interacting opinion and disease dynamics in multiplex networks: Discontinuous phase transition and nonmonotonic consensus times

Opinion formation and disease spreading are among the most studied dynamical processes on complex networks. In real societies, it is expected that these two processes depend on and affect each other. However, little is known about the effects of opinion dynamics over disease dynamics and vice versa, since most studies treat them separately. In this work we study the dynamics of the voter model for opinion formation intertwined with that of the contact process for disease spreading, in a population of agents that interact via two types of connections, social and contact. These two interacting dynamics take place on two layers of networks, coupled through a fraction q of links present in both networks. The probability that an agent updates its state depends on both the opinion and disease states of the interacting partner. We find that the opinion dynamics has striking consequences on the statistical properties of disease spreading. The most important is that the smooth (continuous) transition from a healthy to an endemic phase observed in the contact process, as the infection probability increases beyond a threshold, becomes abrupt (discontinuous) in the two-layer system. Therefore, disregarding the effects of social dynamics on epidemics propagation may lead to a misestimation of the real magnitude of the spreading. Also, an endemic-healthy discontinuous transition is found when the coupling q overcomes a threshold value. Furthermore, we show that the disease dynamics delays the opinion consensus, leading to a consensus time that varies nonmonotonically with q in a large range of the model's parameters. A mean-field approach reveals that the coupled dynamics of opinions and disease can be approximately described by the dynamics of the voter model decoupled from that of the contact process, with effective probabilities of opinion and disease transmission.


I. INTRODUCTION
The formation of opinions and the propagation of an epidemic disease on a population of individuals are among the most studied dynamical processes on complex networks [1,2]. The behavior of each of these two processes has been explored independently of one another for the last decades, and many of their propagation properties on diverse complex topologies are well established already (see [2] and [3] for recent reviews on opinion formation and epidemic spreading, respectively). However, less attention has been paid to a possible case scenario where the dynamics of opinions interact with that of the disease spreading. In fact, it is hardly expected that these two dynamics are isolated in real societies but rather depend on and affect each other, since they both run at the same time on the same population: an individual can transmit a disease to a colleague while having a conversation and exchanging ideas or opinions on a given topic. Then, the following questions arise: does the dynamics of opinion formation have an impact on the extent and prevalence of the epidemic? Does the disease spread facilitate the ultimate dominance of one opinion, or does it rather hinder the consensus of opinions?
In an attempt to explore these questions, we study in this article how opinion formation and disease spreading processes affect each other, using two simple models as a proxy of each process: the voter model (VM) and the contact process (CP). The VM was originally introduced as the simplest system of interacting particles that can be exactly solvable in any dimension [4][5][6], and is one of the most studied models for opinion consensus. In this model, individuals can take one of two possible positions or opinions on a given issue, and are allowed to update them by adopting the opinion of a randomly chosen neighbor. The CP, on its part, has been extensively studied to explore the spread of an infection in a system of interacting agents [7], where infected agents can transmit the infection to susceptible neighbors in a lattice [8] or a complex network [9], and they can also recover at a given rate. The CP exhibits a continuous transition from a healthy to an endemic phase when the infection rate exceeds a threshold value. To model the interaction between the two dynamics we implement the framework of multilayer complex networks [10][11][12] that consists of a set of complex networks interrelated with one another, which allows to study systems composed by many interdependent processes. In the present study we consider that the opinion dynamics takes place on a network of social relations -formed by individuals that influence each other on a social issue, while the disease spreads on a network of physical contacts -formed by people having daily face-to-face contacts. All individuals are in both layers of networks, but the pattern of connections between them may be different in each layer. The overlap of connections is taken as a measure of the coupling between the two networks. The bilayer network system described above may represent a simple case scenario where the social network supports a process that involves peer pressure, like the adoption of new behaviors or opinions, while the contact network supports the spreading of a contagious viral infection like flu, which is transmitted by proximity or direct contact between individuals.
The different combinations of connections may correspond to different types of relationships between two individuals. For instance, two close friends can have both a contact and social tie, as they can see each other at work every day and also interchange ideas on a political issue. But it can also happen that individuals are connected by only one type of tie; e.g., two colleagues having a contact or proximity relation because they work in the same place but never talk about politics; or two friends that never meet but discuss political ideas by electronic means (phone, Facebook, Twitter, email, etc).
Some related works on multilayer networks [13][14][15][16][17][18][19] have also explored the interrelation between two information spreading processes. For instance, in references [13,14] the authors analyzed how the awareness of a disease affects the epidemic spreading on a multiplex network, by using the unaware-aware-unaware and the susceptible-infected-susceptible cyclic models, respectively. The interplay between opinion formation and decision making processes was studied in [15] using two interconnected networks. Another work considered two political parties (two interacting networks) that compete for votes in a political election [16]. In a recent article [17], the authors studied the dynamics of the voter model on bilayer networks with coevolving connections, while in [18] the same authors explored whether is appropriate to reduce the dynamics of the voter model from a two-layer multiplex network to a single layer. A recent work [20] considers a complex threshold dynamics that competes with a simple Susceptible-Infected-Susceptible dynamics on two interconnected networks.
All the works listed above explore the interplay between two social or two epidemiological processes that are alike. However, there is a lack of specific studies on the interplay between opinion and disease dynamics.
In this article we show that the dynamics of opinions has striking consequences on the disease spreading and vice versa. The nature of the healthy-endemic transition observed in the CP, as the infection probability increases, is largely modified by the dynamics of the VM. The transition changes from continuous to discontinuous when the disease and opinions are coupled, showing a jump in the disease prevalence at the transition point, where the magnitude of the jump increases with the coupling. Also, a discontinuous transition from an endemic to a healthy phase is found when the coupling overtakes a threshold value. In addition, we find that the dynamics of the CP has important consequences in the dynamical properties of the VM. The diffusion of opinions is slowed down by the disease in a nontrivial manner as the coupling increases. This leads to consensus times that vary either monotonically or non-monotonically with the coupling, for a large range of the model's parameters. We develop a mean-field approach to study the time evolution of macroscopic quantities, which takes into account state correlations between neighbors in the same network (pair approximation). This approach reveals that the interdependent system of opinions and disease can be thought of as two independent systems, with external parameters that depend on the coupling. Specifically, the opinion dynamics can be approximated as the dynamics of the VM on an isolated network, with an effective probability of opinion transmission that decreases with the coupling and the prevalence. Analogously, the disease spreading is approximately described by the CP dynamics on an isolated network, with an effective infection probability that decreases with the coupling and the fraction of neighbors with different opinions.
The article is organized as follows. In section II, we introduce the multiplex framework and the dynamics of the model on each layer. We present simulations results in section III and develop an analytical approach in section IV. Finally, in section V we give a summary and conclusions.

II. THE MODEL
We consider a bilayer system composed by a contact and a social network layer of mean degree k = µ and N nodes each. These two layers are interrelated through their nodes, which are the same in both networks, while links connecting nodes may not necessarily be the same. That is, both layers have the same number of nodes N and links µN/2, but the configuration of connections can be different in each layer. The overlap of links is measured by the fraction q (0 ≤ q ≤ 1) of links shared by both networks. In our model, the extreme values q = 0 and q = 1 correspond to totally uncoupled and totally coupled networks, respectively. To build this particular topology, we start by connecting the same pairs of nodes at random in both networks until the number of links reaches the overlap value q µN/2. Then, the rest of the links (1 − q)µN/2 are randomly placed between nodes in each network separately, making sure that the chosen pair of nodes in one network is not already linked in the other network.
Social links in this system connect individuals that influence each other on a given issue, while the infection is transmitted through contact links. In Fig. 1 we illustrate the bilayer system composed by a social and a contact network (top and middle layers), and its representation as a single layer with two types of links (bottom layer). We observe in Fig. 1 that nodes i and j are connected by both a social and a contact link, representing individuals that have a daily face-to-face conversation, where they interchange opinions and also one can infect its partner. Nodes j and k are only connected by a social link: they do not have face-to-face contacts but still exchange ideas electronically or by phone. Nodes i and k are only connected by a contact link: they have face-to-face or proximity contacts but they do not discuss and interchange opinions about the given issue.
To mimic the spreading of opinions and the disease we use the voter model (VM) and the contact process (CP) on each layer, respectively. Each node is endowed with an opinion state O that can take two possible values O = +, − (see top layer of Fig. 1), and a disease state D = 0, 1 that represents the susceptible and infected states of an individual, respectively (middle layer of Fig. 1). These two dynamics are coupled through the opinion and disease states of nodes, which affect each other by reducing the flow of information between neighbors, as we describe below with a simple example.
Let's consider a situation in which two individuals have a daily social and physical contact because they see each other at work and talk about politics. On the one hand, we assume that each individual is less influenced by its partner when she/he is sick, because the sick partner normally stays at home or at hospitals, reducing physical contacts between them.
This makes social relations (and the interchange of opinions) less likely when at least one of the two individuals is sick. Thus, we consider that a social relation takes place with probability 1.0 if both social neighbors are healthy, and with a reduced probability p o ≤ 1.0 when one or both are sick. In case they have a social contact but not a physical contact (they do not see each other but discuss ideas by electronic means), the disease state is not supposed to affect the probability of social interactions between them. Therefore, the social interaction probability is not reduced by the disease and, for simplicity, is set to 1.0 as in the case of healthy neighbors.
On the other hand, we consider that physical contacts (and therefore infections) between the two social and contact neighbors are more likely to happen when they share the same opinion. This is a consequence of a sociological mechanism called homophily [21][22][23], i. e., the tendency for individuals to interact with similar others. The effects of homophily in the propagation of cultural attributes in a society were studied by Robert Alxerod using an agent-based model [21], in which the probability that two neighboring agents interact is proportional to their cultural similarity (the number of shared attributes). Following this idea, we assume that the contact probability between the neighbors when they have the same opinion is higher than that when they have different opinions. Therefore, we set to 1.0 the contact probability of same-opinion neighbors and denote by p d ≤ 1.0 the contact probability between opposite-opinion neighbors. Once they have a physical contact the infection is transmitted with probability β, leading to effective infection probabilities β and β p d ≤ β in each respective case. In a situation where there is a contact but not a social connection between two neighbors (they see each other but they don't talk about politics), opinions are not expected to affect (neither increase or decrease) the contact probability.
Therefore, this can be considered as an intermediate situation respect to the two cases mentioned above, where the contact probability should be smaller than 1.0 but larger than p d , leading to an infection probability between β and β p d . However, for simplicity we assume that the contact probability in the absence of a social relation is the same as that in homophilic relations (1.0), and thus the infection probability takes the value β. This approximation and the one mentioned above for the social interaction probability have the advantage of reducing the number of free parameters, allowing for a deeper analysis of the model which already exhibits a very rich behavior as we shall see.
We now define the dynamics of the model according to the interaction properties discussed above. In a single time step ∆t = 1/N an opinion and a disease update attempt take place in each network, as we describe below (see Fig. 2). In case they hold opposite opinions the transmission happens with probability β p d ≤ β.
neighbors, which is reduced to β p d ≤ β only in the case they are attached by a social link and they share different opinions [see Fig. 2

III. NUMERICAL RESULTS
The CP and the VM are two of the most studied dynamical processes [6]. A relevant feature of the CP is the existence of a transition from a healthy phase to an endemic phase as the infection probability overcomes a threshold value β c . The healthy phase is static, as all nodes are susceptible and infection events cannot occur. The endemic phase is active, where each node undergoes an infected-susceptible-infected cycle and the total number of infected nodes fluctuates around a stationary value. The healthy-endemic transition is continuous, and the critical value β c depends on the topological properties of the network [9].
For its part, the VM has been extensively used to explore opinion consensus on different network topologies [24][25][26][27][28][29]. It was found that the diffusion properties of opinions depend on the heterogeneity of the network. This is reflected in the mean consensus time, which is proportional to the ratio µ 2 /µ 2 [28,29], where µ and µ 2 are the first and second moments of the network's degree distribution.
The behavior described above is particular of each model on single isolated networks.
In order to explore how the properties of these two processes are affected when they are coupled through a multiplex network, we run extensive Monte Carlo (MC) simulations of the model described in section II, using two Erdös-Rényi (ER) networks of mean degree k = µ = 10 each. Initially, each node in the system is infected with probability 1/2, and adopts either opinion state + or − with equal probability 1/2. That is, the system starts from a symmetric initial condition with roughly 1/4 of nodes in each of the four possible opinion-infection states: + 0 , + 1 , − 0 and − 1 . In the next two subsections we study separately the effects of one dynamics over the other.

A. Effects of opinion formation on disease prevalence
We start the analysis of the model by describing the results related to the effects of opinion formation on the properties of disease spreading. In Fig. 3 we show the stationary fraction of infected nodes averaged over many independent realizations of the dynamics, ρ stat on q, where it suddenly decays to a value close to zero. The sudden decrease in ρ stat 1 becomes more abrupt as N increases, leading to a discontinuous change of ρ stat 1 at β c q in the thermodynamic limit (N → ∞). This behavior is reminiscent of a discontinuous transition.
We also see that the jump in ρ stat 1 decreases with q and vanishes for the uncoupled case q = 0, where the transition becomes continuous, in agreement with the known behavior of the CP on isolated networks. The critical point β c 0 ≃ 0.53 for q = 0 agrees very well with the one found in previous numerical and analytical works [9].
These results show that the dynamics of opinions has a profound effect on the statistical properties of disease spreading, changing the type of phase transition in the CP from a continuous transition in the absence of coupling (when the two dynamics are independent) to a discontinuous transition when the dynamics are coupled.
In order to achieve a deeper understanding of the nature of this transition we studied the time evolution of the fraction of infected nodes ρ 1 (t) for the case q = 0.4, where the transition point is β c 0.4 ≃ 0.58 (see Fig. 3). Solid lines in Fig. 4 correspond to results for networks of size N = 10 4 . As we can see, for β > β c 0.4 ≃ 0.58 the average value of ρ 1 (t) over many  for β = 0.59). We note that this peculiar non-monotonic temporal behavior is known to induce discontinuous transitions in social models with multiple states and constrains, like the Axelrod model (see for instance [23,30]).
As we explain below, the origin of the non-monotonic behavior of ρ stat with N = 10 4 nodes (solid lines). We observe that curves for N = 10 6 decay monotonically with time to a stationary value denoted by ρ stat 1,∞ (only shown for β = 0.6), which agrees with the minimum of the non-monotonic curves for N = 10 4 . We need to note that these states are not truly stationary, in the sense that ρ 1 exhibits a very long plateau (outside the shown scale) but eventually increases and reaches the same stationary value ρ stat 1 of the curves for N = 10 4 . We have checked that the length of the plateau diverges with N, and thus is infinitely large when N = ∞. Therefore, we take ρ stat 1,∞ as the stationary value when N = ∞. In Fig. 5 we observe that the numerical values ρ stat 1,∞ (open symbols) agree reasonable well with the theoretical approximation from Eq. (11) (dashed curves) for the three values of q, even though the agreement worsens as q gets larger. We also see that ρ stat 1,∞ decays continuously as β decreases and becomes zero at the same value β c q of the transition in the thermodynamic limit corresponding to ρ stat 1 (filled symbols). That is, the healthy-endemic transition is continuous in an infinite system.
Up to here we studied the response of the system when the infection probability is varied, for a fixed coupling. We now explore the effects of having a varying coupling on disease prevalence. In Fig. 6  In this section we explore how the spreading of the disease affects the dynamics of opinions. As the transmission of opinions between neighboring nodes is more difficult when at least one of them is sick, we are particularly interested in studying up to what extent the disease slows down opinion diffusion over the social network, and how that depends on q, β, p o and p d . A way to quantify this is by looking at the time to reach opinion consensus.
In Fig. 8  We observe in Fig. 8 that the q-dependence of τ is quite diverse, showing monotonic as well as non-monotonic behaviors. This is a consequence of the competition between two different mechanisms that directly affect opinion transmission. One is the link overlap between the two networks that is proportional to q, and the other is the disease prevalence that decreases with q, as we explain below. The opinion transmission through a social link that overlaps with a contact link is slowed down when at least one of the two nodes is infected and p o < 1. Therefore, the overall delay in opinion transmission caused by the total overlap tends to increase with q, and so does τ . This effect explains the initial monotonic increase of τ as q increases from 0, in all curves. However, as q becomes larger a second effect becomes important: the fraction of infected nodes decreases with q [see inset of Fig. 8(a)], due to the coupling with the opinion dynamics that reduces the effective infection probability as discussed in section III A. Then, lower disease prevalence translates into fewer social links affected by the disease and, therefore, into a smaller opinion delay. This effect tends to reduce τ with q.
With these two mechanisms at play, the shapes of curves in Fig. 8 for different values becoming dominant for q above 0.8 when τ decays, and leading to a non-monotonic behavior of τ (q). Finally, for p 0 = 0.0 we observe a non-monotonicity similar than that of the p d = 0.5 curve, but with the addition that τ becomes very similar to τ 0 for all values of q > 0.6. This is because ρ stat 1 becomes zero above q ≃ 0.583 and thus the disease has no effect on opinions, leading to consensus times similar to the ones measured in isolated networks.
These In Fig. 9 we plot the normalized mean consensus time τ /τ 0 as a function of the infection probability β obtained from Eq. (17). Panels (a) and (b) correspond to couplings q = 0.5 and q = 1.0, respectively. To analyze these plots we recall that, as explained above, consensus times increase with the level of disease prevalence in the contact network, given that a larger disease prevalence translates into a larger delay in opinion propagation and in the subsequent consensus. A first simple observation is that τ increases with β and also with p d , as we expect from the fact that a larger value of β and p d implies a larger disease prevalence.
A second observation is that τ decreases with the likelihood of opinion transmission p o , as explained before when we compared τ in Fig. 8(a) with Fig. 8(b). A third observation is that τ approaches a value independent of p d when β goes to 1.0. This is because for β = 1.0 (recovery probability equals zero) and a fixed value of q and p o all nodes are infected at the stationary state, independently on the value of p d , and thus consensus times are the same for all p d . As we see in Fig. 9(b), the case q = 1 and p o = 0 is special because τ diverges as β approaches 1.0. This happens because in this situation the transmission of opinions is only possible between connected nodes that are both susceptible, which vanish in the β → 1.0 limit, leading to divergent consensus times. A rough estimation of how τ scales with β can be obtained by assuming that τ is proportional to the time scale associated to the opinion transmission across two given neighboring nodes in the social network, i and j, with opinions + and −, respectively. As q = 1, i and j are also neighbors in the contact network. Starting from a situation where i and j are infected for high β, the opinion transmission happens after both nodes recover. Therefore, τ is determined by the time it takes the 1-1 contact link to become a 0-0 link, which scales as (1 − β) −2 . In section IV we derived a more accurate expression for τ that exhibits this quadratic divergence in the β → 1 limit, shown in the inset of Fig. 9(b) by solid lines.

IV. ANALYTICAL APPROACH
In order to gain an insight into the behavior of the two-layer system described in section III, we develop here a MF approach that allows to study the time evolution and analogously for the contact layer In appendices A, B, C and D we develop a mean-field approach that allows to obtain the following system of coupled differential equations for ρ + , ρ +− , ρ 1 and ρ 10 , respectively: with and with These equations represent an approximate mathematical description of the time evolution of the model on infinitely large networks, where finite-size fluctuations are neglected. We note that Eqs. (3) and Eqs. (5) are coupled through the prefactors ω and γ, which depend on the coupling q and describe the opinion and disease dynamics, respectively. The interested reader can find in the appendices the details of the derivation of these equations. For the sake of simplicity, we assumed in the derivation that all nodes have the same number of neighbors k = µ chosen at random, which is equivalent to assuming that networks are degree-regular random graphs. However, we expect this approximation to work well in networks with homogeneous degree distributions like the ER networks we used in the MC simulations. We also implemented an homogeneous pair approximation [29] that takes into account correlations between the state of neighboring nodes within the same layer (intralayer pair approximation), but neglects correlations between opinion and disease states of both layers (interlayer annealing approximation). That is, we considered that the opinion state of each node is uncorrelated with its own disease state and with its neighbors' disease states and, conversely, that its disease state is uncorrelated with its own and its neighbors' opinion states.
It is instructive to analyze the structure of Eqs. (3) and (5) It seems that the disease dynamics is not able to break the intrinsic symmetry of opinion states induced by the voter dynamics. Equation (3b) for the evolution of ρ +− has an extra prefactor ω compared to the corresponding equation for the VM on isolated networks [29], which reveals that the disease affects the dynamics of opinions through its prevalence level, expressed by ρ 1 and ρ 10 [see Eq. (4)]. As discussed in appendix A, ω can be interpreted as the "effective probability" that a node i adopts the opinion of a chosen social neighbor j with opposite opinion, which depends on the disease state of both i and j. Within a MF approach, we can assume that the probability that i copies j's opinion depends on the disease state of an "average pair" of contact neighbors, and that this probability is the same for all social neighbors. In these terms, ω becomes the average copying probability over Thus, our MF approach assumes that this "effective infection probability" from i to j depends on the opinion states of an "average pair" of neighbors on the social layer, and that is the same for all contact neighbors. We can say that, at the MF level, the disease dynamics follows the standard CP on a single isolated network with homogeneous infection probability γβ and recovery probability 1 − β in each node. Therefore, the dynamics of opinions has an effect on the disease dynamics equivalent to that of an external homogeneous field acting on each node of the contact network, reducing the probability of infection between neighbors by a factor γ, while keeping the same recovery probability.
In the next two subsections we derive analytical expressions for the disease prevalence ρ stat 1 and the mean consensus time τ , from the system of Eqs. (3)(4)(5)(6).

A. Disease prevalence
In order to study how the opinion dynamics affects the disease prevalence, we find the fraction of infected nodes at the stationary state ρ stat 1 from Eqs. (3) and (5). We start by setting the four time derivatives to zero, substituting ρ 10 by 2(1 − β)ρ 1 /γβ from Eq. (5a) into Eq. (5b), and solving for ρ 1 . After doing some algebra we obtain two solutions, but only one is stable depending on the values of the parameters. The non-trivial solution corresponds to the endemic phase, where a fraction of nodes is infected, and is stable only where we used ρ + = ρ + (0) given that ρ + remains constant over time, as mentioned before.
We notice that ω does not affect the stationary value of ρ +− , which remains the same as in the original VM [29]. For a symmetric initial condition on the social layer (ρ + (0) = 1/2), as the one used in the simulations, we have ρ +− stat = (µ − 2)/[2(µ − 1)]. Replacing this last expression for ρ +− stat in Eq. (6) we obtain the following expression for γ: Finally, plugging Eq. (9) into Eq. (7) we arrive to the following approximate expression for the stationary fraction of infected nodes in the endemic phase: which is plotted in Figs. 5 and 6 (dashed curves). As the MF theory is meant to work for infinitely large systems, we also plot for comparison the numerical results obtained from simulations for very large networks (open symbols). We observe that, in all cases, the estimated theoretical value of the fraction of infected nodes from Eq. (11) is larger than that from simulations. As we explain below, this due to the fact that correlations between opinion and disease states are neglected by the MF approach. We first notice that an infection event 0 → 1 between two neighbors connected by a social and a contact link is only possible when the states of nodes are + 1 and + 0 (a ++ 1 0 pair) or − 1 and − 0 ( a −− 1 0 pair), because p d = 0 in Figs. 5 and 6. Then, it is expected that ++ social links are negatively correlated with 10 contact links and positively correlated with 11 and 00 contact links, given that sameopinion neighbors tend to infect each other and thus, at a given time, they are more likely to be either both infected or both susceptible. However, the theoretical approximation assumes that ++ social links are uncorrelated with 10 contact links (see appendix A) and, therefore, the estimated probability of finding a ++ 1 0 pair is larger than that obtained when negative correlations are considered. The same conclusion also holds for −− 1 0 pairs. This leads to a theoretical overestimation of the number of ++ 1 0 and −− 1 0 pairs and, consequently, to a larger rate of infections which increases the disease prevalence respect to numerical results, as we see in Figs. 5 and 6. Figure 5 shows that ρ stat 1 form Eq. (11) continuously decreases and vanishes as β decreases beyond a threshold value, as it happens in the standard CP. This shows that the transition to the healthy state is continuous within the MF approach, which assumes that the system is infinitely large. In Fig. 6 we see that ρ stat 1 decreases with q, reducing the prevalence and inducing a transition to the healthy phase. That is, Eq. (10) predicts a healthy-endemic continuous transition as β and q are varied, which happens at the point where ρ stat 1 vanishes, leading to the relation The transition line obtained from Eq. (12) for µ = 10 is plotted in Fig. 7 for p d = 0 (dashed curve). We can see that the agreement with simulations is good for small values of the coupling q, but discrepancies arise as q increases, where the theoretical prediction from Eq. (13) overestimates numerical values. Another simple observation that follows from Eq. (13) is that for we obtain the nonphysical value q c > 1. This means that, in the network model, it is possible to induce a transition by increasing the coupling only when β is lower than a given value, as we see in Fig. 7 for β < 0.68.
As a final remark we stress that the transitions within this MF approach are continuous, in agreement with simulations in very large networks. This is so because Eqs. (3) and (5) correspond to an infinite system where finite-size fluctuations are neglected.

B. Opinion consensus times
In this section we study the quantitative effects of the disease on the time to reach opinion consensus. For that, we find an analytical estimation of the mean consensus time τ as a function of the model parameters.
As mentioned in section IV, in infinitely large systems ρ + remains constant over time [see Eq. (3a)]. However, in finite systems ρ + fluctuates until it reaches either value ρ + = 1 (+ consensus) or ρ + = 0 (− consensus), with both configurations characterized by the absence of +− social links (ρ +− = 0). A typical evolution of ρ +− towards the absorbing state can be seen in Fig. 4 (b) for q = 0.4 on networks with N = 10 4 nodes. That is, consensus is eventually achieved in finite systems due to the stochastic nature of the opinion dynamics, which leads the social network to a state where all nodes share the same opinion. In a single opinion update ρ + may increase or decrease by 1/N with the same probability ωρ +− /2, calculated as the probability ρ +− /2 that a node and an opposite-opinion neighbor are selected at random, times the probability ω of opinion adoption. Therefore, the stochastic dynamics of the VM can be studied by mapping ρ + into the position of a symmetric one-dimensional random walker on the interval [0, 1], with a jumping probability proportional to ωρ +− stat /2 and a step length of 1/N. Starting from a symmetric configuration with N/2 nodes with + opinion (ρ + (0) = 1/2), the walker reaches either absorbing point ρ + = 1 or ρ + = 0 in an average number of steps that scales as N 2 . Then, given that the walker makes a single step in an average number of attempts that scales as 1/ωρ +− stat , and that the time increases by 1/N in each attempt, we find that the mean consensus time scales as As we see in Eq. (8), ρ +− stat is independent of the disease prevalence ρ 1 and, therefore, the prevalence affects τ only through the effective copying probability ω, which sets the time scale associated to opinion updates. From Eq. (4) we see that ω equals 1.0 when the layers are uncoupled (q = 0) or when p o = 1.0, and thus the dynamics of opinions is exactly the same as that of the original VM. However, ω is smaller than 1.0 in the presence of coupling (q > 0) and p o < 1.0, and thus the evolution of the dynamics is "slowed down" -in averageby a factor 1/ω > 1.0, given that opinions are copied at a rate that is ω times smaller than in the uncoupled case. As a consequence, τ increases by a factor 1/ω respect to the mean consensus time in the uncoupled case To obtain a complete expression for the ratio τ /τ 0 as a function of the model's parameters we express ω in terms of ρ stat 1 , by substituting into Eq. (4) the stationary value of ρ 10 that follows from Eq. (5a), ρ stat 10 = 2(1 − β)ρ stat 1 /γβ. This leads to In the healthy phase is ρ stat 1 = 0, thus ω = 1.0 and τ = τ 0 . In this case, the theory predicts that the disease has no effect on the time to consensus because there are no infected nodes that can affect the opinion dynamics. However, having a value ρ stat 1 > 0 of infected nodes in the endemic phase has the effect of reducing ω or, equivalently, increasing τ respect to τ 0 .
Plugging into Eq. (16) the expressions for γ and ρ stat 1 from Eqs. (9) and (10), respectively, and reordering some terms, we obtain the following expression that relates τ and τ 0 in the endemic phase: In Fig. 8 we plot in solid lines the ratio τ /τ 0 vs q from Eq. (17) for µ = 10.
We observe that the theoretical values of τ /τ 0 are smaller than those obtained from numerical simulations (symbols) for all combinations of p o and p d shown. A possible explanation of these discrepancies can be given by analyzing how correlations affect the estimated number of different types of connected nodes, as we have done in section IV A for disease prevalence. If we take the p 0 = 0 case, we see that an opinion change due to an interaction between two social and contact neighbors happens only if both nodes are susceptible, that is, when they have states + 0 and − 0 . Then, as the theory assumes that +− social links and 00 contact links are uncorrelated (appendix A), the estimated probability of finding a +− 0 0 pair is larger than that in simulations, given that +− social links are expected to be negatively correlated with 00 contact links. This negative correlation is due to the fact that susceptible neighbors tend to align their opinions and, therefore, they are more likely to be in the same opinion state at a given time. This leads to an overestimation of the number of +− 0 0 pairs and, therefore, to a larger rate of opinion transmission. This has the overall effect of speeding up consensus, decreasing the theoretically estimated mean time to reach consensus respect to the mean consensus time measured in simulations, as we see in Fig. 8.
Even though discrepancies with numerical results increase with the coupling q, the analytic expression (17)  is a transition to the healthy phase when q overcomes a value q c < 1 given by Eq. (12) and, therefore, τ = τ 0 for all q > q c [see p d = 0 curves in the main plot and the inset of Fig. 8(a)].
As a consequence, τ /τ 0 exhibits a non-monotonic behavior with q, as we described in section III B. For higher values of p d , the transition to the healthy phase does not happen for the physical values q ≤ 1 used in the model's simulations, given that q c > 1 from Eq. (12).
In this case, τ may either increase monotonically with q for large p d values (see p d = 1.0 curves), or have a maximum at some intermediate value for medium p d values (see p d = 0.5 curves). As explained in section III B, the non-monotonicity is a consequence of the competition between the level of link overlap among the two layers -which increases with q-and the disease prevalence -which decreases with q. This competition can be seen quantitatively in Eq. (16) for ω, which has three factors that depend on q and affect τ . Besides the factor proportional to q, the factor 1/γ also increases with q, as seen from Eq. (9). But these two factors are balanced by ρ stat 1 , which decreases with q. An interesting case is the one for full coupling q = 1.0 and p o = 0, because τ from Eq. (17) diverges as β approaches 1.0. This happens in the model because when β = 1.0 once a node becomes infected it remains infected forever. Then, once all nodes become infected the opinion dynamics stops, as infected neighboring nodes cannot interchange opinions, and thus the social layer freezes in a mixed state of + and − opinions and consensus is never achieved. By doing a Taylor series expansion of expression (17) up to to second order in the small parameter ǫ = 1 − β ≪ 1 we obtain, after some algebra, where we used µ = 10. Equation (18) shows that τ diverges as (1 − β) −2 in the β → 1.0 limit, as shown in the inset of Fig. 9(b). For β = 1.0 and p o = 0, we can check from Eq. (17) that τ /τ 0 ≃ 1/(1 − q), which shows the divergence of τ as the system approaches the fully coupled state q = 1.0.

V. SUMMARY AND CONCLUSIONS
We proposed a bilayer network model to explore the interplay between the dynamics of opinion formation and disease spreading in a population of individuals. We used the voter model and the contact process to simulate the opinion and the disease dynamics running on a social and contact network, respectively. These two networks share the same nodes and they are coupled by a fraction q of links in common. We showed that, when the networks are coupled, the opinion dynamics can dramatically change the statistical properties of the disease spreading, which in turn modifies the properties of the propagation of opinions, as compared to the case of isolated networks.
The VM dynamics is able to change the order of the healthy-endemic phase transition observed in the CP as the infection probability β exceeds a threshold value β c , from a continuous transition for the uncoupled case to a discontinuous transition when the coupling q is larger than zero. The magnitude of the change in the disease prevalence at the transition point β c increases with q. The discontinuity is associated with the non-monotonic time evolution of the fraction of infected nodes. This non-monotonicity is as a consequence of the time-varying nature of the effective infection probability, which varies over time according to the stochastic evolution of the fraction of +− social links. The system also exhibits a discontinuous transition from an endemic to a healthy phase when the coupling overcomes a value q c , for a fixed value of β. The origin of this discontinuity is the same as that of the discontinuous transition with β, that is, the non-monotonicity in the time evolution of the fraction of infected nodes. We also obtained a phase diagram in the β − q space showing the healthy and endemic phases for different values of the probabilities p d and p o . In all cases, we observed that the transition point β c increases with q.
We need to mention that changes in the order of topological and dynamical transitions were already observed in multilayer networks [31][32][33][34][35][36][37]. We developed a mean-field approach that allowed to estimate with reasonable precision the healthy-endemic transition line (β c , q c ) as a function of the model's parameters. This approach reveals that the disease dynamics is equivalent to that of the standard CP on an isolated network, with an effective infection probability that is constant over time and that decreases with the coupling and the stationary fraction of +− social links, for a fixed value of β. This means that, at the mean-field level, the overall effect of the VM on the CP is to decrease the effective infection probability as the coupling increases. Therefore, as q increases, a larger value of β is needed to bring the system to the endemic phase, leading to an increase of the transition point β c with q.
On its part, the CP dynamics has the overall effect of slowing down the propagation of opinions, delaying the process of opinion consensus compared to the one observed in an isolated network. The MF approach reveals that the opinion dynamics corresponds to that of the standard VM model on an isolated network, with a probability of opinion transmission that decreases with q and the disease prevalence. Depending on the parameters values, the mean consensus time τ can show a monotonic increase with q, as well as a non-monotonic behavior. An insight on these results was given by the MF approach, which allowed to obtain an approximate mathematical expression that relates τ with the parameters. This approach shows that the behavior of τ with q is the result of two different mechanisms at play: the overlap of social and contact links that tends to increase τ with q, which is counterbalanced by the fraction of infected nodes that tends to decrease τ with q. Therefore, the non-trivial dependence of τ with q is a consequence of the competition between these two mechanisms.
It is interesting to note that, despite the nontrivial interplay between the CP and the VM, the coupled interdependent system of opinions and disease can be roughly seen as two systems that evolve independently of one another, where each system has an effective parameter that depends on the other dynamics and the coupling. Specifically, the opinion dynamics corresponds to that of the VM with an effective opinion transmission probability that decreases with the disease prevalence and the coupling, while the disease spreading is well described by the dynamics of the CP with an effective infection probability that decreases with the fraction of +− social neighbors and the coupling. However, this is only an approximation that comes from the MF analysis, which neglects correlations between opinion and disease states.
The results presented in this article correspond to a particular initial state that consists on even fractions of + and − opinion states and even fractions of infected and susceptible states, uniformly distributed over the networks. As a future work, it might be worth studying the behavior of the system under different initial conditions, and with uneven fractions of opinion and disease states. For example, one can simulate a population with initial polarized opinions based on the disease, by correlating the opinion of each node with its disease state (for instance by infecting all nodes with − opinion and leaving all + opinion nodes in the healthy state). Finally, it would be interesting to study the behavior of the present model under different update rules. For instance, we have checked a simple rule in which the connection condition -connected or disconnected-between two nodes in one layer is not taken into account for the update in the other layer. This is an ongoing work with some preliminary results that suggest that the critical behavior of this new model is quite different from that of the original model.

ACKNOWLEDGMENTS
We thank Gabriel Baglietto and Didier Vega-Oliveros for helpful comments on the manuscript. We acknowledge financial support from CONICET (PIP 0443/2014).
We denote by O D the state of a given node, where O = +, − and D = 1, 0 are its opinion and disease states, respectively. Thus, there are four possible node states: + 0 , + 1 , − 0 and − 1 . In a single time step of the dynamics, the transitions from state − D to state + D when a node switches opinion from − to + lead to a gain of 1/N in ρ + , while the transitions + D → − D when there is a − → + opinion change lead to a loss of 1/N in ρ + . Considering these four possible transitions, the average change of ρ + in a single time step of time interval ∆t = 1/N is described by the rate equation where for instance the term ∆ρ + and similarly for the − 1 → + 1 transition corresponding to the second term in Eq. (A1) Third and fourth terms in Eq. (A1) are obtained by interchanging symbols + and − in Eqs. (A2) and (A3), respectively, due to the symmetry between + and − opinion states.
We notice that disease states remain the same after the interactions, as only a change in the social layer can lead to a change in ρ + . The first term in Eq. (A2) represents the average change in ρ + due to interactions in which a node i in state − 0 copies the opinion of one its social neighbors j in state + 0 , changing the state of i to + 0 . This interaction is schematically represented by the symbol −+ 0 0 , where the horizontal line over the opinion symbols describes a social link between i and j. In the same way, the symbol −+ 0 0 represents an interaction between a − 0 node and a neighboring + 0 node connected by both a social and a contact link that are indicated by horizontal lines on top of the respective symbols.
The second, third and fourth terms in Eq. (A2) describe, respectively, the transitions due to an interaction of node i with a + 0 social/contact neighbor, a + 1 social neighbor and a + 1 social/contact neighbor. We now illustrate how to build an approximate expression for each term of Eq. (A2) for . For the sake of simplicity, we assume that all nodes have the same number of neighbors k = µ chosen at random, which is equivalent to assuming that networks are degreeregular random graphs. However, we expect this approximation to work well in networks with homogeneous degree distributions like ER networks. The first term in Eq. (A2) can be written as which can be understood as the product of the different probabilities associated to each of the consecutive events that lead to the −+ when O=+,− where the prefactor p o in the last term accounts for the probability of copying the opinion of an infected contact neighbor. Keeping in mind that we aim to obtain a closed system of rate equations for ρ + , ρ 1 , ρ +− and ρ 10 , we now find approximate expressions for the different probabilities of Eq. (A6) in terms of the fractions of nodes and links in each layer. We start by assuming that correlations between opinion and disease states of a given node are negligible, and thus we can write Then, to estimate the conditional probabilities P [ −+ 0 D ] and P [ −+ 0 D ] it proves convenient to split each of them into two conditional probabilities using the relation P (a, b|c) = P (a|c)P (b|a, c) and interpreting the entire event of connecting a given type of link to a [ + D ] node as two separate events. Assuming that the type of link connected to node i is uncorrelated with the state of i, we have Within an homogeneous pair approximation [29], the probability P (+| −) that a social neighbor j of a node i with opinion O i = − has opinion O j = + can be estimated as the ratio between the total number µNρ +− /2 of links from − to + nodes and the total number µNρ − of links connected to − nodes, that is P (+| −) ≃ ρ +− /2ρ − . Similarly, we estimate the probability that a contact neighbor j of a susceptible node has disease state D j = 0 as P (0| 0) ≃ ρ 00 /ρ 0 , and disease state D j = 1 as P (1| 0) ≃ ρ 10 /2ρ 0 . And if j is not a neighbor of i on the contact layer then P (D| 0) ≃ ρ D . Assembling all these factors we obtain from Eqs. (A8) and for P − 0 from Eq. (A7) we arrive to where we have used the conservation relations Eqs. (2a) and (2d).
We now calculate the second gain term in Eq. (A1), ∆ρ + − 1 → + 1 , which represents the average change in ρ + due to − 1 → + 1 transitions, following the same steps as above for the term ∆ρ and using the approximations for the conditional probabilities we arrive to where we have used the conservation relations Eqs. (2a) and (2c).
Adding Eqs. (A9) and (A11) we obtain the following expression for the average gain of a + node in single time step, corresponding to the sum of the first and second terms of with The prefactor ω plays an important role in the dynamics of opinion consensus, by setting the time scale associated to opinion updates, and can be interpreted as an effective probability that a node adopts the opinion of randomly chosen opposite-opinion neighbor. That is, Eq. (A12) for the gain of a + node simply describes the process of selecting a − node i and a + neighbor j, which happens with probability ρ +− /2, and then switching i's opinion with a probability ω that depends on the connection type and disease state of both i and j. This "effective copying probability" ω turns out to be an average copying probability over the entire social network, as shown in section IV.
In order to find the equation for the average loss of a + node in a time step, corresponding to the sum of the third and forth terms of Eq. (A1), we can exploit the symmetry between + and − opinion states and simply interchange signs + and − in Eq. (A12) where we used ρ −+ = ρ +− . Finally, adding Eqs. (A12) and (A14) we obtain quoted in Eq. (3a) of the main text. Therefore, the fractions of + and − nodes are conserved at all times: ρ + (t) = ρ + (0) and ρ − (t) = ρ − (0) = 1 − ρ + (0). Even though the above calculation leads to a very simple result, it serves as an introduction to the methodology used for deriving rate equations for the other fractions ρ +− , ρ 1 and ρ 10 , as we show next.
Appendix B: Derivation of the rate equation for ρ +− In analogy to the calculation for ρ + in the previous section, the average change of the faction of +− social links ρ +− in a time step is given by the rate equation with , which takes into account the specific configuration of links and neighbors connected to i, as depicted in Fig. 10. We obtain where the first and second moments of M N − 0 , µ are Here D i = 1, 0 and D j = 1, 0 are the disease states of nodes i and j, respectively. Replacing the expressions for the moments from Eqs. (B7) in Eq. (B6) and regrouping terms we obtain Plugging the expressions for the probabilities P [ −+ 0 D ] and P [ −+ 0 D ] from Eq. (A8) into Eq. (B8), and after doing some algebra we finally obtain We now follow an approach similar to the one above for ∆ρ +− By adding Eqs. (B9) and (B11) we obtain the following expression for the change in ρ +− due to − → + transitions Then, by interchanging sings + and − in Eq. (B12) we obtain the change in ρ + due to Finally, adding Eqs. (B12) and (B13) we arrive to the following rate equation for ρ +− quoted in Eq. (3b) of the main text The average change of the fraction of infected nodes ρ 1 in a single time step can be written as where each term represents a different transition corresponding to a disease update on the contact layer. The first term of Eq. (C1) corresponds to the recovery of a + 1 node, and can be estimated as That is, with probability P + 1 ≃ ρ + ρ 1 a + 1 node is picked at random, and then recovers with probability 1 − β, decreasing ρ 1 in 1/N. The second term corresponds to the infection of a + 0 node, while the last two terms are equivalent to the first two, but where a node with opinion − is recovered and infected, respectively. By the symmetry of + and − opinions, the last two terms are obtained by interchanging signs + and − in the first two.
We now find an approximate expression for the second term of Eq. (C1). A susceptible node j in state + 0 can be infected by a sick neighbor i with + or − opinion and connected to j by a contact link or by both a social and a contact link. Thus, four possible contact interactions lead to the + 0 → + 1 transition: The symbol O + 1 0 represents a contact interaction between node i in state O 1 (O = +, −) and node j in state + 0 . The state that changes in the interaction is now displayed on the right-hand side of the symbol, instead on the left-hand side as for the case of the social interactions described in the previous sections. This is because the chosen neighbor j of i changes state in the CP, while in the VM is node i who changes state. Taking into account the events and their associated probabilities that lead to each of the four interactions described above, we can write Eq. (C3) as The first and third terms of Eq. (C4) correspond to selecting an O 1 node i and a contact neighbor j with state + 0 at random, which happens with probability P O 1 N [ O+ 1 0 ]/µ, and then i infecting j with probability β, given that they are not connected by a social link. The second and fourth terms are similar to the first and second terms, respectively, but selecting a social/contact neighbor j. As both types of links are present in this case, i infects j we finally arrive to where we have used the conservation relations from Eqs. (1a) and (1c). Now that we estimated the first two terms of Eq. (C1), the last two terms are obtained by interchanging signs + and − in Eqs. (C2) and (C8): Adding Eqs. (C2), (C8), (C9) and (C10), the rate equation (C1) for ρ 1 becomes as quoted in Eqs. (5a) and (6) where the expression in square brackets is the change in the number of 10 links connected to a node i in state + 1 when i recovers, given a specific configuration of node types connected to i [see Fig. 11 Replacing these expressions for the probabilities and using the conservation relations from Eqs. (1a) and (1c) we obtain, after doing some algebra, We now calculate the second term of Eq. (D1) corresponding to the change in ρ 10 after the infection of a node with + opinion: The term called P in Eq. (D5) -the two summations inside the square brackets-is the probability that an O 1 node i infects a + 0 neighbor j, and is the same as the one calculated in Eq. (C4) for ∆ρ 1 + 0 → + 1 , which is estimated in Eq. (C8) as We notice that the extra 1/N prefactor in Eq. (C8) comes from the change in ρ 1 , which for ρ 10 depends on the neighborhood of node j. The subindex i in the term P indicates that the infection probability term depends only on node i and its neighborhood [see Fig. 11(b)].
The term called C corresponding to the summation outside the square brackets expresses the change in ρ 10 when node j gets infected [see Fig. 11(b)]. Here the subindex j refers to node j and its neighborhood. This term carries the information that the infection on j comes from one of its infected neighbors i, and thus it is known already that at least one of j's neighbors has disease state D i = 1. This is taken into account by running the summation on the other µ − 1 unknown neighbors and considering that the number of 10 links connected to j is at least one, which is added to the total number of 10 links inside the parentheses.