# A framework for epidemic spreading in multiplex networks of metapopulations

###### Abstract

We propose a theoretical framework for the study of epidemics in structured metapopulations, with heterogeneous agents, subjected to 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 both the classical Susceptible-Infected-Susceptible and the Susceptible-Infected-Removed epidemic models within this framework, and compare macroscopic and microscopic indicators of the spreading process with extensive Monte Carlo simulations. Our results are in excellent agreement with the simulations. We also derive an exact expression of the epidemic threshold on this general framework revealing a non-trivial dependence on the mobility parameter.

###### pacs:

89.75.Hc, 89.75.Fb## 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 gleam1 (); gleam2 ().

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 hanski (); metapop1 (); metapop2 (); metapop3 (); metapop4 (). For the case of epidemic modeling, the usual metapopulation scenario metapop5 (); pion1 (); pion2 () 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 airline () and the recent advances of network epidemiology review () have motivated a thorough analysis about the impact that the structure of mobility networks has on the onset of global-scale contagions. In the last decade, important steps towards the inclusion of realistic mobility structures have been done colizza1 (); colizza2 (); colizza3 (); colizza4 (). These approaches had to compromise between realism and analytical feasibility. On one side, lengthy mechanistic simulationsgleam1 (); eubank () 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 mean-field 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 hmf (), 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 real-world diseases prediction. As data gathering techniques and epidemic surveillance influenzanet () increase their accuracy, metapopulation models face new challenges ball (). 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 commutes (); prx (); jtb (); epjb (), at the expense of considering either mean-field 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 kivela2014 (); boccaletti2014 (); dedomenico2013 (); battiston (); radicchi (), consist on a set of networks (usually called layers) and a set of nodes. Each node is represented once in each network layer allowing it to share different connectivity patterns in each of the 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 natphys (); scarpino (). 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 nodes (the patches) and a total population of agents. Importantly, each agent has associated one of the patches so that all the movements of agents associated to a node, say , initiate from and return to it. In its turn, a node of the network is the basement (or home) of a number of agents so that . 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 account for the weight of the interaction from node to node .

The dynamical model implemented in our metapopulation involves three different stages at each time step : movement, interaction and return (MIR). First, each agent decides whether moving with probability or remaining in its associated home node with probability . If the agent moves, it goes to any of the nodes connected to as dictated by the adjacency matrix . The probability that a patch is chosen, is proportional to the weight of the corresponding entry of the adjacency matrix:

(1) |

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 . 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 in order to construct the transition rates matrix . This way, the model has as control parameters the displacement probability 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:

(2) |

for the SIS model, and:

(3) |

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 EPL (); Guerra (); PRE () 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.

### iii.1 SIS model

For the SIS model, we have a set of variables denoting the fraction of infected agents associated to patch at time . It is important to stress that, according to the MIR model, an agent whose associated patch is can be in other node at time . The time evolution of can be written as:

(4) |

where the first term denotes the fraction of infected agents associated to that do not recover at time . The second term instead accounts for the fraction of healthy agents associated to that pass to infected at time . In this second term, is the probability that a healthy agent associated to node becomes infected at time . This probability reads:

(5) |

where the first term denotes the probability that a susceptible agent associated to patch becomes infected when remaining at its home node and the second one accounts for the probability that this agent catches the disease when moving to any neighbor of .

Finally, the probability in Eq. (5) denotes the probability that a healthy agent in (but not necessarily associated to) node at time becomes infected after the contact with any of the infected agents present inside at the same time. Then, probability reads:

(6) |

where:

(7) |

being when and 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 , and . In addition, matrix 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, , can be also set according to the local census of the population under study.

### iii.2 SIR dynamics

The formulation of the Markovian equations for a metapopulation under a SIR spreading dynamics demands to add another set of variables: (), i.e., the fraction of recovered agents associated to patch . Thus, the set of equations (4) for the SIS model is now substituted by the following set of equations:

(8) | |||||

(9) |

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 becomes infected at time , , 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).

### iii.3 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 and average connectivities and respectively. The nodes of these networks are homogeneously populated, , so that the total population of our systems is . The weights between the nodes of the graphs are randomly assigned following a homogeneous distribution within the range . Once all the weights are set, we construct the transition matrix [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 so that, on average, there is infected agent per node at time . This initial configuration corresponds to set as initial conditions of the Markovian equations (and 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 realizations for each combination of the parameters (, 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, , as a function of the infection probability, , for different movement probabilities . The points denote the results of Monte Carlo simulations for each value of and 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, , as function of for the same set of movement probabilities 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 (SIS) and (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 natphys (); scarpino (). 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 spatio-temporal evolution of the infections in the population. To this aim we now fix , and , 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 .

The results of the agent based simulation are shown in the top panel in Fig. 4 where we have used an SF network of nodes and a total population of agents where each patch contains a number of agents proportional to its out-strength (). We have ranked the nodes (from to ) according to the order in which infections occur, so that node is the one in which the initial infectious seed is placed. As shown, after a transient time of roughly 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:

(10) |

where and are, respectively, the number of recovered and infected agents associated to node at time observed in Monte Carlo simulation while and 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 () just after the avalanche of contagions across the metapopulation starts. After this peak, reaches a stationary value around 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 types of mobility patterns coexist within each patch. This way, the population of a patch is the sum of the number of agents of each type and the probability that an agent of patch and type visits another patch is now written as the generalization of Eq. (1):

(11) |

where is associated to the number of observed trips of agents of type in patch to patch .

To analyze this situation, it is natural to make use of a multiplex formulation kivela2014 (); boccaletti2014 (); dedomenico2013 (); battiston (); radicchi () 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 () and the architecture of each layer is described by a different matrix . Each patch of the system is represented as one node in each network layer and the corresponding 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 with respect to the networked metapopulation. In particular for the SIS (and SIR) model, the variables are (and ), which denote the fraction of infected (and recovered) individuals of layer associated to node . In this case, SIR equations become:

(12) | |||||

(13) |

while for the SIS model we only have Eq. (12) with . The term , which denotes the probability that an agent of type placed in patch at time becomes infected, reads:

(14) |

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 that travel to a different patch is given by:

(15) |

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.

### iv.1 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, . We have studied three types of multiplex of layers, namely ER-ER, SF-SF, and ER-SF, of nodes and each node has an identical population of agents. The weights of each link is randomly assigned following an homogeneous distribution in the range .

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 cozzo (); arruda (), where coupling 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 movements occur according to an SF graph. The number of patches is and inside each patch there are 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 and that for the agents moving in the SF layer is set to (recall that is the epidemic threshold for a well mixed population of agents). In its turn, we have set the infection probability between agents of different type to . 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 is far less than the threshold () in a closed population of 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.

### iv.2 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 different socio-economic classes coexist. Specifically, these social classes range from 1, which gather those inhabitants with the lowest incomes, to class 6, corresponding to the wealthiest individuals. The separation into socioeconomic classes in Colombia colombia () and, in particular, in large cities such as Medellín (the second largest city in Colombia with around inhabitants) leads to a different demographic distribution across towns and, equally important, to different mobility patterns, as previously presented in socio1 (); socio2 ().

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 source (); data (), a multiplex network of 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 nodes, where 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, , 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 . This counter-intuitive behaviour, which we have already reported for monolayer configurations natphys (), 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 diseases within each social class. For this purpose, we represent in Fig. 9, for several values of the mobility , the epidemic diagrams for a SIS disease, . 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 natphys ()) 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 () 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 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 (, , ) 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: in Fig. 10.a. and 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 , , and . 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 , , fulfills:

(16) |

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 (). This way, we can neglect second order terms in in Eq (14), so that is given by:

(17) |

Introducing this expression into Eq.(16), the stationary state of the epidemics can be written as:

(18) |

By using the value of from Eq. (15) and keeping up to first order in we finally obtain the expression:

(19) |

At this point, it becomes clear that Eq. (19) defines an eigenvalue problem for the feasible solutions . Indeed, there are feasible solutions of corresponding to the eigenvalues of the supra-matrix . However, since we are interested in the minimum value for which Eq. (19) is fulfilled, the epidemic threshold is thus associated to the largest eigenvalue of as:

(20) |

Let us now describe the entries of the matrix , see Eq. (19), since they allow us to quantify the microscopical interactions among agents across the multiplex metapopulations. In fact, the elements correspond, close to the epidemic threshold, to the probability that an agent of type associated to patch contacts with another one of type from patch . 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 ], from agents from a different patch, either at one of the two patches they are associated to [weighed by ], and from agents with whom she contacts inside a third place different from their associated nodes (weighted by ).

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 , the value of the epidemic threshold for diseases which propagates over single metapopulations. Indeed, for , Eq. (19) turns into:

(21) |

so that the epidemic threshold is given by:

(22) |

where is now an matrix. In the same fashion as supra-matrix , each term of matrix encodes the probability that an agent associated to patch contacts with another from patch .

We have checked the validity of Eq. (20) by computing the largest eigenvalue of for the three synthetic multiplexes under study in Fig. (6) for a range of values of . This way, through Eq. (20) we obtain a curve , see Fig. 11, that reproduces the onset observed in Monte Carlo simulations. The monotonous decreasing behavior of curve 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 in which the different mobility patterns, local census and the degree of mobility interplay. Remarkably, the structure of these supra-matrix captures three basic contagion processes for a healthy individual.

On more general grounds, dynamical processes on multiplexes have been a research focus in the recent years diffusion (); games (); patterns (); synchronization () and, in particular, their application to epidemics epid1 (). 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 cozzo (); buono (); arruda (), the cooperative spreading of different diseases sanz (); sahneh (); nahid (); fakhteh () or the coevolution of different contagion processes granell1 (); granell2 () 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 natphys () 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 analyze 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.

## Vii Acknowledgements

We acknowledge useful discussions and suggestion by S. Meloni. Financial support came from MINECO (through projects FIS2015-71582-C2 and FIS2014-55867-P), from the Departamento de Industria e Innovación del Gobierno de Aragón y Fondo Social Europeo (FENOL group E-19). AA acknowledges ICREA Academia and the James S. McDonnell Foundation.

## References

- (1) D. Balcan, B. Gonçalves, H. Hu, J. J. Ramasco, V. Colizza and A.Vespignani. Modeling the spatial spread of infectious diseases: The GLobal Epidemic and Mobility computational model. Journal of Computational Science 1, 132–145 (2010).
- (2) M. Tizzoni et al. Real-time numerical forecast of global epidemic spreading: case study of 2009 A/H1N1pdm. BMC Medicine 10 165 (2012).
- (3) I. Hanski. Metapopulation dynamics. Nature 396, 41–49 (1998).
- (4) I. Hanski, and M.E. Gilpin. Metapopulation Biology: Ecology, Genetics, and Evolution. (Academic Press, 1997).
- (5) D. Tilman and P. Kareiva. Spatial Ecology. (Princeton University Press, 1997).
- (6) J. Bascompte, and R.V. Solé. Modeling Spatiotemporal Dynamics in Ecology. (Springer, 1998).
- (7) I. Hanski, and O.E. Gaggiotti. Ecology, Genetics, and Evolution of Metapopulations.(Elsevier and Academic Press, 2004).
- (8) M.J. Keeling, and P. Rohani. Modeling Infectious Diseases in Humans and Animals. (Princeton University Press, 2008).
- (9) L. Sattenspiel, and K. Dietz. A structured epidemic model incorporating geographic mobility among regions. Mathematical Biosciences 128, 71–91 (1995).
- (10) B. Grenfell, and J. Harwood. (Meta)population dynamics of infectious diseases. Trends in Ecology & Evolution 12, 395–399. (1997).
- (11) R. Guimerá, S. Mossa, A. Turtschi, and L. A. N. Amaral The worldwide air transportation network: Anomalous centrality, community structure, and cities’ global roles. Proc. Nat. Acad. Sci (USA) 102 7794–7799 (2005).
- (12) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani. Epidemic processes in complex networks. Rev. Mod. Phys. 87, 925–979 (2015).
- (13) V. Colizza, R. Pastor-Satorras, and A. Vespignani. Reaction–diffusion processes and metapopulation models in heterogeneous networks. Nature Phys. 3, 276–282 (2007).
- (14) V. Colizza, and A. Vespignani. Invasion threshold in heterogeneous metapopulation networks. Phys. Rev. Lett. 99, 148701 (2007).
- (15) V. Colizza, and A. Vespignani. Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: Theory and simulations. Journal of theoretical biology 251, 450–467 (2008).
- (16) D. Balcan, V. Colizza, B. Gonï¿½alves, H. Hao, J.J. Ramasco, and A. Vespignani. Multiscale mobility networks and the spatial spreading of infectious diseases- Proc. Nat. Acad. Sci (USA) 106, 21484–21489 (2009).
- (17) S. Eubank, et al. Modelling disease outbreaks in realistic urban social networks. Nature 429, 180–184 (2004).
- (18) R. Pastor-Satorras, and A. Vespignani. Epidemic Spreading in Scale-free networks. Phys. Rev. Lett. 86, 3200 (2001).
- (19) C. Guerrisi et al. Participatory Syndromic Surveillance of Influenza in Europe. Journal of Infectious Diseases 214 S386-S392 (2016).
- (20) F. Ball et al. Seven challenges for metapopulation models of epidemics, including households models. Epidemics 10, 63–67 (2015).
- (21) D. Balcan, and A. Vespignani. Phase transitions in contagion processes mediated by recurrent mobility patterns. Nature Phys. 7, 581–586 (2011).
- (22) V. Belik, T. Geisel and D. Brockmann. Natural Human Mobility Patterns and Spatial Spread of Infectious Diseases. Phys. Rev. X 1, 1, 011001 (2011).
- (23) D Balcan, and A Vespignani. Invasion threshold in structured populations with recurrent mobility patterns. J. Theor. Biol. 293, 87–100 (2011).
- (24) V. Belik, T. Geisel and D. Brockmann. Recurrent host mobility in spatial epidemics: beyond reaction-diffusion. Eur. Phys. J. B 84, 579–587 (2011).
- (25) M. Kivelä, A. Arenas, M. Barthelemy, J.P. Gleeson, Y. Moreno, and M.A. Porter. Multilayer Networks. J. Complex Networks 2, 203-271 (2014)
- (26) S. Boccaletti, G. Bianconi, R. Criado, C.I. del Genio, J. Gómez-Gardeñes, M. Romance, I. Sendiña-Nadal, Z. Wangk, and M. Zanin. The structure and dynamics of multilayer networks. Phys. Rep. 544, 1–122 (2014)
- (27) M. De Domenico, A. Solé-Ribalta, E. Cozzo, M. Kivelä, Y. Moreno, M.A. Porter, S. Gómez and A. Arenas. Mathematical formulation of multi-layer networks. Phys. Rev. X 3, 041022 (2013).
- (28) F. Battiston, V. Nicosia, and V. Latora. Structural measures for multiplex networks. Phys Rev E 89, 032804 (2014).
- (29) F. Radicchi, and A. Arenas. Abrupt transition in the structural formation of interconnected networks. Nature Phys. 9, 717–720 (2013).
- (30) J. Gómez-Gardeñes, D. Soriano-Paños, and A. Arenas. Critical regimes driven by recurrent mobility patterns of reaction-diffusion processes in networks. Nature Phys. (doi: 10.1038/s41567-017-0022-7) (2018).
- (31) S. V. Scarpino. Don’t close the gates. Nature Phys. (doi: 10.1038/s41567-017-0028-1) (2018).
- (32) S. Gómez, A. Arenas, J. Borge-Holthoefer, S. Meloni, and Y. Moreno. Discrete-time Markov chain approach to contact based diseases spreading in complex networks. EPL 89, 38009 (2010).
- (33) B. Guerra, and J. Gómez-Gardeñes, Annealed and mean-field formulations of disease dynamics on static and adaptive networks. Phys. Rev. E 82, 035101 (2010).
- (34) S. Gómez, J. Gómez-Gardeñes, Y. Moreno and A. Arenas. Nonperturbative heterogeneous mean-field approach to epidemic spreading in complex networks. Phys. Rev. E 84 036105 (2011).
- (35) E. Cozzo, R.A. Baños, S. Meloni, and Y. Moreno. Contact-based social contagion in multiplex networks. Phys. Rev E 88, 050801 (2013).
- (36) G.F. de Arruda, E. Cozzo, T.P. Peixoto, F.A. Rodrigues, and Y. Moreno. Disease Localization in Multilayer Networks. Phys. Rev. X 7, 011014 (2017).
- (37) C. Medina, L. Morales, R. Bernal, and M. Torero. Stratification and Public Utility Services in Colombia: Subsidies to Households or Distortion of Housing Prices? Economia 7, 41–99 (2007).
- (38) L. Lotero, R.G. Hurtado, L.M. Floría, and J. Gómez-Gardeñes. Rich do not rise early: Spatio-temporal patterns in the Mobility Networks of different Socio-economic classes. Royal Society Open Science 3, 150654 (2016).
- (39) L. Lotero, A. Cardillo, R.G. Hurtado, and J. Gómez-Gardeñes. Several Multiplexes in the same city: The role of Wealth differences in urban mobility, in Interconnected Networks (Springer, 2016), pp. 149–164.
- (40) Universidad Nacional de Colombia and AREA Metropolitana del Valle de Aburrá 2006. Encuesta origen destino de viajes 2005 del Valle de Aburrá, estudios de tránsito complementarios y validación. (Technical Report).
- (41) L. Lotero, R.G. Hurtado, L.M. Floría, and J. Gómez-Gardeñes. Rich do not rise early: Spatio-temporal patterns in the Mobility Networks of different Socio-economic classes. Dryad Digital Repository. (http://dx.doi.org/10.5061/dryad.hj1t4).
- (42) S. Gómez, et al. Diffusion Dynamics on Multiplex Networks. Phys. Rev. Lett. 110, 028701 (2013).
- (43) J. Gómez-Gardeñes, I. Reinares, A. Arenas, LM Floría. Evolution of cooperation in multiplex networks. Sci. Rep. 2, 620 (2012).
- (44) N.E. Kouvaris, Shigefumi Hata, and Albert Díaz-Guilera. Pattern formation in multiplex networks. Sci. Rep. 5, 10840 (2015).
- (45) C.I. Del Genio, J. Gómez-Gardeñes, I. Bonamassa, and S. Boccaletti. Synchronization in networks with multiple interaction layers. Science Adv. 2, e1601679 (2016).
- (46) A. Saumell-Mendiola, M.A. Serrano, and M. Boguñá. Epidemic spreading on interconnected networks. Phys. Rev. E 86, 026106 (2012).
- (47) C. Buono, L.G. Alvarez-Zuzek, P.A. Macri, and L.A. Braunstein. Epidemics in partially overlapped multiplex networks. PloS One 9, e92200 (2014).
- (48) J. Sanz, Ch.-Y. Xia, S. Meloni, and Y. Moreno. Dynamics of Interacting Diseases. Phys. Rev. X 4, 041005 (2014).
- (49) F.D. Sahneh, and C. Scoglio. Competitive epidemic spreading over arbitrary multilayer networks. Phys. Rev. E 89, 062817 (2014).
- (50) W. Cai, L. Chen, F. Ghanbarnejad, and P. Grassberger. Avalanche outbreaks emerging in cooperative contagions. Nature Phys. 11, 936 (2015).
- (51) N. Azimi-Tafreshi. Cooperative epidemics on multiplex networks. Phys. Rev. E 93, 042303 (2016).
- (52) C. Granell, S. Gómez, and A. Arenas. Dynamical Interplay between Awareness and Epidemic Spreading in Multiplex Networks. Phys. Rev. Lett. 111, 128701 (2013).
- (53) C. Granell, S. Gómez, and A. Arenas. Competing spreading processes on multiplex networks: Awareness and epidemics. Phys. Rev. E 90, 012808 (2014).