Spreading Processes in Multiplex Metapopulations Containing Different Mobility Networks

We propose a theoretical framework for the study of spreading processes in structured metapopulations, with heterogeneous agents, subjected to different recurrent mobility patterns. We propose to represent the heterogeneity in the composition of the metapopulations as layers in a multiplex network, where nodes would correspond to geographical areas and layers account for the mobility patterns of agents of the same class. We analyze classical epidemic models within this framework and obtain an excellent agreement with extensive Monte Carlo simulations. This agreement allows us to derive analytical expressions of the epidemic threshold and to face the challenge of characterizing a real multiplex metapopulation, the city of Medellín in Colombia, where different recurrent mobility patterns are observed depending on the socioeconomic class of the agents. Our framework allows us to unveil the geographical location of those patches that trigger the epidemic state at the critical point. A careful exploration reveals that social mixing between classes and mobility crucially determines these critical patches and, more importantly, it can produce abrupt changes of the critical properties of the epidemic onset.


I. INTRODUCTION
During the last decades we have witnessed the onset of several major global health threats such as the 2003 spread of SARS, the H1N1 influenza pandemic in 2009, the western Africa 2014 Ebola outbreaks and more recently the Zika epidemics in the Americas and Caribbean regions.These outbreaks are increasingly characterized by the small elapsed time between initial infections in a single region to the global epidemic state affecting different cities, regions, countries and, in some cases, continents.Thus, in the recent years a great effort has been devoted to understand the fast unfold of emergent diseases and to design both local and global contention strategies.The most common avenue to tackle this problem is to adapt classical epidemic models taking into account the multiscale nature of diseases propagation [1,2].
It is clear that the spread of an emergent infectious disease is the result of human-human interactions in small geographical patches.However, in order to understand the geographical diffusion of diseases, one has to combine these microscopic contagion processes with the long-range disease propagation due to human mobility across different spatial scales.To tackle this problem, epidemic modeling has relied on reaction-diffusion dynamics in metapopulations, a family of models first used in the field of population ecology [3][4][5][6][7].For the case of epidemic modeling, the usual metapopulation scenario [8][9][10] is as follows.A population is distributed in a set of patches, being the size (number of individuals) of each patch in principle different.The individuals within each patch are well-mixed, i.e., pathogens can be transmitted from an infected host to any of the healthy agents placed in the same patch with the same probability.The second ingredient of metapopulation frameworks concerns the mobility of agents.
Each host is allowed to change its current location and occupy another patch, thus fostering the spread of pathogens at the system level.Mobility of agents between different patches is usually represented in terms of a network where nodes are locations while a link between two patches represents the possibility of moving between them.
The non-trivial mobility patterns observed in real populations [11] and the recent advances of network epidemiology [12] have motivated a thorough analysis about the impact that the structure of mobility networks has on the onset of globalscale contagions.In the last decade, important steps towards the inclusion of realistic mobility structures have been done [13][14][15][16].These approaches had to compromise between realism and analytical feasibility.On one side, lengthy mechanistic simulations [1,17] provide fair predictions on realistic scenarios while, on the other, theoretical frameworks allowing for analytical results usually rely on strong assumptions limiting their applicability to real-world threats.For instance, it is usual to assume simplified mobility patterns and meanfield approximations for hosts and patches behavior to be able to predict the onset of an outbreak.In this class of models, random diffusion of agents between the nodes is often used as proxy of human mobility while, as in the case of the heterogeneous mean field approach in scale-free networks [18], subpopulations with identical connectivity patterns are assumed to be equally affected by the disease.
These mean-field like approximations for patches having identical properties, while useful for deriving analytical results, add important limitations for their applicability in realworld diseases prediction.As data gathering techniques and epidemic surveillance [19] increase their accuracy, metapopulation models face new challenges [20].In an effort of relaxing the assumption about random diffusion of hosts and approaching realistic mobility patterns used in agent based simulations, recently the recurrent and spatially constrained nature of most human movements, such as daily commutes, has been addressed [21][22][23][24], at the expense of considering either meanfield assumptions or simple mobility networks.Therefore, a theoretical framework of general metapopulations of arbitrary structure, able at incorporating real mobility patterns, remains as an open theoretical challenge.
In addition, very recently, network science has tackled the formulation of networked systems in which different types of interactions between a given set of nodes coexist and interplay.These systems, termed as multiplex networks [25][26][27][28][29], consist on a set of L networks (usually called layers) and a set of N nodes.Each node is represented once in each network layer allowing it to share different connectivity patterns in each of the L layers.In terms of metapopulation models, for which nodes account for geographical locations, the multiplex formalism captures the coexistence, within each subpopulation, of different of agents with different mobility preferences (such as different species or social classes).This way, each layer represents the mobility network of each type of agent, while each subpopulation is represented in each layer.
In an attempt to increase the realism of epidemiological models without compromising the possibility of a theoretical analysis, here we propose a mathematical framework in which the dynamical variables of each patch forming the metapopulation are treated independently.Our framework can accommodate any mobility multiplex network from real commuting datasets containing different types of individuals and is amenable to any particular distribution of the population across the patches, then generalizing previous findings on monolayer networks [30,31].We will analyze the classical Susceptible-Infected-Susceptible (SIS) and the Susceptible-Infected-Recovered (SIR) models, achieving an excellent agreement with intensive Monte Carlo simulations.In addition, we derive an exact expression of the epidemic threshold and show its nontrivial dependence with the different mobility patterns represented in the multiplex.

II. METAPOPULATION MODEL
For the sake of clarity, we start considering a metapopulation framework consisting of one single type of agents.In this case we have a network composed of N nodes (the patches) and a total population of P agents.Importantly, each agent has associated one of the patches so that all the movements of agents associated to a node, say i, initiate from and return to it.In its turn, a node i of the network is the basement (or home) of a number n i of agents so that P = N i=1 n i .For the sake of generality, we consider that the network connecting the patches of the population is a weighted and directed graph, encoded in an adjacency matrix whose entries W ij account for the weight of the interaction from node i to node j.
The dynamical model implemented in our metapopulation involves three different stages at each time step t: movement, interaction and return (MIR).First, each agent decides whether moving with probability p or remaining in its associated home node i with probability (1 − p).If the agent moves, it goes to any of the nodes connected to i as dictated by the adjacency matrix W. The probability that a patch j is chosen, is proportional to the weight of the corresponding entry W ij of the adjacency matrix: Once all the agents have been placed in the nodes, the interaction stage takes place.Each agent updates its dynamical state according to the epidemic model at work (see below) by interacting with the agents that are placed in the same patch at time t.Finally, agents come back to their corresponding residence node and another time step starts.These stages are depicted in Fig. 1 where, for clarity, we have considered that the states of agents are either Healthy or Infectious as in the SIS model.This metapopulation model captures the commuting nature of most of human displacements within cities (at the level of neighborhoods) or countries (at the level of cities).Interestingly, let us remark that empirical data about real recurrent mobility patterns can be incorporated straightforward in the MIR model by considering the number of observed trips between two locations W ij in order to construct the transition rates matrix R.This way, the model has as control parameters the displacement probability p and those controlling the epidemic model under study.

III. POPULATION-BASED MARKOVIAN DYNAMICS IN COMPLEX NETWORKS
In the following we will focus on the two most paradigmatic epidemic models, SIS and SIR.The reaction laws of these models are given by two parameters: (i) the probability λ that a Susceptible (healthy) agent catches the diseases after the contact with a single infected individual and (ii) the probability µ that an infected overcomes the disease and turns to be susceptible again (SIS) or becomes immunized (SIR).These reactions can be expressed as: for the SIS model, and: for the SIR.
As in any metapopulation model on large complex networks, we face the problem of computationally expensive simulations.A useful avenue to analyze these models, with the byproduct of obtaining analytical estimations for the impact of the epidemic, is to formulate coarse-grained models that reduce significantly the complexity of the problem.Typically, heterogeneous mean field (HMF) techniques have been applied in a number of works related to epidemic spreading in contact networks and metapopulations.As anticipated above, the main assumption of HMF is to correlate the relevant parameters of nodes and patches with their number of connections to other nodes, i.e. their degree.This way, two distant patches that are connected to the same number (but not the same set) of locations are considered to have the same static and dynamical properties such as, for instance, the number of habitants and the fraction of infected agents.This assumption, although being strong, has been shown to be valid for small epidemic sizes, thus allowing quite good predictions of epidemic thresholds.
Here we formulate the mathematical equations of the MIR model by following a similar avenue as in [32][33][34] for contact networks, thus generalizing the Markovian approach to complex metapopulations.This way, we will consider both static and dynamical variables of each individual patch as independent, allowing us to compare directly with the findings of Monte Carlo simulations at the microscopic level and, more importantly, to derive theoretical results for any kind of particular mobility networks.

A. SIS model
For the SIS model, we have a set of N variables ρ i (t) denoting the fraction of infected agents associated to patch i at time t.It is important to stress that, according to the MIR model, an agent whose associated patch is i can be in other node j at time t.The time evolution of ρ i (t) can be written as: where the first term denotes the fraction of infected agents associated to i that do not recover at time t + 1.The second term instead accounts for the fraction of healthy agents associated to i that pass to infected at time t + 1.In this second term, Π i (t) is the probability that a healthy agent associated to node i becomes infected at time t.This probability reads: where the first term denotes the probability that a susceptible agent associated to patch i becomes infected when remaining at its home node i and the second one accounts for the probability that this agent catches the disease when moving to any neighbor of i.Finally, the probability P i (t) in Eq. ( 5) denotes the probability that a healthy agent in (but not necessarily associated to) node i at time t becomes infected after the contact with any of the infected agents present inside i at the same time.Then, probability P i (t) reads: where: being δ ij = 1 when i = j and δ ij = 0 otherwise.The expressions in Eqs. ( 4)-( 7) compose the closed set of equations covering the evolution of an SIS disease spreading in the MIR metapopulation model with parameters p, µ and λ.In addition, matrix R is given by the topology of the mobility network, that can be constructed from the observed flows between the patches, and the set of node populations, {n i }, can be also set according to the local census of the population under study.

B. SIR dynamics
The formulation of the Markovian equations for a metapopulation under a SIR spreading dynamics demands to add another set of N variables: {r i (t)} (i = 1, ..., N ), i.e., the fraction of recovered agents associated to patch i.Thus, the set of N equations ( 4) for the SIS model is now substituted by the following set of 2 • N equations: On the other hand, since the infection processes within each of the patches in the SIR model follow identical rules as those of the SIS one, the expression in Eq. ( 8) for the probability that a healthy agent associated to node i becomes infected at time t, Π i (t), has the same form as in the SIS case.This way, the SIR metapopulation dynamics is fully described by Eqs. ( 8) and ( 9) with the addition of Eqs. ( 5)- (7).

C. Validation of the Markovian equations
To check the accuracy of the Markovian equations we have considered Erdös-Rényi (ER) and Scale-free (SF) synthetic networks having the same number of nodes N = 10 3 and average connectivities k = 5.5 and k = 7.3 respectively.The nodes of these networks are homogeneously populated, n i = 5 • 10 3 ∀i, so that the total population of our systems is P = 5•10 6 .The weights W ij between the nodes of the graphs are randomly assigned following a homogeneous distribution within the range W ij ∈ [1,50].Once all the weights are set, we construct the transition matrix R [see Eq. ( 1)] for each graph.
Monte Carlo simulations start by infecting a small fraction of agents in each of the nodes.In particular, we infect each agent with probability 10 −3 so that, on average, there is 1 infected agent per node at time t = 0.This initial configuration corresponds to set as initial conditions of the Markovian equations ρ i (0) = 10 −3 ∀i (and r i (0) = 0 ∀i in the SIR case).For Monte Carlo simulations, due to the stochastic nature of the initial configuration and the disease models, we have averaged the results over 10 2 realizations for each combination of the parameters (p, λ and µ) considered.
First we analyze the SIS model.In Fig. 2 (top and bottom panels correspond to ER and SF networks respectively) we plot the number of infected agents in the steady state, I, as a function of the infection probability, λ, for different movement probabilities p.The points denote the results of Monte Carlo simulations for each value of λ and p while solid curves correspond to the solution of the Markovian equations.The agreement between simulations and the equations is almost exact, capturing with high accuracy the macroscopic state of the metapopulation both in the disease-free and epidemic regimes.For the SIR model in Fig. 3 (top and bottom panels correspond to ER and SF networks respectively) we plot the number of recovered agents, R, as function of λ for the same set of movement probabilities p as for the SIS model.Again, we observe the exact agreement between Monte Carlo simulations and the solution of Markovian equations both before and after the epidemic threshold.
The high accuracy of the solution of Markovian equations shown in Figs. 2 and 3 allows us to overcome the computational costs associated to large scale Monte Carlo simulations.However, it is clear that both I (SIS) and R (SIR) are macroscopic indicators of the outreach of the disease in the whole population.The examples shown above assume that the populations across nodes are homogeneously distributed.However, in real metapopulations, such as cities, each patch contains a different number of agents.These demographic heterogeneity may lead to interesting effects, such as the increase of the epidemic threshold with the increase of mobility [30,31].Here, due to the homogenous distribution of agent across patches, mobility always leads to a decrease of the epidemic onset (as shown in Figs. 2 and 3).
To check further the accuracy of the Markovian equations we now increase the level of resolution and monitor the spatiotemporal evolution of the infections in the population.To this aim we now fix λ, µ and p, and set a small infected fraction of agents placed at the same node.Then we run an SIR epidemics and follow how this infectious seed spreads across the patches of the metapopulation by monitoring r i (t).
The results of the agent based simulation are shown in the top panel in Fig. 4 where we have used an SF network of N = 200 nodes and a total population of P = 7 • 10 5 agents where each patch contains a number of agents proportional to its out-strength (s out i = N j=1 W ij ).We have ranked the nodes (from 1 to 200) according to the order in which infections occur, so that node 1 is the one in which the initial infectious seed is placed.As shown, after a transient time of roughly 50 time steps the system reaches the frozen state and the density of recovered agents at each patch remains stationary.Interestingly, the final pattern clearly shows that the local significance of the disease is not homogeneous.
The solution of the Markovian evolution equations (second panel of Fig. 4) fairly reproduces the spatio-temporal pattern from Monte Carlo simulations.To quantify the agreement, we have computed the error in quantifying the number of recovered and infected agents by the Markovian equations per time step as: where R i (t) and I i (t) are, respectively, the number of recovered and infected agents associated to node i at time t observed in Monte Carlo simulation while r i (t) and ρ i (t) are, respectively, the fractions of recovered and infected individuals as obtained when solving the Markovian equations.The bottom panel in Fig. 4 shows that the Error reaches its maximum value (8%) just after the avalanche of contagions across the metapopulation starts.After this peak, E(t) reaches a stationary value around 1% pointing out the high accuracy in reproducing the stationary pattern shown above.

IV. MULTIPLEX METAPOPULATIONS
We address now the case of metapopulations in which different types of agents interplay.In particular we will focus on systems in which agents displaying L types of mobility patterns coexist within each patch.This way, the population of a patch i is the sum of the number of agents of each type n i = L α=1 n α i and the probability that an agent of patch i and type α visits another patch j is now written as the generalization of Eq. ( 1): where W α ij is associated to the number of observed trips of agents of type α in patch i to patch j.
To analyze this situation, it is natural to make use of a multiplex formulation [25][26][27][28][29] of the metapopulation, as it is illustrated in Fig. 5.In our case, the number of layers of the multiplex is equal to the number of types of agents (L) and the architecture of each layer is described by a different matrix R α .Each patch of the system is represented as one node in each network layer and the corresponding L nodes are virtually connected (dotted lines) as they mix their agents when the contagion processes take place.
The number of Markovian equations of the multiplex are now multiplied by L with respect to the networked metapopulation.In particular for the SIS (and SIR) model, the variables are ρ α i (t) (and r α i (t)), which denote the fraction of infected (and recovered) individuals of layer α = 1, ..., L associated to node i.In this case, SIR equations become: while for the SIS model we only have Eq.( 12) with r α i (t) = 0.The term P α i (t), which denotes the probability that an agent of type α placed in patch i at time t becomes infected, reads: where λ βα is the probability that an diseased agent of type β infects a healthy agent of type α.In addition the number of agents of type α associated to patch j that travel to a different patch i is given by: The set of Eqs. ( 12)-( 15) conforms the Markovian model of the multiplex metapopulation.For the sake of simplicity, we will now restrict to the case λ αβ = λ ∀α, β, so that the infection probability between healthy and infected agents does not depend on their types.

A. Validation or the Markovian equations
To validate the Markovian equations for the multiplex metapopulation we proceed in the same fashion as we did for networked ones.First, we compute the impact that SIR and SIS diseases have as a function of the infectivity of the disease, λ, and the degree of mobility, p.We have studied three types of multiplex of L = 2 layers, namely ER-ER, SF-SF, and ER-SF, of N = 10 3 nodes and each node has an identical population of 500 agents.The weights of each link W α ij is randomly assigned following an homogeneous distribution in the range [1,50].
In Fig. 6 we show the diagrams for the SIR (top) and the SIS (bottom) where dots represent the results obtained for Monte Carlo simulations of the epidemic processes and the solid lines are for the solution of the Markovian equations.As in the case of networked metapopulations, we observe a perfect agreement between simulations and the numerical solution of Eqs. ( 12)- (15).From the physical point of view we observe that, while for all the cases mobility enhances the anticipation of the epidemic onset, the multiplex composed of an ER and a SF topologies yields an intermediate anticipation effect compared to those observed for ER-ER and SF-SF.This is an interesting result that differentiates what has been recently observed in epidemic processes in multiplex contact networks [35,36], where coupling L layers yields an overall epidemic threshold that is equal to the smallest threshold of the isolated layers or, in other words, the epidemic onset is driven by the largest of the maximum eigenvalues of the set of adjacency matrices that define the layers.It is clear that the case of metapopulations the situation is more complicated as we show in the following section.
We now focus on the general scenario in which λ αβ = λ, i.e., the contagion probability between two agents depends on their corresponding types.To this aim, we consider one population of agents whose movements are described by an ER mobility network and another population whose move-ments occur according to an SF graph.The number of patches is N = 10 3 and inside each patch there are 500 agents of each type (ER and SF).We consider the situation in which λ αβ λ αα (α = β).In particular contagion between agents moving in the ER layer occur with probability λ ER = 1.5µ/500 and that for the agents moving in the SF layer is set to λ SF = 1.1µ/500 (recall that µ/500 is the epidemic threshold for a well mixed population of 500 agents).In its turn, we have set the infection probability between agents of different type to λ ER−SF = λ SF −ER = 0.025µ/500.Finally, to work with a more heterogeneous setup, we study the case of an SIR dynamics in which a small seed of initial infected agents is set in a single patch and affect only agents of one type (here those moving across the ER layer).
To analyze the accuracy of Eqs. ( 12)-( 15) in capturing the spatio-temporal evolution of epidemics, we first consider the temporal evolution of the fraction of infected individuals of each type (layer).In panel (a) of Fig. 7, we show this evolution comparing the solution of the Markovian equations (solid lines) with the result obtained from Monte Carlo simulations (points).It is clear that the Markovian equations ( 12)-( 15) fairly reproduce the output of the numerical simulation, capturing the delay of the onset of the epidemics in the SF layer with respect to that in the population moving across the ER.This delay is a clear consequence of (i) the localization of the initial infected individuals in the ER layer and (ii) the small contagion probability between agents of different type (layer).Interestingly, the fact that λ ER−SF is far less than the threshold (µ/500) in a closed population of 500 agents does not prevent the disease from invading the SF layer.
Finally, in Fig. 7.(b)-(e) we show the temporal evolution of the fraction of recovered individuals for each patch in each of the layers (ER top and SF bottom) obtained from numerical simulations (left panels) and solving Eqs. ( 12)-( 15) (right panels).The fair agreement between left and right panels indicates the great spatio-temporal accuracy of the Markovian model.Here, in addition to the delay in the onset of the epidemics in the SF population already observed in (a), it is remarkable that two different stationary regimes are obtained in each layer.Namely, the fraction of recovered individuals in the ER layer is nearly identical for all the patches.However, in the SF population the stationary pattern points out a far more heterogeneous distribution of recovered individuals across the different patches.

B. Real Multiplex Metapopulations
To shed more light into the validity of the Markovian equations and, as a byproduct, to illustrate one relevant scenario where multiplex metapopulations capture contagion processes in real setups, we now study the SIS and SIR dynamics in an urban system (Medellín) where 6 different socio-economic classes coexist.Specifically, these social classes range from FIG. 8: Epidemic diagrams, I(λ) (Left) and R(λ) (Right) for SIR and SIS dynamics respectively in a real multiplex metapopulations.Solid lines denotes the predictions of our model about the incidence of a disease whereas black dots show the results obtained from averaging 20 realisation of numerical simulations.The colour code denotes the value of the mobility of the agents p.The recovery rate is set to µ = 0.2.FIG.9: Impact of an SIS disease, I(λ), on each of the layers of a real multiplex metapopulation.Note that the value of λ has been re-scaled by the critical value λc at p = 0.The mobility of the agents has been set from left to the right to p = 0, p = 0.2 and p = 1.Solid lines correspond to theoretical predictions obtained by iterating Eqs.12-15, whereas black dots are the result from averaging 20 realisations of numerical simulations.The recovery rate of both dynamics has been set to µ = 0.2.
1, which gather those inhabitants with the lowest incomes, to class 6, corresponding to the wealthiest individuals.The separation into 6 socioeconomic classes in Colombia [37] and, in particular, in large cities such as Medellín (the second largest city in Colombia with around 5 • 10 6 inhabitants) leads to a different demographic distribution across towns and, equally important, to different mobility patterns, as previously presented in [38,39].
To study the evolution of diseases while preserving the information related to the existence of different socio-economic classes, we make use of the former formalism by constructing, from the data presented in [40,41], a multiplex network of 6 layers.Each layer accounts for the mobility patterns as well as the distribution of the population which belong to each of the socioeconomic classes.Thus, the result is a multiplex network of 413 × 6 nodes, where 413 corresponds to the number of neighborhoods in which Medellin was divided for this study.
As in the case of synthetic networks, we first check the accuracy of our equations in predicting the total epidemic incidence.In Fig. 8 we plot the epidemic diagrams corresponding to different values of the degree of mobility, p, for both SIS and SIR diseases.These diagrams are obtained via Monte Carlo simulations (points) and by solving the Markovian equations (lines), showing an excellent agreement.Interestingly, we can also notice that, in Medellín, the agents mobility has a detrimental effect on the onset of epidemics, since the epidemic threshold increases with p.This counterintuitive behaviour, which we have already reported for monolayer configurations [30], emerges from the homogenization of the demographic distribution while increasing the mobility.
The particular architecture of the Medellin multiplex allows us to assess the effect of the mobility on the impact of dis-eases within each social class.For this purpose, we represent in Fig. 9, for several values of the mobility p, the epidemic diagrams for a SIS disease, I(λ).There we can find again a fair agreement between theory and simulations, what reveals that the multiplex Markovian equations allow to capture the incidence of epidemics at the layers' level.Remarkably, we observe that the critical point of the full multiplex moves towards high values (detrimental effect on the epidemic spreading [30]) as mobility increases, but this effect is not reported on the individual layers by themselves.
Remarkably, some features about the underlying multiplex network can be inferred from these graphs.For instance, the results corresponding to the static case (p = 0) unveil the demographic distribution of the layers.On one hand, in Fig. 9.a, we can observe that agents from socioeconomic classes 1, 2 and 3 occupy the most populated nodes, since the epidemic onset associated to these layers is the smallest one.On the other hand, it becomes clear that individuals from social class 6 reside practically isolated from the rest of the classes, occupying sparsely populated neighbourhoods.Besides, from Fig. 9.b-c, we can also notice that mobility promotes the social mixing, since by increasing p we obtain a more homogeneous impact of the disease across the layers.
To get more insight about the interaction among the different layers and to further validate our formalism, we now address the spatio-temporal propagation of diseases whose initial seed is localised inside one of the layers.For this purpose, we have fixed the parameters of our model (p, λ, µ) and represented in Fig. 10 the time evolution of the number of infected agents according to a SIR disease for each socioeconomic class when the seed is localised in classes 1 (Fig. 10.a) and 5 (Fig. 10.b).The solution of the Markovian equations captures the non-trivial interaction patterns between the different classes.In particular, it can be noticed that contagion processes take place mainly among close classes (in terms of incomes) since they show a cascade-like structure: 1 → 2, 3 → 4 in Fig. 10.a. and 5 → 4 → 3 → 4, 2 → 1 in Fig. 10.b.Finally, the nontrivial nature of the time evolution of infections is captured by the existence of a feedback phenomenon when looking to the sequence of local outbreaks for classes 2, 3, and 4. The observed correlations between layers' outbreaks reveals the closeness between the individuals in these middle class layers.

V. DEDUCTION OF THE EPIDEMIC THRESHOLD
The fair agreement between agent-based simulations and the solution of the Markovian equations allow us to make use of them in order to derive the analytical expression of the epidemic threshold.To get a first insight about the behavior of the epidemic threshold and its relation with the topology of the layers composing the multiplex network, we again consider that all the agents, regardless of their layers, interact in the same way so that λ αβ = λ ∀(α, β).For the sake of simplicity, we also focus on the SIS case (similar results are obtained for the SIR model).In this case, see Eq. ( 12), the stationary solution for the fraction of infected agents of type α associated to patch i, ρ α i , fulfills: As usual for calculating the threshold, we linearize the above expression by considering that the fraction of infected people in the stationary state is very small (ρ α: i = α i << 1 ∀α ∀i).This way, we can neglect second order terms in α i in Eq (14), so that P α * i is given by: Introducing this expression into Eq.(16), the stationary state of the epidemics can be written as: (18) By using the value of n β j→i from Eq. ( 15) and keeping up to first order in α i we finally obtain the expression: At this point, it becomes clear that Eq. ( 19) defines an eigenvalue problem for the feasible solutions α i .Indeed, there are N • L feasible solutions of λ corresponding to the eigenvalues of the N • L × N • L supra-matrix M.However, since we are interested in the minimum value λ c for which Eq. ( 19) is fulfilled, the epidemic threshold is thus associated to the largest eigenvalue of M as: Let us now describe the entries of the matrix M, see Eq. ( 19), since they allow us to quantify the microscopical interactions among agents across the multiplex metapopulations.
In fact, the elements M αβ ij correspond, close to the epidemic threshold, to the probability that an agent of type α associated to patch i contacts with another one of type β from patch j.Specifically, each element contains three contributions accounting for the three potential sources of infections that a healthy agent can find: from agents associated to the same node inside this node [weighted by (1 − p) 2 ], from agents from a different patch, either at one of the two patches they are associated to [weighed by p(1 − p)], and from agents with whom she contacts inside a third place different from their associated nodes (weighted by p 2 ).
To round off this derivation, let us remark that the generality of the expression for the epidemic threshold of multiplex metapopulations allows us to recover, by setting L = 1, the value of the epidemic threshold for diseases which propagates over single metapopulations.Indeed, for L = 1, Eq. ( 19) turns into: so that the epidemic threshold is given by: where M is now an N × N matrix.In the same fashion as supra-matrix M, each term M ij of matrix M encodes the probability that an agent associated to patch i contacts with another from patch j.
We have checked the validity of Eq. ( 20) by computing the largest eigenvalue of M for the three synthetic multiplexes under study in Fig. (6) for a range of values of p ∈ [0, 1].This way, through Eq. ( 20) we obtain a curve λ c (p), see Fig. 11, that reproduces the onset observed in Monte Carlo simulations.The monotonous decreasing behavior of curve λ c (p) corroborates that mobility enhances the spread of the disease.

VI. CONCLUSIONS
In this work, we have elaborated a theoretical formalism to analyze spreading processes in multiplex metapopulations characterized by recurrent mobility patterns.Our framework gets rid of the assumptions about the correlations between the node attributes and epidemic variables introduced in heterogeneous mean field formulations.This way, the formalism introduced here is general enough so to accommodate any origin-destination (weighted and directed) matrix containing different commuting patterns within a population and to cast the information about the local census of each patch.
First, we have introduced the Markovian evolution equations for the monoplex (single layer multiplex) case under the SIR and SIS dynamics.We have tested their validity by solving these equations and comparing their solution with the results obtained from Monte Carlo simulations in synthetic ER and SF networks.The agreement obtained is remarkable both at the macroscopic and the microscopic level.In particular, by solving the Markovian equations we can reproduce the spatio-temporal epidemic patterns capturing the onset of epidemics at the local level of patches.The second step has been to generalize the former formalism to address metapopulations composed of several types of agents whose mobility patterns are different.To this aim, we have made use of the multiplex formalism, thus constructing a multiplex metapopulation.We have again checked the validity of the Markovian formalism showing a great accuracy for both macroscopic and microscopic indicators.
The validity of the Markovian equations has allowed us to derive analytical expressions for the global epidemic threshold of multiplex metapopulations.Again, the analytical prediction is in complete agreement with numerical simulations.Interestingly, the onset is related to the maximum eigenvalue of a supra-matrix M in which the different mobility patterns, local census and the degree of mobility interplay.Remarkably, the structure of these supra-matrix M captures three ba-sic contagion processes for a healthy individual.
On more general grounds, dynamical processes on multiplexes have been a research focus in the recent years [42][43][44][45] and, in particular, their application to epidemics [46].As usual in the multiplex literature, the scenario considered is that of coupled contact networks, so that a node is an individual that interact in different ways (i.e. through different interaction layers) with the rest of the nodes.Under this setting, different problems such as the diffusion of a disease through different contagion channels [35,36,47], the cooperative spreading of different diseases [48][49][50][51] or the coevolution of different contagion processes [52,53] have been addressed.Here, at variance, the two interaction levels (epidemics and mobility) of the metapopulation yield interesting results related to the interplay of the architecture of layers.We have shown the applicability of the formalism to a real case study (city of Medellin) where we have gathered data of the mobility patterns for different socio-economical classes (layers).The results present an interesting behaviour of epidemic detriment [30] for the full multiplex structure while there is no epidemic detriment for all individual layers.
In a nutshell, the formalism introduced here provides with a reliable and computationally time-saving platform to ana-lyze the epidemic risk of systems displaying recurrent mobility patterns.This way, the formalism can be used to readily identify those critical areas that spur the unfolding of diseases.In addition, the possibility of handling analytical equations can be further exploited beyond the derivation of the epidemic threshold and combined together with control techniques to test in an efficient way different contention policies.We expect as well that our Markovian formalism can be further extended in the future to accommodate more sophisticated commuting patterns and more refined epidemic models, thus approaching more to real epidemic scenarios.

FIG. 1 :
FIG. 1: (color online).Schematic representation of one time step of the Movement-Interaction-Return (MIR) metapopulation model.The network is composed of N = 3 patches.At the movement stage some of the local agents decide to move to the other patches according to the probabilities encoded in matrix R. Once agents have moved interact in a well-mixed way and change their epidemic status (Healthy or Infected) according to an SIS model.Finally, the agents come back to their home patches and a new time step starts.

1 FIG. 2 :
FIG. 2: (color online) I(λ) for the SIS dynamics in ER (top) and SF (bottom) networks of 10 3 nodes and k = 5.5 and k = 7.3 respectively.The population of each node is 5 • 10 3 individuals.The solid curves indicate the solution obtained by solving the Markovian evolution equations (the color of each curve indicates the value of p as shown in the color bars), whereas the points correspond to the results obtained by using MC simulations (20 realizations for each value of λ).Note that the value of λ has been re-scaled by the critical value at p = 0, i.e., that of a well-mixed population of n = 5000 individuals: λc(p = 0) = µ/n = 4 • 10 −4 .The recovery rate is µ = 0.2.

1 FIG. 3 :
FIG. 3: R(λ) for the SIR dynamics in ER (top) and SF (bottom) networks of 10 3 nodes and k = 5.5 and k = 7.3 respectively.The population of each node is 5000 individuals.The solid curves indicate the solution obtained by solving the Markovian evolution equations (the color of each curve indicates the value of p as shown in the color bars), whereas the points correspond to the results obtained by using MC simulations (10 2 realizations for each value of λ).Note that the value of λ has been re-scaled by the critical value at p = 0, i.e., that of a well-mixed population of n = 5 • 10 3 individuals: λc(p = 0) = µ/n = 4 • 10 −4 .The recovery rate is µ = 0.2.

FIG. 4 :
FIG. 4: Microscopic evaluation of the Markovian equation.The top panel shows the time evolution, according to the agent based simulation, of the fraction of recovered agents of each node, ri(t) in an SF network of N = 200 when a initial infectious seed is placed at a single node.The population of the system is P = 7 • 10 5 agents while the size of a node is proportional to its out strength (see text).The weights of the links are randomly assigned following an homogeneous distribution in the range Wij ∈ [1, 50].The panel in the middle of the figure shows the same scenario but solved using the Markovian equations.Finally, to evaluate the accuracy of the Markovian equations we plot in the bottom panel the time evolution of the error [see Eq. (10)] with respect to the Monte Carlo simulation.The parameters used are: p = 0.1, λ = 2λc(p = 0) and µ = 0.2.

FIG. 5 : 1 FIG. 6 :
FIG. 5: Schematic representation of a metapopulation multiplex composed of L = 3 layers.Each of the N = 3 patches (nodes) is represented in each of the layers.The layers highlight that individuals of type α associated to patch i move to another patch j with probability R α ij which, in general, is different from the rate of transitions of agents of type β = α associated to the same node.This way each layer α presents a topology captured by a different matrix R α .

FIG. 7 :
FIG. 7: Spatio-temporal patterns of the SIR dynamics in a metapopulation multiplex composed of an ER and a SF layer.Each layer has 10 3 patches and 500 individuals are associated to each patch.Theoretical prediction are shown by lines whereas dots represent Monte Carlo results.The initial infected agents are placed in a single patch of the ER layer.This, together with the small contagion probability between agents of different layers (see the text for details), causes the time difference between the epidemic onsets in each layer as observed from panel (a).Panels (b)-(d) show the time evolution of the fraction of recovered agents for each patch.The top panels show this evolution in the ER layer obtained for Monte Carlo simulations [panel (b)] and the solution of the Markovian model [panel (c)].On the other hand, bottom panels show the same evolution in the SF layer as obtained again from simulations [panel (d)] and by solving the Markovian equations (12)-(15) [panel (e)].

FIG. 10 :
FIG. 10: Temporal evolution of a SIR disease whose seed is initially localised inside layer 1 (Left) and layer 5 (Right).Solid lines correspond to theoretical predictions according to Eqs.(12-15) whereas black dots are the output of Monte Carlo simulations.The mobility of the agents p, the contagion rate λ and the recovery rate µ have been set to (p, λ, µ)= (0.05,4λc(p=0),0.2).

FIG. 11 :
FIG. 11: Epidemic diagrams, I(λ, p) for SIS dynamics over those multiplex metapopulations used in Fig.6.The color-code (see color bar) denotes the fraction of infected individuals in the steady state as obtained from Monte Carlo simulations.The solid curves indicate the function λc(p) as dictated from Eq. (20) for the multiplex metapopulations.The recovery rate of the SIS dynamics has been set to µ = 0.2.From left to the right we have used SF-SF, ER-ER and ER-SF architectures to simulate the evolution of a disease.