Modeling epidemic in metapopulation networks with heterogeneous diffusion rates

1.   Introduction

Metapopulation model is not only used to describe population reproduction, migration, competition and death ^[1,2,3,4], but also describe the spread of disease in real networks ^[5,6,7,8,9]. Colizza et al. have done a lot of representative results on metapopulation models with heterogeneous degree distribution, which mainly consider the spread of epidemic influenced by individual movement ^{[10,11,12,13,14,15,16]}. Juher et al. have given threshold for metapopulation epidemic model with uncorrelated network ^[17,18]. The reaction-diffusion equation is usually used to describe the propagation process, which assumed that the processes of reaction and diffusion occur simultaneously ^{[14,17,19,20]}. Furthermore, there is becoming up-front trend that concerns not only disease spreading but also human decision making; whether he is committing vaccination or not ^{[21,22,23,24,25]}.

In an SIS dynamic system, there are two kinds of processes between individual states: $I\xrightarrow{\mu}S, I+S\xrightarrow{\beta}2I.$ Here $\mu$ is the recovery rate and $\beta$ is the transmission rate across an infective contact. The influence of individuals diffusion was considered on disease transmission in papers ^[26,27]. Here the node represents a city or an airport, and the edge represents the connection between cities. For a system with $N$ nodes, which include two types of individuals $S$, $I$. The diffusion rate of infected individuals and susceptible individuals in different subpopulations are $D_{I}, D_{S}(D_{I}, D_{S}$ are all nonnegative). Thus the dynamic system is as follows:

where $\rho_{S, i}(\rho_{I, i})$ is the average density of susceptible individuals (infected individuals) in node $i$, $A_{ji} (1\leq j, i\leq N)$ is an adjacency matrix, and $k_{j}$ is the degree of node $j$. If there is a connection between node $j$ and $i$, then $A_{ji} = 1$; otherwise $A_{ji} = 0$. In the first equation of above system, the first term represents the number of susceptible individuals becoming infected, and the second term is the number of infected individuals recovering susceptible, and the third term is the number of individuals that diffuse away from the node; the last term is the number of individuals which diffuse into the node.

The data of the International Air Transport Association was analyzed in ^[26,28,29]. It is shown that there is a strong heterogeneity about the airport connectivity and traffic capacity. The different carrying capacity of route is considered in ^[11,12,30], which is reflected in the size of the traffic flow. In paper ^[11], the metapopulation epidemic model is established in the mean field, and the average traffic flow in the subpopulation with degree $k$ in per unit time is $T_{k} = k\sum_{k^{'}}p(k^{'}/k)\omega_{kk^{'}},$ where $p(k^{'}/k)$ is the conditional probability that the node with degree $k$ connects the node with degree $k^{'}$(in the uncorrelated network $p(k^{'}/k) = k^{'}p(k^{'})/\langle k\rangle)$. Here $\omega_{kk^{'}} = \omega_{0}(kk^{'})^{\nu}$ is the average weight, where $\omega_{0}$ is the coefficient of a particular system, and $\nu$ is the heterogeneity parameter of traffic flow ($0\leq\nu\leq1$). In the uncorrelated network $T_{k}$ can be wrote in the form of $T_{k} = \omega_{0}k^{1+\nu}\langle k^{1+\nu}\rangle /\langle k\rangle.$ The diffusion rate of nodes with degree $k$ is $D_{k} = T_{k}/\rho_{k}$. The effects of diffusion on the disease due to traffic driving have also been studied in ^[31].

In this paper we consider an epidemic model in metapopulation network with heterogeneous diffusion rate. The heterogeneity includes the difference of urban scale, the discrepancy of traffic conditions and so on ^{[11,32,33,34]}. We use $D_{S}^{ij}, D_{I}^{ij}$ to express the diffusion rate of susceptible and infected individuals from node $i$ to node $j$, which is more in line with the heterogeneity of metapopulation network. In fact, when the traffic of a city is more developed, the number of people flowing into and out of the city is larger. So, for the convenience of studying the problem, we assume that $D_{S}^{ij} = D_{S}^{ji} = D_{S}^{i}$. From the view of travel, infected individuals will reduce the amount of travel, thus the relation is taken as follows: $D_{I}^{i} = r D_{S}^{i}$, where $0 < r\leq 1$ is constant. In this paper, we also give the relationship among the diffusion rate, connectivity and the heterogeneity parameter of traffic flow. The results show that, if considering the heterogeneity of the degree and the heterogeneity of the traffic flow at the same time, the relationship can be obtained among the diffusion rate, connectivity and the heterogeneity of traffic flow. It can be got a more real conclusion than before.

The paper is organized as follows. The model is formulated in section 2 and the disease-free equilibrium is obtained. Some mathematical analysis are given in section 3 and the stability of the disease-free equilibria of the model is investigated. The existence and stability of the endemic equilibria of the model are studied in sections 4 and 5 respectively. In section 6, numerical simulations are illustrated.

2.   The model

Considering an SIS model with nonlimited transmission, the master equation is obtained as follows

In uncorrelated networks, $A_{ji}$ can be approximately expressed in the form $A_{ji}\simeq k_{j}k_{i}/(N\langle k\rangle)$ ^[26], where $\langle k\rangle = \sum_{i = 1}^{N}k_{i}/N$ is the average degree of the network, thus one obtains the following equations for the epidemic spread in metapopulation networks:

It is easy to see that the total density of individuals $\rho(t) = \rho_{S}(t)+\rho_{I}(t)$ remains constant and equal to $\rho_{0}$, the initial average number of individuals per cities, where $\rho_{S} = \sum_{j = 1}^{N}\rho_{S, j}/N$, $\rho_{I} = \sum_{j = 1}^{N}\rho_{I, j}/N$. In these metapopulation networks, the connectivity matrix $\Phi$ is given by

3.   The stability of disease-free equilibrium

The disease-free equilibrium of system (2.2) is $E_{0} = \left(k_{1}\rho_{0}/\langle k\rangle, \cdots, k_{N}\rho_{0}/\langle k\rangle, 0, \cdots, 0\right)$. The local stability of the disease-free equilibrium can be determined by Jacobian matrix. The Jacobian matrix of system (2.2) at the disease-free equilibrium is

where $A, B, O, C$ are matrix blocks of $N\times N$, $A = diag(D_{S}^{i})(\Phi-I_{d})$, $\Phi$ is the connectivity matrix, $I_{d}$ is the identity matrix, $B = -diag(\beta k_{i}\rho_{0}/\langle k\rangle-\mu)$, $O$ is the null matrix, and $C = diag(\beta k_{i}\rho_{0}/\langle k\rangle-\mu)+diag(D_{I}^{i})(\Phi-I_{d})$. Then, the characteristic polynomial of $J_{_{E_{_{0}}}}$ are the product of the characteristic polynomial of diagonal matrix blocks $P_{J_{_{E_{_{0}}}}}(\lambda) = P(\lambda^{A})P(\lambda^{C})$. Taking $D_{S}^{\min} = \min\{D_{S}^{i}, i = 1, 2, \cdots, N\}$ and we consider another matrix $\bar{A} = D_{S}^{\min}(\Phi-I_{d})$, which has the maximum eigenvalue $\lambda_{\max}^{\bar{A}} = 0$. Noticed that $\lambda_{\max}^{A}\leq\lambda_{\max}^{\bar{A}} = 0$, it can be proved that the eigenvalues of the matrix block $A$ are all non-positive.

By using the eigenvalue perturbation theorem in ^[35], the relationship between the eigenvalues of $C$ and $diag(\beta k_{i}\rho_{0}/\langle k\rangle-\mu)$ is given by

Therefore, the maximum eigenvalue of $J_{_{E_{_{0}}}}$ is $\lambda_{\max} = \max\{0, \lambda_{\max}^{C}\}$. If $\lambda_{\max} > 0$, the disease-free equilibrium of system (2.2) is unstable, and the system has an endemic equilibrium. Thus we have the following theorem:

Theorem 3.1. The sufficient condition of the disease-free equilibrium to be unstable is

Theorem 3.2. It can be seen from Eq (3.2), for fixed $\mu$, even if $\beta$ is sufficient small, if $\rho_{0}$ is large enough, the disease-free equilibrium is also unstable. While for fixed $\rho_{0}$, the diffusion rates $D_{I}^{i}$ will affect the spread of the disease.

As mentioned in the first section, the traffic flow is the physical quantity in the mean field. Here, we assume that all nodes with the same degree are one class ^[14,15]. In this paper $T_{k_{i}}$ represent traffic flow with degree $k_{i}$ in node $i$. According to the previous discussion in the first section, the diffusion rate of node $i$ can be expressed as

where $\rho_{i} = \frac{k_{i}\rho_{0}}{\langle k\rangle}$ is the average density of node $i$. From Eq (3.3), it is not difficult to see that $D_{i}$ is a function of $k_{i}$. Similarly, we can get the diffusion rates $D_{I}^{i}$, $D_{S}^{i}$. Based on $D_{i} = \frac{T_{I, k_{i}}+T_{S, k_{i}}}{\rho_{I, i}+\rho_{S, i}}$, the relationship can be got about $D_{i}, D_{I}^{i}, D_{S}^{i}:$

Substituting $\rho_{S, i} = \rho_{i}-\rho_{I, i}$ into the Eq (3.4), it can be simplified

(where $h_{i} = \frac{\rho_{I, i}}{\rho_{i}}$). Substituting (3.3) into (3.5), we can get

It can be seen that $D_{I}^{i}$ is a function of $k_{i}$ and parameter $\nu$.

Theorem 3.3. Let $D_{I}^{i} = r D_{S}^{i}, (0 < r \leq 1)$, then

When $r = 1$, then $D_{I}^{i} = D_{S}^{i} = \frac{\langle k^{1+\nu}\rangle \omega_{0}k_{i}^{\nu}}{\rho_{0}}$, which is satisfying (3.4).

Theorem 3.4. The unstable condition of the disease-free equilibrium is $\beta > \beta_{c}$, where

Here $k_{\max}$ is the maximum degree of all nodes in the metapopulation network. When $\nu = 0$, $\beta_{c} = \frac{\langle k\rangle (\mu\rho_{0}+\langle k\rangle\omega_{0})}{\rho_{0}^{2} k_{\max}}$; And when $\nu = 1$, there is $\beta_{c} = \frac{\langle k\rangle (\mu\rho_{0}+\langle k^{2}\rangle\omega_{0}k_{\max})}{\rho_{0}^{2}k_{\max}}$.

Comparing the above two thresholds, it is shown that the epidemic threshold increases with the traffic flow heterogeneity parameter increase.

4.   The existence and uniqueness of endemic equilibrium

The endemic equilibrium $E^{*} = (\rho_{S, 1}^{*}, \rho_{S, 2}^{*}, \cdots, \rho_{S, N}^{*}, \rho_{I, 1}^{*}, \rho_{I, 2}^{*}, \cdots, \rho_{I, N}^{*})$ of system (2.2) is satisfied

where $\rho_{S}^{*} = \frac{1}{N}\sum_{i = 1}^{N}\rho_{S, i}^{*}$, $\rho_{I}^{*} = \frac{1}{N}\sum_{i = 1}^{N}\rho_{I, i}^{*}$, and $\rho_{S, i}^{*}$ can be expressed

Putting Eq (4.1) into the second equation of system (2.2) and solving a quadratic equation, then we get

where $\delta = \beta\frac{k_{i}}{\langle k\rangle}\rho_{I}^{*}(\frac{D_{I}^{i}}{D_{S}^{i}}-1), \; \theta = 4\beta\frac{k_{i}}{\langle k\rangle}\rho_{I}^{*}D_{I}^{i}$. The negative root of equation has been taken out.

Now, $\rho_{I, i}^{*} > 0$ need to be proven. Taking the summation over $i$ and multipling $\frac{1}{N}$ for (4.2), it can be obtained as follows

where $0\leq\rho_{I}^{*}\leq\rho_{0}$. We define function

where $\delta^{'} = \beta\frac{k_{i}}{\langle k\rangle}\rho_{I}(\frac{D_{I}^{i}}{D_{S}^{i}}-1), \; \theta^{'} = 4\beta\frac{k_{i}}{\langle k\rangle}\rho_{I}D_{I}^{i}.$

When $D_{I}^{i} = rD_{S}^{i}\; (0 < r \leq 1)$, we get

where $\delta^{'} = \beta\frac{k_{i}}{\langle k\rangle}\rho_{I}(r-1), \; \theta^{'} = 4\beta\frac{k_{i}}{\langle k\rangle}\rho_{I}D_{I}^{i}$. Especially, when $r = 1$, then

$F(\rho_{0}) < 0$ is always satisfied.

When $\beta = \beta_{c}$, we can obtain $F(0) = 0$;

$F^{''}(\rho_{I}) < 0$. So, there is a $\rho_{I}^{*}$ that makes $F(\rho_{I}^{*}) = 0.$

When $\beta > \beta_{c}$, we can get $F(0) > 0$. Thus, there is also a $\rho_{I}^{*}$ that makes $F(\rho_{I}^{*}) = 0.$

Theorem 4.1. There is a unique endemic equilibrium $E^{*} = (\rho_{S, 1}^{*}, \rho_{S, 2}^{*}, \cdots, \rho_{S, N}^{*}, \rho_{I, 1}^{*}, \\ \rho_{I, 2}^{*}, \cdots, \rho_{I, N}^{*})$ of system (2.2), if the infection rate is larger than epidemic threshold $\beta \geq \beta_{c}$.

Theorem 4.2. In fact it can be proved that the existence of $\rho_{I}^{*}$ that make (4.3) satisfied, when $0 < r < 1$. Therefore, the existence and uniqueness of the solution are independent of the relation about $D_{I}^{i}, D_{S}^{i}.$

5.   The stability of endemic equilibrium

The Jacobian matrix of system (2.2) at the endemic equilibrium is

When $r = 1$, the above matrix can be transformed into a new matrix as follows:

Thus the characteristic polynomial of $J^{'}_{E^{*}}$ can be given by $P_{J^{'}_{E^{*}}}(\lambda^{'}) = P(\lambda^{'}_{A^{'}})P(\lambda^{'}_{C^{'}}).$

For the matrix $A^{'} = D_{S}^{i}(\Phi-I_{d})$, $\lambda^{'}_{A^{'}}\leq0$ is also satisfied according to section 3. $C^{'} = diag[\beta(\rho_{S, i}^{*}-\rho_{I, i}^{*})-\mu]+D_{I}^{i}(\Phi-I_{d})$ is equivalent to the sum of a diagonal matrix and a perturbation matrix. So, using the general formula of secular equation ^[35], we can get

where

are the eigenvalues of $diag[\beta(\rho_{S, i}^{*}-\rho_{I, i}^{*})-\mu-D_{I}^{i}]$, and $\lambda^{'}_{C^{'}}$ is an eigenvalue of $C^{'}$. Even $\alpha_{i} < \lambda^{'}_{C^{'}}$, but $\lambda^{'}_{C^{'}} < 0$ can not be judged. Then, we prove $\lambda^{'}_{C^{'}} < 0$ is satisfied.

Making right side of the second formula of system (2.2) is 0, we can obtain $\rho_{I, i}^{*}(\beta \rho_{S, i}^{*}-\mu) = D_{I}^{i}(\rho_{I, i}^{*}-\frac{k_{i}}{\langle k\rangle})\rho_{I}^{*}$, which is simplified to

Summing over $i$ and dividing by $N$ for Eq (5.2), we obtain $\langle k\rangle N = \sum_{i = 1}^{N}\frac{D_{I}^{i}k_{i}}{(\mu+D_{I}^{i})-\beta\rho_{S, i}^{*}}.$ So, it can be judged that

Eq (5.3) can be rewritten as

From Eqs (5.1) and (5.4), it can be inferred that $\lambda^{'}_{C^{'}} < 0$. So the eigenvalues of $J_{E^{*}}$ is negative. According to the above statement, we have the following theorem.

Theorem 5.1. The endemic equilibrium $E^{*}$ of system (2.2) is always locally asymptotically stable, if it exists.

6.   Numerical simulation and sensitivity analysis

In the above section, we get the equilibrium existence and stability of system (2.2) through theoretical analysis, and get the epidemic threshold. The following is the numerical simulations.

6.1.   On scale-free network

We simulate the spread of disease among $N$ interacting nodes, here $N = 100$. It is found that there are two main factors that affect the diffusion rate and epidemic threshold: One is the connectivity of nodes, and another is the heterogeneity parameter of the traffic flow. On scale-free network, the adjacency $A$ of order $N\times N$ is generated randomly ($A$ is a symmetric matrix with the diagonal elements are all 0, the other elements are 0 or 1). Then, the degree of node $i$ is the sum of the elements of all the rows in line $i$ of $A_{ij}$. On scale-free network, the degree distribution obeys power-law distribution $p(k)\sim k^{-\gamma}$. Therefore, the effect of the parameter $\nu$ for the diffusion rate and epidemic threshold should be explored by some sensitivity analysis about $\nu$. In the real word network, the traffic flow have a certain heterogeneity. It can be seen from Figure 1a, $D_{I}^{i}$ increase with the increase of $\nu$, and the heterogeneity of the diffusion rate also increase. At the same time, the threshold $\beta_{c}$ is also increase in Figure 1b. That is to say, the parameter $\nu$ increases the epidemic threshold and suppresses the outbreak of diseases. In the case of $\beta > \beta_{c}$, the endemic equilibrium of the system is locally asymptotically stable. When $\nu = 0.6$, the number of susceptible individuals tends to be stable; At the same time, the number of infected individuals gradually tends to be stable, thus the disease is prevalent. It is shown that the more developed of city traffic is, the more people will be infected (Figure 1a). And for the whole system, the change of average density of overall $S, I$ population is similarity (Figure 1b).

For one city, the heterogeneity of traffic flow will also affect the infection rate. For the cities with larger degree, the greater the $\nu$ is, the more the infected individuals are. That is to say, the heterogeneity of traffic flow has an influence on the spread of the disease. For the cities with smaller degree, with the increase of $\nu$, the disease will go extinct (Figure 2). In general, the traffic heterogeneity of cities is about $\nu \approx 0.5$ ^[11], which can inhibit the prevalence of the disease.

6.2.   On small-world network

In this section, the spread of disease on the small-world network is studied. The influence of $\nu$ for the diffusion rate and epidemic threshold are shown in Figure 3. In Figure 4, it is the change of individuals (susceptible and infected) and the overall average density over time. From Figure 5, it is shown that the influence of $\nu$ for the infected individuals and the infection rate. In Figure 7 the comparison of the epidemic between on the small-world and scale-free networks is given.

It can be seen that for the appropriate $\nu$, the greater the degree of the city is, the greater the diffusion rate (Figure 4). With the increase of $\nu$, $D_{I}^{i}$ and the epidemic threshold $\beta_{c}$ all increase. That is also to say, the parameter $\nu$ increases the epidemic threshold and suppresses outbreak of the disease (Figure 5). In the case of $\beta > \beta_{c}$, the endemic equilibrium of the system is locally asymptotically stable. If the epidemic threshold is reached, the prevalence of the disease will be promoted with increase of parameter $\nu$ (Figure 6). From Figure 7, the system is locally asymptotically stable at the disease-free equilibrium on the small-world network, while the disease is spreading on the scale-free network. It can be seen that the disease is more prevalent on the scale-free network than on the small-world network at the same parameters. The results show that the heterogeneity of traffic flow has a greater impact on the disease, and the epidemic is easier to be controlled on the small-world network.

7.   Conclusion

In this paper, an SIS model in metapopulation networks with heterogeneous diffusion rates is established. According to the qualitative analysis of dynamics, the existence and stability of the disease-free equilibrium are analyzed, and the epidemic threshold is obtained. When the epidemic threshold is reached, the disease-free equilibrium is unstable, and the system has an endemic equilibrium. It is proved that it is locally asymptotically stable, if the endemic equilibrium is existing.

Due to there is a big difference in the traffic level and population density in each city, the diffusion rate of each city also have heterogeneity. In this paper, we consider the spread of disease with heterogeneity of the degree and heterogeneity of diffusion rate, which is more suitable for the real world. Based on this study, the relationship among the diffusion rate, connectivity and the heterogeneity of traffic flow are given. There are two conclusions: On the one hand, for the larger degree nodes, the number of infected individuals is higher than the smaller one. When the epidemic threshold is reached, the more developed of city is, the higher the diffusion rate is, which result in the large number of individuals enter from other nodes. Therefore, this increases the spread of disease; on the other hand, from the view of disease control, heterogeneity of the traffic flow by increase the epidemic threshold, can improve the ability to control the disease. Finally, numerical simulations are compared on the scale-free network and the small-world network.

In this paper, because the eigenvalues of the Jacobian matrixes at the equilibria could not be directly calculated, we estimated them by using other methods. In addition, we consider that the diffusion rate from different cities is the same, and the diffusion rate, which is proportion to the traffic flow, will be discussed in the future.

Acknowledgments

The authors would like to thank the reviewers for their helpful comments and valuable suggestions, and the support of the National Sciences Foundation of China (11571324, 61403393), the Fund for Shanxi "1331KIRT", Shanxi Scholarship Council of China and the Scientific Activities of Selected Returned Overseas Professionals in Shanxi Province.

Conflict of interest

The authors declare there is no conflict of interest.