# Exploring the Epidemic Spreading in a Multilayer Metapopulation Network by considering Individuals’ Periodic Travelling

Authors: Dun Han; Qi Shao; Dandan Li
Journal: Complexity
 (2020)

---

## Abstract

The convenience of transportation brings the diversity of individuals’ travelling modes; in this paper, we present an improved epidemic diffusion model in a multilayer metapopulation network. Firstly, we construct the metapopulation network with different travelling ways, and then, the epidemic spreading threshold is calculated by means of the mean-field method. Taking the periodicity of individuals’ travelling into account, we further explore the epidemic diffusion model with individuals’ periodic travelling and deduce the epidemic spreading threshold using the Perron–Frobenius theorem. Our results show that if all individuals in each area decide to move, the epidemic threshold can be effectively raised while each individual chooses an unbiased region to arrive. In addition, with the increase of individuals’ mobility rate or regional heterogeneous infection coefficient, the fluctuation range of the density of infected becomes larger, while the fluctuation period is almost unchanged. However, the change of individuals’ periodic motion could cause the change of the fluctuation period of infected density. We try to provide a new perspective for the research of metapopulation.


---

## Body

## 1. Introduction

The study of epidemic diffusion based on dynamic model has a history of more than 200 years. It can be traced back to 1760 when Bernoulli, a Dutch mathematician and physicist, researched the vaccinated smallpox in order to prolong people’s life span [1]. Since the small-world [2] and scale-free networks [3] are presented, many scholars have studied the epidemic spreading using dynamic analysis method based on complex network model [4–9] and obtained a lot of better prediction results for the spread of infectious diseases. The use of mathematical models in the study of epidemic dynamics and their mechanisms and their prevention and control provided valuable insights and contributed to draw up the global picture of epidemics occurrence, patterns, and management [10–13]. However, considering the impact of mathematical models in epidemiology, it is important to build them on realistic assumptions [14–16].The metapopulation model has been developed in ecological modeling [17], which consists of spatially separated subpopulations of a species [18], and individuals can migrate between subpopulations [19]. A metapopulation model is a useful mathematical framework to understand the dynamics of epidemic spreading in spatially separated subpopulations, corresponding to countries, cities, or local communities, which are connected by pathways for human migration and population mobility [20]. Recently, metapopulation dynamics on networks, which incorporates mobility over the subpopulations, local interaction at the subpopulations, and a complex network structure, has gained great research attention [19, 21–23]. Ni and Weng studied the impact of travel patterns on epidemic dynamics in heterogeneous spatial metapopulation networks. Their results showed that the occurrence probability of global outbreaks is significantly dependent on the characteristic travel distance, the characteristic waiting time, and the memory effects of human travel [24]. Wang et al. developed a theoretical framework to investigate the safety-information-driven human mobility with metapopulation epidemic dynamics. Their findings uncovered that the critical threshold depends on mobility in such a way that personal execution of safe moves unexpectedly promotes the global spread of a disease, while unsafe moves counterintuitively cause a locally, relatively small outbreak size [25]. Halvor et al. introduced a generalized reaction-diffusion approach to epidemic spreading in a metapopulation modeled as a complex network. They found that the heterogeneous conditions cause the epidemic threshold to be a nontrivial function of the reaction rates, the network’s topology, and the crossover population size that separates “village dynamics” from “city dynamics” [26]. Naoki analyzed the effects of diffusion rates on epidemic spreads in metapopulation networks. The results showed that diffusion can suppress epidemics contrasts with that for diffusive SIS dynamics and its variants when individuals are fixed at nodes on static networks [27].At present, although there are many studies on the diffusion of infectious diseases in the metapopulation networks, few studies have considered the individuals’ travelling patterns [28–31]. Especially, the convenience of transportation brings the diversity of individual travel. For example, when people in the region A decide to travel to the region B, they could choose a train, a car, or a plane. In a result, we could regard different travelling modes as different transportation networks. However, the multilayer system is a very effective tool to analyze this case. In this paper, according to the difference of individuals’ travelling modes, we propose an epidemic spreading model on the multilayer metapopulation network and then calculate the epidemic spreading threshold by means of mean field method. The theoretical analysis and numerical simulation results shows that if an individual arrives to other subpopulations without bias, the diffusion of diseases can be effectively suppressed. Moreover, empirical observations in biological and social systems show that individuals’ travelling behaviors are regular; for example, the individuals in region A may periodically travel to region B. Inspired by the above fact, we analyze the epidemic diffusion model by considering individuals’ periodic movement. The results show that the change of individuals’ periodic movement leads to the change of fluctuation cycle of infected density.The rest of this paper is organized as follows. In Section2, preliminaries and model formulation for the epidemic spreading in the multilayer metapopulation network is proposed. In Section 3, the multilayer metapopulation network with individuals’ periodic movement is explored. In Section 4, simulation examples are given to illustrate the theoretical analysis. The conclusion is given in Section 5.


## 2. Epidemic Spreading in the Multilayer Metapopulation Network

### 2.1. Preliminaries and Model Formulations

For two real-valued column vectorsu=u1,u2,…,uNT, v=v1,v2,…,vNT∈RN, we say u≺¯v if ui≤vi for all i∈1,2,…,N. For a vector u=u1,u2,…,uNT, if ui>0,i=1,2,…,n, we note u>0. The classical SIS epidemic spreading model includes two types of individuals: the susceptible individual S and the infected individual I. Each susceptible individual becomes infected with probability λ when it comes into contact with an infected individual, and the infected nodes are cured and become susceptible again with probability β. The infection dynamics can be identified by the following set of reaction equations: S+I⟶λ2I, I⟶βS. In SIR epidemic diffusion model, a class of recovered individuals R is added. The corresponding epidemic diffusion model can be described as S+I⟶λ2I, I⟶βR. For the listed two kinds of epidemic spreading systems, in the mixed homogeneous populations, it is easy to calculate the epidemic spreading threshold: λC=λ/β, when λC=λ/β>1, infectious diseases will always exist, but when λC=λ/β>1, infectious diseases will disappear.In our model, the total population is divided inton subpopulations, each of which is represented by a node of the metapopulation network. We assume that each node i is initially occupied by ni individuals, so the total number of individuals in the whole system is N=∑i=1nni. Individuals belonging to node i decide whether or not to move at each time step. Without loss of generality, we assume that individuals belonging to region i move to other regions and return to region i after a time step. The probability of moving is p, the probability of remaining is 1−p. If an individual has L≥0 travelling modes (walking, car, train, ship, airplane, etc.), each travelling mode could be regarded as a layer of metapopulation network (see Figure 1). For an individual in subpopulation i and deciding to move to region j, he will choose a travelling mode l∈0,L with probability Pijl=Wijl/∑l=1L∑j=1nWijl.Figure 1
Individuals’ travelling map with three travelling modes in two subpopulations. Initially, there are six individuals belonging to regioni and region j, respectively. The individuals with the same color own the same travelling mode.


### 2.2. SIS Epidemic Spreading Model

When all individuals in the metapopulation network have finished their moving, the infected densityρi belonging to region i (including the primitive inhabitants in region i and those who move to other regions) varies as follows:(1)ρit+1=1−βρit+1−ρitϕi,where β represents recovery rate and the value of ϕi indicates the infected probability who belonging to region i, and it is specifically expressed as follows:(2)ϕi=1−pψi+p∑l=1L∑j=1nPijlψj,where ψi indicates the probability that the susceptible individuals staying in the region i are infected, specifically as follows:(3)ψi=1−∏j=1n1−λjiρjtnj⟶i,where λji denotes the probability that the susceptible individuals belong to region i and infected by the infected individuals of region j; nj⟶i denotes the total number of individuals who move from region j to region i, which is described as follows:(4)nj⟶i=1−pniδij+p∑l=1LPjilnj,where δij represents the Dirac function, if i=j, δij=1, and i≠j, δij=0.We incorporate the comprehensive factors, including the relationship (cooperation, affection, friendship, etc.), the interaction frequency, and the impact of individual environment (for example, one may easy to be in the same room or classroom with the infected one, etc.); therefore, we assume that individuals are easy to be infected by infected individuals in their areas. In this paper, the general expression of infection probabilityλji could be presented as follows:(5)λji=fλ0,α,i,j,where λ0 represents the underlying infection rate, α∈0,1 indicates the regional heterogeneous infection coefficient caused by regional differences. The expression of λji we adopted in this paper is as follows:(6)λji=λ01−α1−δji.Equation (6) shows that λji≤λjj, i.e., the infected individuals are more likely to infect the susceptible individuals in the same region.Theorem 1.
The epidemic spreading threshold in the metapopulation network considering multiple travelling modes is as follows:λ0c=β/ΛmaxM, where ΛmaxM is the spectral radius of matrix M, M=Mijn×n, where Mik=1−pniδik+p∑l=1LPkilnk1−p1−α1−δki+∑l=1L∑j=1nPijlp1−pnjδjk+p∑l=1LPkjlnk1−α1−δkj.Proof.
According to equations (1) and (2), we could get that at the equilibrium point:(7)βρi∗=1−ρi∗1−pψi+p∑l=1L∑j=1nPijlψj.
In order to calculate the epidemic threshold, we linearize equation (7) and neglect the higher order terms of rhoi∗, we obtain(8)βρi∗=1−pψi+p∑l=1L∑j=1nPijlψj.
Based on equation (3), we can get(9)ψi=1−∏j=1n1−λjiρj∗nj⟶i=∑j=1nnj⟶iλjiρj∗.
Combining equations (4), (6), (8), and (9), we obtain(10)βλρi∗=∑k=1n1−pniδik+p∑l=1LPkilnk1−p1−α1−δki+∑l=1L∑j=1nPijlp1−pnjδjk+p∑l=1LPkjlnk1−α1−δkjρk∗.
Letρ∗=ρ1∗,ρ2∗,…,ρn∗T,M=Mijn×n, where(11)Mik=1−pniδik+p∑l=1LPkilnk1−p1−α1−δki+∑l=1L∑j=1nPijlp1−pnjδjk+p∑l=1LPkjlnk1−α1−δkj.
Further, we can get(12)βλ0ρ∗=Mρ∗.
So we can calculate the epidemic spreading threshold as follows:(13)λ0c=βΛmaxM.Corollary 1.
Ifp=0, that is, all individuals stay in their respective subpopulation without moving, we get if λ/β<1/maxni, i=1,2,…,n, infectious diseases do not break out. However, once λ/β<1/maxni, infectious diseases will break out.Proof.
Letx=x1,x2,…,xn, y=y1,y2,…,yn, constructing the Lagrangian Functions Lx,y,λ1,λ2=∑i=1nxiyi+λ1∑i=1nxi−1+λ2∑i=1nyi−1, let ∂L/∂xi=0, ∂L/∂yi=0, ∂L/∂λ1=0, ∂L/∂λ2=0, 1≤i≤n, we can get xi=yi=1/n,i=1,2,…,n.Corollary 2.
Ifp=1 and the number of individuals in each region is the same at the initial time, i.e., ni≡nc,i=1,2,…,n, the epidemic spreading threshold is λ0c=β/nc1−a.Proof.
According to equation (11), we can get(14)Mik=∑l=1L∑j=1nPijl∑l=1LPkjlnk1−α1−δkj.
LetPij=∑lPijl, so we can get(15)Mik=∑l=1L∑j=1nPijl1−α1−δkjnkPkj=∑j=1nPijPkjnk1−α1−δkj.
If the number of individuals at the beginning of each region is the same, then we obtain(16)Mik=nc∑j=1nPijPkj1−α1−δkj.
Further, we can get(17)Mik=nc1−α∑j=1nPijPkj+αPikPkk.
As∑jPij=1, according to Lemma 1, we can get(18)∑k=1nMik≥∑k=1nnc1−α∑j=1nPijPkj≥nc1−a,i=1,2,…,n.
If and only ifPij=1/n, 1≤i,j≤n, the second equal sign was established. Because ΛmaxM≥min∑j=1nMij, 1≤i≤n; therefore, the transmission threshold of infectious diseases is λ0c=β/nc1−a.It can be concluded from Corollary2 that in a multilayer metapopulation network, if p=1 and when the individuals move to other regions with equal probability, the epidemic spreading threshold can be effectively increased, i.e., if each individual chooses an unbiased region to arrive, the diffusion of infectious diseases can be effectively suppressed.Remark 1.
In this part, for the convenience of analysis, we replace∑l=1LPijl by Pij, 1≤i,j≤n. However, in the actual situation, we need to consider different travelling modes and set different value of Pij, 1≤i,j≤n.


### 2.3. SIR Epidemic Spreading Model

In theSIR epidemic spreading model, a new state R is included, so that we can make corresponding improvements to equation (1):(19)ρit+1=1−βρit+1−ρit−ritϕi,rit+1=rit+βρit,where rit is the density of recovered individuals in region i at the time t. In order to calculate the epidemic spreading threshold, we linearize equation (19) and neglect the high order of ρi∗ and ri∗, then we obtain(20)βρi∗=ϕi.According to the threshold calculation method inSIS model, we can get a similar result; that is, the epidemic spreading threshold is λ0c=β/ΛmaxM.


## 2.1. Preliminaries and Model Formulations

For two real-valued column vectorsu=u1,u2,…,uNT, v=v1,v2,…,vNT∈RN, we say u≺¯v if ui≤vi for all i∈1,2,…,N. For a vector u=u1,u2,…,uNT, if ui>0,i=1,2,…,n, we note u>0. The classical SIS epidemic spreading model includes two types of individuals: the susceptible individual S and the infected individual I. Each susceptible individual becomes infected with probability λ when it comes into contact with an infected individual, and the infected nodes are cured and become susceptible again with probability β. The infection dynamics can be identified by the following set of reaction equations: S+I⟶λ2I, I⟶βS. In SIR epidemic diffusion model, a class of recovered individuals R is added. The corresponding epidemic diffusion model can be described as S+I⟶λ2I, I⟶βR. For the listed two kinds of epidemic spreading systems, in the mixed homogeneous populations, it is easy to calculate the epidemic spreading threshold: λC=λ/β, when λC=λ/β>1, infectious diseases will always exist, but when λC=λ/β>1, infectious diseases will disappear.In our model, the total population is divided inton subpopulations, each of which is represented by a node of the metapopulation network. We assume that each node i is initially occupied by ni individuals, so the total number of individuals in the whole system is N=∑i=1nni. Individuals belonging to node i decide whether or not to move at each time step. Without loss of generality, we assume that individuals belonging to region i move to other regions and return to region i after a time step. The probability of moving is p, the probability of remaining is 1−p. If an individual has L≥0 travelling modes (walking, car, train, ship, airplane, etc.), each travelling mode could be regarded as a layer of metapopulation network (see Figure 1). For an individual in subpopulation i and deciding to move to region j, he will choose a travelling mode l∈0,L with probability Pijl=Wijl/∑l=1L∑j=1nWijl.Figure 1
Individuals’ travelling map with three travelling modes in two subpopulations. Initially, there are six individuals belonging to regioni and region j, respectively. The individuals with the same color own the same travelling mode.


## 2.2. SIS Epidemic Spreading Model

When all individuals in the metapopulation network have finished their moving, the infected densityρi belonging to region i (including the primitive inhabitants in region i and those who move to other regions) varies as follows:(1)ρit+1=1−βρit+1−ρitϕi,where β represents recovery rate and the value of ϕi indicates the infected probability who belonging to region i, and it is specifically expressed as follows:(2)ϕi=1−pψi+p∑l=1L∑j=1nPijlψj,where ψi indicates the probability that the susceptible individuals staying in the region i are infected, specifically as follows:(3)ψi=1−∏j=1n1−λjiρjtnj⟶i,where λji denotes the probability that the susceptible individuals belong to region i and infected by the infected individuals of region j; nj⟶i denotes the total number of individuals who move from region j to region i, which is described as follows:(4)nj⟶i=1−pniδij+p∑l=1LPjilnj,where δij represents the Dirac function, if i=j, δij=1, and i≠j, δij=0.We incorporate the comprehensive factors, including the relationship (cooperation, affection, friendship, etc.), the interaction frequency, and the impact of individual environment (for example, one may easy to be in the same room or classroom with the infected one, etc.); therefore, we assume that individuals are easy to be infected by infected individuals in their areas. In this paper, the general expression of infection probabilityλji could be presented as follows:(5)λji=fλ0,α,i,j,where λ0 represents the underlying infection rate, α∈0,1 indicates the regional heterogeneous infection coefficient caused by regional differences. The expression of λji we adopted in this paper is as follows:(6)λji=λ01−α1−δji.Equation (6) shows that λji≤λjj, i.e., the infected individuals are more likely to infect the susceptible individuals in the same region.Theorem 1.
The epidemic spreading threshold in the metapopulation network considering multiple travelling modes is as follows:λ0c=β/ΛmaxM, where ΛmaxM is the spectral radius of matrix M, M=Mijn×n, where Mik=1−pniδik+p∑l=1LPkilnk1−p1−α1−δki+∑l=1L∑j=1nPijlp1−pnjδjk+p∑l=1LPkjlnk1−α1−δkj.Proof.
According to equations (1) and (2), we could get that at the equilibrium point:(7)βρi∗=1−ρi∗1−pψi+p∑l=1L∑j=1nPijlψj.
In order to calculate the epidemic threshold, we linearize equation (7) and neglect the higher order terms of rhoi∗, we obtain(8)βρi∗=1−pψi+p∑l=1L∑j=1nPijlψj.
Based on equation (3), we can get(9)ψi=1−∏j=1n1−λjiρj∗nj⟶i=∑j=1nnj⟶iλjiρj∗.
Combining equations (4), (6), (8), and (9), we obtain(10)βλρi∗=∑k=1n1−pniδik+p∑l=1LPkilnk1−p1−α1−δki+∑l=1L∑j=1nPijlp1−pnjδjk+p∑l=1LPkjlnk1−α1−δkjρk∗.
Letρ∗=ρ1∗,ρ2∗,…,ρn∗T,M=Mijn×n, where(11)Mik=1−pniδik+p∑l=1LPkilnk1−p1−α1−δki+∑l=1L∑j=1nPijlp1−pnjδjk+p∑l=1LPkjlnk1−α1−δkj.
Further, we can get(12)βλ0ρ∗=Mρ∗.
So we can calculate the epidemic spreading threshold as follows:(13)λ0c=βΛmaxM.Corollary 1.
Ifp=0, that is, all individuals stay in their respective subpopulation without moving, we get if λ/β<1/maxni, i=1,2,…,n, infectious diseases do not break out. However, once λ/β<1/maxni, infectious diseases will break out.Proof.
Letx=x1,x2,…,xn, y=y1,y2,…,yn, constructing the Lagrangian Functions Lx,y,λ1,λ2=∑i=1nxiyi+λ1∑i=1nxi−1+λ2∑i=1nyi−1, let ∂L/∂xi=0, ∂L/∂yi=0, ∂L/∂λ1=0, ∂L/∂λ2=0, 1≤i≤n, we can get xi=yi=1/n,i=1,2,…,n.Corollary 2.
Ifp=1 and the number of individuals in each region is the same at the initial time, i.e., ni≡nc,i=1,2,…,n, the epidemic spreading threshold is λ0c=β/nc1−a.Proof.
According to equation (11), we can get(14)Mik=∑l=1L∑j=1nPijl∑l=1LPkjlnk1−α1−δkj.
LetPij=∑lPijl, so we can get(15)Mik=∑l=1L∑j=1nPijl1−α1−δkjnkPkj=∑j=1nPijPkjnk1−α1−δkj.
If the number of individuals at the beginning of each region is the same, then we obtain(16)Mik=nc∑j=1nPijPkj1−α1−δkj.
Further, we can get(17)Mik=nc1−α∑j=1nPijPkj+αPikPkk.
As∑jPij=1, according to Lemma 1, we can get(18)∑k=1nMik≥∑k=1nnc1−α∑j=1nPijPkj≥nc1−a,i=1,2,…,n.
If and only ifPij=1/n, 1≤i,j≤n, the second equal sign was established. Because ΛmaxM≥min∑j=1nMij, 1≤i≤n; therefore, the transmission threshold of infectious diseases is λ0c=β/nc1−a.It can be concluded from Corollary2 that in a multilayer metapopulation network, if p=1 and when the individuals move to other regions with equal probability, the epidemic spreading threshold can be effectively increased, i.e., if each individual chooses an unbiased region to arrive, the diffusion of infectious diseases can be effectively suppressed.Remark 1.
In this part, for the convenience of analysis, we replace∑l=1LPijl by Pij, 1≤i,j≤n. However, in the actual situation, we need to consider different travelling modes and set different value of Pij, 1≤i,j≤n.


## 2.3. SIR Epidemic Spreading Model

In theSIR epidemic spreading model, a new state R is included, so that we can make corresponding improvements to equation (1):(19)ρit+1=1−βρit+1−ρit−ritϕi,rit+1=rit+βρit,where rit is the density of recovered individuals in region i at the time t. In order to calculate the epidemic spreading threshold, we linearize equation (19) and neglect the high order of ρi∗ and ri∗, then we obtain(20)βρi∗=ϕi.According to the threshold calculation method inSIS model, we can get a similar result; that is, the epidemic spreading threshold is λ0c=β/ΛmaxM.


## 3. Epidemic Spreading in Multilayer Metapopulation Network with Individuals’ Periodic Travelling

In Section2, we analyze the epidemic spreading model in a multilayer metapopulation network; we next expand the model by considering individuals’ periodic travelling. In general, individual travelling is periodic, i.e., individuals may arrive at other regions regularly. For simplicity, we assume that all the individuals have the same travelling period. Therefore, the following changes need be modified to equations (1)–(4):(21)ρit+1=1−βρit+1−ρitϕi∗,where ϕi∗=1−p∗ψi+p∗∑l=1L∑j=1nPijlψj, ψi=1−∏j=1n1−λjiρjtnj⟶i, nj⟶i=1−p∗niδij+p∗∑l=1LPjilnj, p∗=pHsinπt/ts. H∗ denotes the Heaviside function, and Hx=1,x>00,x≤0, ts is the half period of interpopulation communication. So when the half period ts⟶+∞, 0<t<ts, p∗=p, at this point, there will always be interindividual flow between subpopulations. However, if ts=1, and t∈N+, there will be no inter-individual flow between subpopulations.Theorem 2.
(See Perron-Frobenius theorem [32]). If M=Mijn×n is an n×n nonnegative and irreducible matrix, then one of its eigenvalues is positive and greater than or equal to (in absolute value) all other eigenvalues. Such an eigenvalue is called the “dominant eigenvalue” or Perron-Frobenius eigenvalue of the matrix; there is a positive eigenvector corresponding to that eigenvalue; and ρM is equal to the dominant eigenvalue of the matrix.Lemma 1.
The matrix of strong connection graph satisfies Perron-Frobenius theorem. That is, for a strongly connected graph matrixM, there is and only one eigenvalue λ=ρM>0. The corresponding eigenvector u>0.Theorem 3.
For system (21), the epidemic spreading threshold is λ0c=β/ΛmaxM˜, where ΛmaxM˜=argtmaxρM˜, ρ∗ represents the spectral radius of a matrix, M˜=M˜ijn×n, M˜ij=1−pHsinπt/ts1−α1−δkink⟶i+pHsinπt/ts∑j=1nPij1−α1−δkink⟶i.Proof.
According to equation (21), it can be obtained that at the disease-free equilibrium point:(22)βρi1−ρi=1−pHsinπt/ts∑j=1nλ01−α1−δjinj⟶iρjt+pHsinπt/ts∑l=1L∑j=1nPijl∑k=1nλ01−α1−δkink⟶iρkt.
By derivingλ0 at both ends of equation (22), one can obtain(23)dρi/dλ01−ρi−ρidρi/dλ01−ρi2=1−pHsinπtts∑j=1n1−α1−δjinj⟶iρjt+λ01−α1−δjinj⟶idρjdλ0=pHsinπtts∑l=1L∑j=1nPijl∑k=1n1−α1−δkink⟶iρkt+λ01−α1−δkink⟶idρkdλ0.
Because at the point of epidemic spreading thresholdλ0c, there must be ρiλ0=λ0c=0, ∂ρi/∂λ0c>0, so when λ0=λ0c, we can simplify equation (23):(24)βdρidλ0=1−pHsinπtts∑j=1nλ01−α1−δjinj⟶idρjdλ0+pHsinπtts∑l=1L∑j=1nPijl∑k=1nλ01−α1−δkink⟶idρkdλ0.
Letdρi/dλ0=ξi, Pij=∑lPijl. According to simplification of equation (24), we obtain(25)ξi=1−pHsinπtts∑k=1nλ0β1−α1−δkink⟶iξk+pHsinπtts∑k=1n∑j=1nλ0βPij1−α1−δkink⟶iξk.
LetM˜=M˜ijn×n, M˜ij=1−pHsinπt/ts1−α1−δkink⟶i+pHsinπt/ts∑j=1nPij1−α1−δkink⟶i. According to equation (25), we can obtain(26)ξi=λ0β∑kM˜ikξk.
BecauseM˜ik>0, and the multilayer metapopulation network is strongly interconnected, let ΛmaxM˜=argtmaxρM˜, ρ∗ represents the spectral radius of the matrix. According to Theorem 2 and Corollary 2, the epidemic spreading threshold is as follows:(27)λ0c=βΛmaxM˜.Corollary 3.
When the half periodts⟶+∞, i.e., there are always interindividual flows between subpopulations, the epidemic spreading threshold is λ0c=β/ΛmaxM.Proof.
Because whents⟶+∞, equation (21) can be equivalent to equation (1), so it can be proved.


## 4. Numerical Simulations

Letft=1/n∑i=1nρitni denote the infected density at time t, and f∞=1/n∑i=1nρi∞ni denote the final infected density. We take n=10 regions. The number of people in each region follows a uniform distribution, i.e., ni∼U100,500, i=1,2,…,n, each individual has L=3 travelling modes. At each time, the individual randomly chooses one mode to travel. The probability of individuals in region i move to neighbouring region j by means of trips l=1,2,3 is Pijl∼U0,1 and ∑j∑lPijl=1. In order to eliminate the effect of randomness on the results, the following experimental results are averaged over 100 times.Figure2 depicts the change of infected density ft in the whole metapopulation network over time. It can be seen from Figure 2 that when λ0 is larger enough, the value of ft increases gradually. However, while the infection rate λ0 is smaller, the value of ft decreases gradually until no infected individual exists. Besides, although the interindividual flow exists in the metapopulation network, the value of ft does not show periodic oscillation.Figure 2
The infected densityft changes with time t. With the increase of the infection rate λ0, the value of ft gradually increases. Parameters are set to be p=0.5, β=0.2, α=0.1.We adopt a small infection rateλ0. On one hand, the number of individuals in the whole metapopulation is huge; on the other hand, the contact between individuals is frequent. So even at such a small infection rate, there will still be a lot of infected individuals. Figure 3 depicts how the infected density changes with the infection rate λ0. The theoretical epidemic spreading threshold λ0c=β/ΛmaxM can be calculated by equation (13). With the increase of infection rate λ0, the density of infected individuals in the whole system increases gradually, while with the increase of recovery rate β, the density of infected persons in the whole system decreases gradually. In addition, based on the theoretical epidemic spreading threshold λ0c=β/ΛmaxM, it can be concluded that for the same ΛmaxM, the ratio of epidemic spreading threshold under different parameters is λ0c1/λ0c2=β1/β2. This conclusion is proved by the numerical simulation results in Figure 3.Figure 3
The final infected densityf∞ changes with λ0. With the increase of λ0, the final infected density f∞ increases gradually. Parameters are set to be p=0.5, α=0.1.Figure4 depicts how the final infected density f∞ changes with the regional heterogeneous infection coefficient. One can observe from Figure 4 that the density of the infected individuals decreases gradually with the increase of and p. If the probability of infection in the original residence is greater than that in other regions, i.e., λij≤λii, 1≤i,j≤n, moving to other regions could effectively reduce the infected risk.Figure 4
The final infected densityf∞ changes with the regional heterogeneous infection coefficient. With the increase of α, the final density of infected f∞ in the whole system decreases gradually. Parameters are set to be β=0.1, λ=0.001.Figure5 presents the evolution of the infected density ft over time under different half period ts. According to Figure 5, we can see that in the multilayer metapopulation network with individuals’ periodic travelling, the infected density ft also changes periodically with time. In addition, the larger the individuals’ travelling period 2ts, the greater the fluctuation period of the density ft. This interesting result reveals the diffusion of recurrent infectious diseases to a certain extent.Figure 5
The evolution of infected densityft over time under different half period ts. The parameters are set as follows: p=0.5, α=0.5, β=0.1. (a) The infection rate λ0=0.001; (b) the infection rate λ0=0.002.
(a)
(b)Figure6 depicts the evolution of the infected density ft over time at different mobility rates p. By comparing Figures 6(a) and 6(b), we find the fluctuation range of the infected density increase with p and ts, whereas even if the values of p and ts change, the maximum infected density is almost the same. Therefore, under different individuals’ mobility modes and travelling period, the peak of infected density generally does not change.Figure 6
The evolution of the infected densityft with time under different value of p. We find the fluctuation range of the infected density increase with p and ts. The parameters are set as follows: α=0.5, β=0.1, λ0=0.001. (a) The half period ts=10; (b) the half period ts=30.
(a)
(b)Figure7 illustrates the change of the value of ft over time under different regional heterogeneous infection coefficient α. By comparing and analyzing Figures 7(a) and 7(b), it can be found that with the increase of the regional heterogeneous infection coefficient α, the fluctuation range of the infected density ft becomes larger and larger. In addition, when the infection rate is high λ0=0.002, it is interesting to find that the fluctuation period of the infected density ft is almost the same under different α.Figure 7
The infected densityft evolves with time under different regional heterogeneous infection coefficient α. It is interesting to find that the fluctuation period of the infected density ft is almost the same under different α. The parameters are set to be ts=50, p=0.5, β=0.1. (a) The infection rate α0=0.001; (b) the infection rate λ0=0.002.
(a)
(b)


## 5. Conclusion

The convenience of transportation brings the diversity of individual travel; that means we may not describe the individuals’ mobility just by using a single layer network. However, the multilayer system and metapopulation model are the very effective tools to model it. Some researchers use airport traffic data, train traffic data, or other public transport data to simulate the spread of infectious diseases in the metapopulation network. However, some data are still unavailable. In particular, if the data of infected people cannot be obtained, it will lead to inaccurate results, so the validity and accuracy of these data need to be further improved. At present, it is more difficult to obtain real and accurate data of different travel modes. In this paper, we propose a multilayer metapopulation network model to study the epidemic spreading process. In addition, we consider the periodicity of individual commuting, exploring the epidemic spreading dynamic behavior in the periodic metapopulation network. The epidemic spreading threshold is theoretically calculated based on the Perron-Frobenius theorem. Further, we perform a large number of numerical simulations, and the simulation results are consistent with the theoretical deduction. The results show that if each individual arrives in an unbiased regions, the diffusion of infectious diseases can be effectively suppressed. Moreover, with the increase of mobility rate or regional heterogeneous infection coefficient, the fluctuation range of infected density becomes larger and larger, while the fluctuation period is almost unchanged. However, the change of individual travelling period leads to the change of fluctuation period of infected density. At present, the researches for the epidemic spreading in the metapopulation network mostly consider the single-layer network or the individual’s behavior as nonperiodic movement, and our study attempts to provides a theoretical and analytical method for the study of epidemic spreading in the multilayer metapopulation network with individual periodic movement.



---

*Source: 6782018-2020-04-21.xml*