Metapopulation epidemic models with a universal mobility pattern on interconnected networks.

The global pandemic of the coronavirus disease 2019 (COVID-19) exemplifies the influence of human mobility on epidemic spreading. A framework called the movement-interaction-return (MIR) model is a model to study the impact of human mobility on epidemic spreading. In this paper, we investigate epidemic spreading in interconnected metapopulation networks. Specifically, we incorporate the human mobility pattern called the radiation model into the MIR model. As a result, the proposed model is more realistic in comparison to the original MIR model. We use the tensorial framework to develop Markovian equations that describe the dynamics of the proposed model on interconnected metapopulation networks. Then we derive the corresponding epidemic thresholds by converting tensors into matrices. Comprehensive numerical simulations confirm our analysis.

Introduction

The advent of various complex network models [1], [2] fills researchers with enthusiasm to investigate epidemic spreading. In the classical models of epidemic spreading in complex networks, nodes represent populations and edges represent human–human interactions. Then either mean-field approaches or generating function approaches can be used to describe epidemic processes and derive epidemic thresholds [3]. In some situations, nodes, instead of representing an individual, represent entities that host multiple individuals. Thus a framework using reaction–diffusion processes in metapopulation networks was introduced where each node contains a non-negative number of individuals and each individual in a node can diffuse along the edges [4]. Later, metapopulation epidemiological models were widely explored [5], [6], [7], [8], [9], [10], [11], [12]. Recently, Paños et al. [13] improved the reaction–diffusion-based models and then introduced the MIR model to characterize the dynamics of epidemic spreading in metapopulation networks. To further improve the MIR model and generalize the study of Paños et al. [13] to the classes of complex networks other than multiplex networks, it is necessary to focus on the following two key factors.

The first critical ingredient is human mobility or migration. Human mobility may significantly impact contagion propagation in reality. For example, human mobility is a crucial factor that leads to the current outbreak of coronavirus disease 2019 (COVID-19) [14] worldwide. The studies of human mobility or migration patterns attract the attention of researchers [15], [16], [17], [18]. In the past, the gravity model was generally used to describe mobility and migration patterns such as human mobility patterns. However, limitations exist for the gravity model. To avoid the limitations of the gravity model, Simini et al. [19] introduced the radiation model. The radiation model has superior performance to predict mobility and transport processes in the real world. The mobility pattern in the MIR model is only related to the weights of edges. Then a suitable human mobility pattern should be incorporated into the MIR model. To tackle this problem, we replace the original mobility pattern in the MIR model with the radiation model.

The other key factor is the network topology. Besides the classical complex network models, researchers pay special attention to multilayer networks [20]. Epidemic spreading in multiplex networks [21], [22], [23], interconnected networks [24], [25], [26], temporal networks [27], and networks with the core–periphery structure [28] was extensively studied. For multilayer networks, it is inappropriate to use adjacency matrices to represent connections. Then a tensorial framework was introduced to describe the structure of multilayer networks [29] and epidemic processes in multilayer networks [30]. Paños et al. [13] mainly investigate metapopulation epidemic models in a multiplex network, which is a special multilayer network. To generalize the previous research, we apply the tensorial framework to epidemic models on interconnected metapopulation networks. The presented work is straightforward to extend to other types of multilayer networks.

In summary, our contributions in the paper are as follows. To improve the MIR model and develop a more realistic metapopulation epidemic model, we incorporate the radiation model into the MIR model. To generalize the study of Paños et al. [13], we use the tensorial framework to develop Markovian equations that describe the dynamics of the proposed metapopulation epidemic model in interconnected networks. Then we derive the epidemic threshold using the presented Markovian equations.

The organization of the paper is as follows. In Section 2, we study epidemic spreading in metapopulation networks. In Section 3, we investigate epidemic spreading in interconnected metapopulation networks. Then we present simulation results in Section 4. Finally, we present the conclusion in Section 5.

Epidemic processes in metapopulation networks

We list important symbols introduced in Section 2 in Table 1.

Table 1

Summary of important symbols introduced in Section 2.

Symbol	Description
n_i	Number of agents in node i
p	Probability that an agent decides to relocate
R_ij	Probability that an agent moves from node i to node j
λ	Infection rate
μ	Recovery rate
I_i(t)	Fraction of infected agents of node i at time t
Δ_i(t)	Probability that a susceptible agent in home node i is infected at time t
P_i(t)	Probability that a susceptible agent in node i at time t is infected at the same time
n_ji	Expected number of agents that move to node i from node j
λ_c	Epidemic threshold

The model

Let be a network with the set of nodes and the set of edges . In addition, let be the number of nodes in . In the metapopulation framework, each node contains agents, and node is regarded as the home node of the agents in . We incorporate the radiation model that describes human mobility into the MIR model. The description of our model at each time step is as follows. At first, each agent either relocates with probability or stays in its associated home node with probability . If an agent of node relocates, then it moves to node () with probability . We denote as the number of agents in the circle of radius centered at (excluding agents in the source and destination). Then is expressed as where is the distance from node to node . The quantity is an important factor that is used to compute the number of agents moving from node to node per unit time in the radiation model. Based on the radiation model, probability is defined as where is the normalization constant. Then let be an matrix with entries . After each agent has completed the movement, the infection process takes place such that any agent can only interact with other agents in the same node. Each infected agent tries to infect each susceptible agent in the same node with probability . Moreover, each infected agent recovers to the susceptible state with probability . At last, agents come back to their associated home node after the completion of the infection process at time . The above process at each time repeats until the process enters into the stationary state.

Markovian equations

The presented Markovian equations in this section is similar to that in [13]. For the completeness, we present all relevant equations. Let be the fraction of infected agents of node at time , then it satisfies where is the probability that a susceptible agent in home node is infected at time . Here is expressed as where () is the probability that a susceptible agent in node () at time is infected at the same time. The probability satisfies where is the expected number of agents that move to node from node . The quantity is expressed as where represents the Kronecker delta.

The epidemic threshold based on Markovian equations is derived as follows. Let be in the stationary state. Then Eq. (3) becomes where represents in the stationary state. Replacing with the right hand side of Eq. (4) in the stationary state, Eq. (7) becomes where denotes in the stationary state. If is close to the threshold, then every is close to zero. Thus we ignore the second-order terms related to and then obtain Based on Eqs. (5), (9), we have Replacing with the right hand side of Eq. (6), Eq. (10) becomes Based on Eq. (11), Eq. (8) becomes (ignore second order terms) We define matrix whose entries satisfy where is the transpose of matrix . Thus epidemic threshold can be expressed as where is the largest eigenvalue of .

Epidemic processes in interconnected metapopulation networks

We list important symbols introduced in Section 3 in Table 2.

Table 2

Summary of important symbols introduced in Section 3.

Symbol	Description
n^il	Number of agents in home node i of layer l
Rjqil	Probability that an agent moves from node i of layer l to node j of layer q
λ_lq	Infection rate for an agent from layer l infects an agent from layer q
μ_l	Recovery rate for an agent in layer l
I^il(t)	Fraction of infected agents associated to home node i in layer l at time t
Δ^il(t)	Probability that a susceptible agent associated to home node i in layer l is infected at time t
P^il(t)	Probability that a susceptible agent in node i of layer l is infected in the same node at time t
niljq	Expected number of agents from node j in layer q move to node i in layer l

An interconnected network consists of a collection of layers where each layer ( is the set of nodes and is the set of edges) is a network. In an interconnected network, if node in layer is connected to a node in the same layer, the corresponding edge is called an intralayer edge; if in layer is connected to a node in layer with , then the corresponding edge is called an interlayer edge. In this section, we extend previous studies of epidemic spreading in multiplex metapopulation networks. Specifically, we investigate epidemic spreading in interconnected metapopulation networks. Fig. 1 illustrates the difference between interconnected networks and multiplex networks.

Fig. 1

Schematic illustration of interconnected networks and multiplex networks. (a) An interconnected network with two layers. (b) A multiplex network with two layers.

We can use tensors and tensor products to represent interconnected networks. Let be the size of for . Then we define as Let and be vector spaces where is the set of real numbers. The dual spaces and are identical to and respectively. Further, let be the canonical basis of and let be the canonical basis of . Here, and with satisfy Similarly, let be the canonical basis of and let be the canonical basis of . Then and with satisfy The tensor product is a vector space and each element of a basis in has the form where and . We use tensor of order on to represent the weights of edges. In particular, tensor is an element of the tensor product where is the dual space of . Then tensor is defined as where 1 is the coordinate of . Notice that we adopt the Einstein summation convention for the representation of tensors and vectors by the basis. Here, represents the weight of the edge between node in layer and node in layer provided that the edge exists in the network. Otherwise, we set .

Suppose that each node in an interconnected network contains a non-negative number of agents and let be the total number of agents in the network. Further, let be the number of agents in home node of layer and we set if . Then we define vector in space as Based on Eq. (19), we have .

The description of the model on interconnected networks at each time step is as follows. At first, each agent either relocates with probability or stays in its associated home node with probability . If an agent of node in layer relocates, then it moves to node in layer with probability . Given a network, we define as where is the distance from node in layer to node in layer . In Eq. (20), means that node in layer is neither node in layer nor node in layer . Based on the radiation model, is then defined as where is the normalization constant. We define associated tensor as After all the agents have completed the movements, the infection process takes place such that any agent can only interact with other agents in the same node. For agents in the same node, each infected agent whose home node is in layer tries to infect each susceptible agent whose home node is in layer with probability . Besides, each agent whose home node is in layer recovers to the susceptible state with probability . At last, agents come back to their associated home node after the completion of the infection process at time . The above process at each time repeats until the process enters into the stationary state.

In this section, we use the tensorial framework and Markovian equations to describe the dynamics of epidemic spreading. To avoid complication, we restrict to the case for any and in the analysis. Let be the fraction of infected agents associated to home node in layer at time . We define vector in as The equation for is expressed as where is the probability that a susceptible agent associated to home node in layer is infected at time . The probability is expressed as where is the probability that a susceptible agent in node of layer is infected by any infected agent in the same node at time . The vector corresponding to Eq. (25) in is defined as The probability is expressed as where is the expected number of agents, which move to node in layer , associated to home node in layer . The vector corresponding to Eq. (27) is defined as The quantity in Eq. (27) satisfies Similar to the Kronecker delta, the quantity in Eq. (29) is defined as The tensor corresponding to Eq. (29) is defined as The tensor corresponding to Eq. (30) is defined as

Epidemic thresholds

To simplify the analysis, we assume for each . Before deriving the epidemic threshold, we transform tensors into matrices. To begin with, we define a special function from the set ( is the set of natural numbers) to non-negative integers. In particular, function is defined as For convenience, we use to represent the value of . Then we define (corresponding to tensor ) as a matrix in whose entries satisfy Similarly, we define (corresponding to tensor ) as a matrix in whose entries satisfy Based on Eq. (30), matrix is the identity matrix in . We use , and to represent the vectors , and in the stationary state respectively. Substituting the right hand side of Eq. (25) into Eq. (24) in the stationary state, we have where and are the entries of vectors and respectively. If is close to the threshold, the entries of vector are close to zero. Then ignoring the second order terms related to the entries of , we have Inserting the right hand side of Eq. (37) into Eq. (27) and ignoring the second order terms, we have Based on Eqs. (29), (38), Eq. (36) becomes (we ignore second order terms) We replace the coordinates of the tensors with the corresponding entries of the matrices defined in Eqs. (34), (35), then we obtain where (or ) represents an entry of matrix (or ). Let be a matrix whose entries are defined as Based on Eqs. (40), (41), we have Then epidemic threshold satisfies where is the largest eigenvalue of .

Experimental results

Simulation results for single layer networks

To verify our theoretical analysis in Section 2, we run numerical simulations on metapopulation networks. We investigate two widely used complex networks including BA networks with parameter ( represents the number of edges to add each time) [31] and WS networks with parameter ( represents the number of neighbors of each node) and ( represents the rewire probability) [32]. In the experiment, the number of nodes of any network is fixed to and recovery rate is fixed to 0.5. Further, each node of the given network contains the same number of agents. Then we randomly choose an agent in every node as the seed node initially. In other words, every node initially contains an infected agent. Next, we run the agents-based simulation corresponding to the model in Section 2.1. Here we generate networks and then simulate the epidemic process times for each generated network. The final result is the average fraction of infected agents obtained from the simulations. We also compute the average fraction of infected agents using the Markovian equations through 10 realizations of a network model and compute the epidemic thresholds presented in Section 2. Fig. 2 and Fig. 3 present the results. Fig. 2, Fig. 3 indicate that the average density of infected agents obtained from the Markovian equations is close to that obtained from the agents-based simulations. Moreover, Fig. 2, Fig. 3 show that the thresholds computed by Eq. (13) are close to the epidemic thresholds derived from simulations.

Fig. 2

Average fraction of infected agents as a function of infection rate for epidemic spreading in BA networks with . The solid curves represent the solutions derived from the Markovian equations, and the points represent the results obtained from the agents-based simulations. The vertical dashed lines correspond to the epidemic thresholds computed by the Markovian equations. (a) Each node contains 10 agents. (b) Each node contains 50 agents. (c) Each node contains 100 agents.

Fig. 3

Average fraction of infected agents as a function of infection rate for epidemic spreading in WS networks with and . The solid curves represent the solutions derived from the Markovian equations, and the points represent the results obtained from the agents-based simulations. The vertical dashed lines correspond to the epidemic thresholds computed by the Markovian equations. (a) Each node contains 10 agents. (b) Each node contains 50 agents. (c) Each node contains 100 agents.

Simulation results for interconnected networks

We investigate epidemic spreading in interconnected metapopulation networks through numerical simulations. At first, we verify our theoretical analysis through Monte Carlo simulations. Specifically, we simulate the epidemic process in 3 classes of interconnected networks with two layers, namely, WS–WS networks, BA–WS networks, and BA–BA networks. Here, a BA–WS network is an interconnected network whose first layer is a BA network and whose second layer is a WS network (the interpretations of other classes of interconnected networks are similar). Given an interconnected network, we restrict to the case that two layers have the same number of nodes and all nodes contain the same number of agents. Moreover, we generate interlayer connections of an interconnected network as follows. For each node in the first layer, we randomly choose a node that has not been chosen before in the second layer and connect to . The connection is the only interlayer connection for both and . In the experiments, the number of nodes of each layer in any given interconnected network is fixed to . The simulation of epidemic spreading is described as follows. We initially randomly choose an agent as a seed in each node. In other words, every node initially contains an infected agent. Then the simulation of epidemic spreading based on the model in Section 3.1 runs until the spreading process enters into the stationary state. We generate networks and run the simulation of epidemic spreading times for each generated network. The average fraction of infected agents is returned as the final result. We also compute the average fraction of infected agents using the Markovian equations through 10 realizations of an interconnected network model and compute the epidemic thresholds using the method presented in Section 3.3.

Figs 4, 5, and 6 present the results for BA–BA networks, BA–WS networks, and WS–WS networks respectively. Figs 4, 5, and 6 indicate that the average fraction of infected agents in the stationary state computed by the Markovian equations is close to that obtained from the agents-based simulations. It is clear from the insets of Figs 4, 5, and 6 that the epidemic thresholds computed by Eq. (43) are close to the thresholds derived from simulations when . Notice that it is impossible to compute the thresholds using Eq. (43) when . For example, Fig. 4(d) does not contain vertical dashed lines.

Fig. 4

Average fraction of infected agents as a function of infection rate for epidemic spreading in BA–BA networks with . The solid curves represent the solutions derived from Markovian equations, and the points represent the results obtained from agents-based simulations. The vertical dashed lines correspond to the epidemic thresholds computed by Markovian equations. In (a), (b) and (c), . In (d), (e) and (f), and . (a) and (d) Each node contains 10 agents. (b) and (e) Each node contains 50 agents. (c) and (f) Each node contains 100 agents.

Fig. 5

Average fraction of infected agents as a function of infection rate for epidemic spreading in BA–WS networks. The first layer is a BA network with and the second layer is a WS network with and . The solid curves represent the solutions derived from the Markovian equations, and the points represent the results obtained from the agents-based simulations. The vertical dashed lines correspond to the epidemic thresholds computed by the Markovian equations. In (a), (b) and (c), . In (d), (e) and (f), and . In (g), (h) and (i), and . In (a), (d) and (g), each node contains 10 agents. In (b), (e) and (h), each node contains 50 agents. In (c), (f) and (i), each node contains 100 agents.

Fig. 6

Average fraction of infected agents as a function of infection rate for epidemic spreading in WS–WS networks with and . The solid curves represent the solutions derived from the Markovian equations, and the points represent the results obtained from the agents-based simulations. The vertical dashed lines correspond to the epidemic thresholds computed by the Markovian equations. In (a), (b) and (c), . In (d), (e) and (f), and . (a) and (d) Each node contains 10 agents. (b) and (e) Each node contains 50 agents. (c) and (f) Each node contains 100 agents.

Finally, We investigate the influence of human mobility on epidemic spreading in interconnected metapopulation networks. In reality, it is likely that a contagion initially outbreaks in a city and then spreads globally. Further, infection rates and relocation probabilities may vary in different states or countries. Thus the process of simulations is modified as follows. We assume that an infected agent in the first (second) layer infects a susceptible agent in the second (first) layer with probability () and infects a susceptible agent in the same layer with probability (). Further, we assume that each agent in the first layer decides to move to another node in the first (second) layer with probability () and each agent in the second layer decides to move to another node in the first (second) layer with probability (). Initially, we randomly choose a seed node and infect every agent in the chosen node. Then we run the agents-based simulation. In the simulation, each layer consists of nodes, and each node contains agents. To simplify the simulation, we set , , , and . Fig. 7 presents the results. Fig. 7 indicates that the average fraction of infected agents does not decrease when or increases. Thus migration or transportation processes facilitate epidemic outbreaks. We also observe that the epidemic patterns are different between Fig. 7(c) and (d) for BA–WS networks. BA networks are heterogeneous networks and WS networks are homogeneous networks. Then the result implies that the network topology of the layers in interconnected networks impacts epidemic spreading.

Fig. 7

Average fraction of infected agents as a function of infection rates and relocation probabilities for epidemic spreading in interconnected metapopulation networks. In (a), (b) and (c), , and the initial seed is in the first layer. In (d), , and the initial seed is in the second layer. (a) BA–BA networks. (b) WS–WS networks. (c) and (d) BA–WS networks.

Conclusion

Human mobility is a critical ingredient in the metapopulation epidemiological models. To deeply understand the influence of human mobility patterns on epidemic spreading, we have developed an epidemiological model that combines the radiation model with the MIR model. Here the radiation model is a superior model to predict population movement in the real world. Thus the proposed model is a more realistic metapopulation epidemiological model in comparison to the original MIR model.

We have also introduced Markovian equations using the tensorial framework to describe the dynamics of the proposed model on interconnected networks. Then we have derived the corresponding epidemic threshold by solving the Markovian equations in the stationary state. It is straightforward to generalize the presented analysis for interconnected networks to that for general multilayer networks based on the tensorial framework. We have tested the validity of the Markovian equations through comprehensive agent-based Monte Carlo simulations. The solutions derived from the Markovian equations are in good agreement with the results obtained from Monte Carlo simulations.

Finally, we have investigated the influence of human mobility on epidemic spreading in interconnected metapopulation networks through Monte Carlo simulations. The results of the simulations indicate that migration or transportation processes can facilitate epidemic propagation in interconnected metapopulation networks.

CRediT authorship contribution statement

Jinyu Huang: Conceptualization, Methodology, Experiments, Funding acquisition, Writing and editing. Chao Chen: Writing – original draft, Supervision, Software, Validation, Funding acquisition.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.