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.

Popular Summary

The spread of disease is strongly driven by how humans move and segregate across cities, regions, and countries. Understanding how epidemics arise from the interplay of these behaviors is crucial for implementing efficient containment and prevention policies. In this work, we propose a theoretical framework for the spread of pathogens that incorporates the real patterns of mobility, demographics, and human diversity observed in urban environments.

Our framework represents various populations in a community as layers in a multiplex network, where nodes correspond to geographical areas and layers account for movement patterns. With this framework, we analyze a real urban system, the city of Medellin in Colombia, considering the interplay between six socioeconomic classes displaying disparate demographic segregation and mobility habits. Our framework can identify those neighborhoods and classes that trigger epidemic outbreaks.

Remarkably, we observe that small changes to both mobility and social mixing among these subpopulations can trigger abrupt changes. As a consequence, containment strategies targeting a certain neighborhood can quickly change from efficient to useless. Finally, we provide a physical interpretation of the abrupt changes of those patches driving the unfolding of epidemics based on the effective number of contacts of agents.

Our framework provides a reliable, time-saving platform for analyzing the spread of pathogens among various populations, allowing researchers to identify areas critical to the unfolding of diseases. With further improvements, our formalism could be extended to accommodate more sophisticated commuting patterns and epidemic models.

Figure 1

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 $\mathbf{R}$. Once agents have moved they interact in a well-mixed way and change their epidemic status (healthy or infected) according to a SIS model. Finally, the agents come back to their home patches and a new time step starts.

Figure 2

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 $\alpha$ associated to patch $i$ move to another patch $j$ with probability $R_{ij}^{\alpha }$ which, in general, is different from the rate of transitions of agents of type $\beta \ne \alpha$ associated to the same node. This way, each network layer $\alpha$ presents a topology captured by a different matrix ${\mathbf{R}}^{\alpha }$.

Figure 3

Epidemic diagrams for the SIR (top), $R(\lambda )$ and SIS (bottom), $I(\lambda )$ dynamics of three different multiplexes with $L=2$ layers. From left to right we have ER-ER, ER-SF, and SF-SF. In all the cases each network layer has $N={10}^3$ nodes and each node contains 500 individuals per layer. 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 agent-based simulations (50 realizations for each value of $\lambda$ and $p$). Note that the value of $\lambda$ has been rescaled by the critical value ${\lambda }_c$ at $p=0$, i.e., that of a well-mixed population of $n={10}^3$ individuals: ${\lambda }_c=\mu /{10}^3$ at $p=0$. The recovery rate is $\mu =0.2$.

Figure 4

Spatiotemporal 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 (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 (b) and the solution of the Markovian model (c). On the other hand, bottom panels show the same evolution in the SF layer as obtained again from simulations (d) and by solving the Markovian equations (11)–(14) (e).

Figure 5

Epidemic diagrams $I(\lambda ,p)$ for SIS dynamics over those multiplex metapopulations used in Fig. 3. 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 ${\lambda }_c(p)$ as dictated from Eq. (20) for the multiplex metapopulations. The values of interlayer coupling are $\gamma =1$ (top) and $\gamma =0$ (bottom). The recovery rate of the SIS dynamics has been set to $\mu =0.2$. From left to right we have used ER-ER, ER-SF, and SF-SF architectures to simulate the evolution of a disease.

Figure 6

Mobility networks of each socioeconomic class in the city of Medellín. Each panel shows the geographical location of each subpopulation (node) in the city of Medellín. The size of each node is proportional to the number of agents of the corresponding socioeconomic class in the subpopulation. The connections between two nodes denotes the existence of back-and-forth movements between two subpopulations for a given socioeconomic class.

Figure 7

Impact of a SIS disease $I(\lambda )$ on each of the layers of a real multiplex metapopulation. Note that the value of $\lambda$ has been rescaled by the critical value ${\lambda }_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. (11)–(14), whereas black dots are the result from averaging 20 realizations of numerical simulations. The recovery rate of both dynamics has been set to $\mu =0.2$.

Figure 8

Temporal evolution of a SIR disease whose seed is initially localized inside layer 1 (top) and layer 5 (bottom). Solid lines correspond to theoretical predictions according to Eqs. (11)–(14), whereas black dots are the output of Monte Carlo simulations. The mobility of the agents $p$, the contagion rate $\lambda$, and the recovery rate $\mu$ have been set to ($p,\lambda ,\mu )=[0.05,4{\lambda }_c(p=0),0.2]$.

Figure 9

(a) Epidemic threshold (color code) as a function of the agents’ mobility $p$ and the interlayer coupling $\gamma$. (b) Epidemic threshold as a function of the mobility for several values of the interlayer coupling. (c) Epidemic threshold as a function of the interlayer coupling for several values of the mobility. Note that the epidemic threshold has been rescaled by the epidemic threshold corresponding to the static case and decoupled layers, so that ${\tilde{\lambda }}_c={\lambda }_c(p,\gamma )/{\lambda }_c(0,0)$.

Figure 10

(a),(b) Magnitude of the components of the leading eigenvector (color coded) as a function of the social mixing $\gamma$ for every subpopulation (patch); see Eq. (23). The mobility of the agents has been set to (a) $p=0.0$ and (b) $p=0.6$. (c) Same as (a) and (b), but now fixing the interlayer coupling to $\gamma =0.1$ and monitoring the evolution of the leading eigenvector with $p$. (d) Inverse participation ratio (IPR) (color code) according to Eq. (24) as a function of $p$ and $\gamma$. Note that there are abrupt changes in the IPR for certain values of ($p$, $\gamma$). These strong variations encode the delocalization processes that take place when the dominant patch, which triggers the epidemic onsets, changes.

Figure 11

(a),(b) Average number of contacts $C_i^{\alpha }$ in which agents from patch $i$ at layer $\alpha$ participate as a function of $\gamma$. For the sake of clarity, only the largest values of this indicator, which are relevant for identifying the critical patches, have been represented. The chosen values of the mobility are $p=0$ in (a) and $p=0.6$ in (b). (c) Average number of contacts $C_i^{\alpha }$ in which agents from patch $i$ at layer $\alpha$ participate as a function of the mobility for $\gamma =0.1$.

Figure 12

$I(\lambda )$ 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\times {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 Monte Carlo simulations (20 realizations for each value of $\lambda$). Note that the value of $\lambda$ has been rescaled by the critical value at $p=0$, i.e., that of a well-mixed population of $n=5000$ individuals: ${\lambda }_c(p=0)=\mu /n=4\times {10}^{-4}$. The recovery rate is $\mu =0.2$.

Figure 13

$R(\lambda )$ 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 Monte Carlo simulations ($10^2$ realizations for each value of $\lambda$). Note that the value of $\lambda$ has been rescaled by the critical value at $p=0$, i.e., that of a well-mixed population of $n=5\times {10}^3$ individuals: ${\lambda }_c(p=0)=\mu /n=4\times {10}^{-4}$. The recovery rate is $\mu =0.2$.

Figure 14

Epidemic diagrams $I(\lambda )$ (top) and $R(\lambda )$ (bottom) for SIR and SIS dynamics, respectively, in a real multiplex metapopulation. Solid lines denote the predictions of our model about the incidence of a disease, whereas black dots show the results obtained from averaging 20 realizations of numerical simulations. The color code denotes the value of the mobility of the agents $p$. The recovery rate is set to $\mu =0.2$.