Moment closure of infectious diseases model on heterogeneous metapopulation network

The global transmission of infectious diseases poses huge threats to human. Traditional heterogeneous mean-field models on metapopulation networks ignore the heterogeneity of individuals who are in different disease states in subpopulations with the same degree, resulting in inaccuracy in predicting the spread of disease. In this paper, we take heterogeneity of susceptible and infectious individuals in subpopulations with the same degree into account, and propose a deterministic unclosed general model according to Markov process on metapopulation networks to curve the global transmission of diseases precisely. Then we make the general model closed by putting forward two common assumptions: a two-dimensional constant distribution and a two-dimensional log-normal distribution, where the former is equivalent to the heterogeneous mean-field model, and the latter is a system of weighted ordinary differential equations. Further we make a stability analysis for two closed models and illustrate the results by numerical simulations. Next, we conduct a series of numerical simulations and stochastic simulations. Results indicate that our general model extends and optimizes the mean-field model. Finally, we investigate the impacts of total mobility rate on disease transmission and find that timely and comprehensive travel restriction in the early stage is an effective prevention and control of infectious diseases.

In order to study the global spread of infectious diseases, researchers applied a metapopulation model to infectious diseases [7,8]. The concept of metapopulation was put forward by Levins [9,10] for the first time to investigate the processes of local extinction, recolonization and regional persistence of populations, which means "the population of populations". With the development of network transmission dynamics, the metapopulation model has been successfully applied to understand the transmission dynamics of mobility. The model is composed of a heterogeneous network of subpopulations, connected by mobility processes. Individuals in each subpopulation stay one of the two states: one is susceptible; the other is infectious. Individuals can move from a subpopulation to another along links of network. Once there exist infectious individuals in a subpopulation, it becomes infectious. For infectious diseases, a susceptible individual can be infected at rate β and become infectious, while an infectious individual may recover at rate μ and be susceptible spatially structured populations with well-defined social units (on a scale of countries, regions, cities, small as towns, villages, families) connected through individuals' mobility on networks, named "metapopulation networks" (see Fig. 1).
In general, we consider susceptible-infectious-susceptible (SIS) transmission process on a metapopulation network with V nodes, label the nodes with the elements of integer set A = {1, . . . , V }, and denote by S k j and I k j the number of susceptible and infectious individuals in a node which gets label j and degree k, respectively. The total individuals of node whose label is j and degree is k is N k j and N k j = S k j + I k j . Then the sum of susceptible and infectious individuals of all nodes with degree k are j∈A S k j , j∈A I k j , respectively. Let V * k denote the number of nodes with degree k. According to heterogeneous mean-field (HMF) theory, assuming that subpopulations with the same degree are statistical equivalence, that is to say, the average number of individuals in all nodes with degree k is Further, the average numbers of susceptible and infectious individuals in nodes with degree k are respectively. Based on the HMF assumption, Colizza and Vespignani [11] proposed some models to address the transmission of diseases on heterogeneous metapopulation networks under two different mobility patterns. They assumed that the process of diseases spreading was in the first, and mobility process next. For the mobility process, at each time step, an individual in a subpopulation with degree k, whether susceptible or infectious, travels to other subpopulations at total mobility rate δ. Along a link, an individual in a subpopulation whose degree is k travels to another subpopulation of degree k at rate d kk , and δ = k k P(k |k)d kk . When the mobility rate depends on traffic (passengers along the link), where w kk represents the average traffic (passengers along the link node with degree k and node with degree k per day) on the link between two subpopulations with degrees k and k , and behaves as here w 0 and θ are positive constants, T k is the total average traffic (passengers move out of node) of subpopulations with degree k (per day), and Then Colizza and Vespignani established a model as follows: As is known in the case of complex networks, the notation P(k |k) represents the conditional probability that any given edge departing from a node of degree k is pointing to a node of degree k . Moreover, this work provided a good tool to calculate a global invasion threshold.
For the results reported so far, the vast majority of researches are based on HMF theory [11]. According to (1.1), under the assumption that size of population in nodes with the same degree k stays the same, the equations are equivalent to the case where the number of susceptible and infectious individuals in nodes whose degree are k are the same with each other, respectively. This conclusion holds when V * k = 1. For V * k = 2, with loss of generality, these two nodes get label 1 and 2, respectively. The number of susceptible and infectious individuals in node label 1 are S k 1 and I k 1 . And there are S k 2 susceptible individuals and I k 2 infectious individuals in a node whose label is 2. From (1.6), we have (1.7) Substituting (1.1) into (1.7), and multiplying by 2 on the two sides of the above equations lead to To make the equations hold, either S k 1 = S k 2 or I k 1 = I k 2 , according to the assumption S k 1 + I k 1 = N k 1 = N k 2 = S k 2 + I k 2 , so both S k 1 = S k 2 and I k 1 = I k 2 are satisfied. For the case of V * k > 2, for simplicity, label these V * k nodes with integers 1, . . . , V * k , respectively, according to Eq. (1.6) and Eq. (1.1), we have We use the assumption that all nodes with the same degree k have the same number of individuals, i.e., N k i = N k j , i, j = 1, . . . , V * k . Replacing S k j by N k j -I k j , and after some algebra, we have From what has been discussed above, under the assumption that the number of individuals of all nodes with degree k keeps consistent, the number of susceptible individuals and infectious individuals of all nodes with degree k, respectively, being the same is a necessary condition for Eq. (1.6). Once I k i = I k j for some i and j, Eq. (1.6) does not hold. It is obvious that the HMF assumption in metapopulation network ignores heterogeneities in subpopulations with the same degree. Under this assumption, the total individuals, the numbers of susceptible individuals and infectious individuals in subpopulations with the same degree are always the same at any time, respectively, which is hard to achieve. Even though the size of a population in nodes with the same degree is identical, there may be a great difference for individuals in different states among nodes in the course of infectious disease spreading, resulting in error for forecasting infectious disease spreading. Motivated by reducing even eliminating this error, we take heterogeneities in subpopulations with the same degree into account, define the number of susceptible and infectious individuals in subpopulations with the same degree as a two-dimensional random variable and we propose a deterministic SIS model according to a continuous-time Markov chain (CTMC) on heterogeneous metapopulation networks to curve the global transmission of infections more precisely. The results show that our model extends and optimizes the HMF model in [11].
The paper is organized as follows. In Sect. 2, we firstly present necessary assumptions and propose a deterministic general model based on CTMC. Next, two moment closure models are given based on a two-dimensional constant distribution and a two-dimensional log-normal distribution in Sects. 3 and 4, respectively. Furthermore, in Sect. 5 we conduct numerical simulations on models and stochastic realizations. Conclusions and a discussion are given in Sect. 6. For simplicity, the definitions of abbreviations presented in the article are given in Table 1.

The SIS infectious disease model based on CTMC
In this section, following the idea of stochastic process, we derive a deterministic SIS model on metapopulation network. The following assumptions will be used throughout the paper.
(H1) Global spread of infectious diseases is relatively fast timescales such that changes in demographics (e.g., births, aging, deaths) is negligible. (H2) Network is connected and the minimal and maximal degrees are k min and k max , respectively. (H3) The memory of individuals for the origin is not taken into account. (H4) For the Markov process, at each time step t, the following four events occur simultaneously. Two-dimensional random variable with the numbers of susceptible and infectious individuals in a subpopulation with degree k.
The number of nodes (subpopulations) with degree k, in which susceptible and infectious individuals are s and i at time t, respectively, s, i ∈ 0, 1, . . . , M. p k s,i (t) Probability that the numbers of susceptible and infectious individuals in a subpopulation with degree k are s and i at time t, respectively, that is, The transition probability from state (s, i) to state (s + s, i + i) in the interval t in a homogeneous Markov process.  For a subpopulation with degree k, in emigration process, each individual in a subpopulation leaves this subpopulation at rate δ to its neighbor subpopulations; while in an immigration process, each individual in neighbor subpopulation with degree k travels to it at rate d k k .

Model derivation
Before using CTMC to derive the deterministic model, in Table 2 we present the necessary notations for model parameters and model formation. To clarify the derivation of our models, the probability function is elucidated in what follows. The probability function with respect to susceptible and infectious individuals in a subpopulation with degree k is where s = 0, 1, . . . , M, i = 0, 1, . . . , Ms, and s,i p k s,i (t) = 1. If (s, i) lies outside of this range, the probability is assumed to be zero.

The stochastic process of two-dimensional random variable (S k , I k )
For this stochastic process, firstly, we derive the forward Kolmogorov differential equation about two-dimensional random variable (S k , I k ).

The forward Kolmogorov differential equation about
Note that, for a homogeneous Markov process, the transition probability does not depend on the starting time t, thus it is expressed as p k s+ s,i+ i ( t). In the paper, we assume the Markov process is homogeneous. Therefore, the transition probabilities about (S k , I k ) satisfy in which the time step t must be chosen sufficiently small Hence, we obtain the forward Kolmogorov equation as follows: Subtracting p k s,i (t) on both sides, dividing by t and letting t → 0 give rise to the forward Kolmogorov differential equation, where si k represents the average total number of all susceptible individuals contacting with infectious individuals in a subpopulation with degree k.
Applying the same approach, we give the deterministic equation about i k as follows: We therefore obtain the following deterministic model: For the sake of convenience, we call system (2.7a)-(2.7b) the general model. Obviously, the general model (2.7a)-(2.7b) is not closed. In order to make it closed, we need to derive the rate of changes of si k . Multiplying Eq. (2.2) by the product si, summing over s and i and simplifying yield here s 2 i k = s,i s 2 ip k s,i , and si 2 k = s,i si 2 p k s,i . The derivation above brings new parameters into play, i 2 k , the average number of all infectious individuals contacting with infectious individuals in a subpopulation with degree k. In the same way, the rate of change of So far, the model is still not closed because of the appearance of the third-order cumulants: s 2 i k and si 2 k . Notice that calculation of third-order cumulants may bring new variables into play: fourth-order cumulants, and so on. Approaching the third-order cumulants by first-and second-order cumulants becomes necessary. It is clear that with different twodimensional quasi-stationary distributions of p k s,i , the closed models vary. We will present two main moment closure equations in Sects. 3 and 4.

Model analysis
Although the general model (2.7a)-(2.7b) is not closed, there exists a common property. Notice that the average total number of individuals in a subpopulation with degree k sat-isfies N k = s k + i k .
Summing (2.7a) and (2.7b), we have In the paper, we overlook the degree correlations, that is to say, P(k |k) = k P(k )/ k . Let whereN , which is a positive constant, represents the average number of individuals in a subpopulation. It is obvious that Eqs. (2.9) have a unique globally asymptotically stable equilibrium N * k = k 1+θ k 1+θ N , k = k min , . . . , k max .

Moment closure based on two-dimensional constant distribution
We first present a closed model based on a two-dimensional constant distribution.

Derivation of moment closure model based on two-dimensional constant distribution
For a fixed degree k, when the quasi-stationary of p k s,i approaches a two-dimensional constant distribution, we have si k = s k i k , so the closed model is The assumption of two-dimensional constant distribution is equivalent to HMF assumption. Under the supposition that mobility process and disease transmission process occur simultaneously, our model (3.1a)-(3.1b) and the models in [11] are equivalent.
In order to study the asymptotic behavior of model, we consider the limiting system of (3.1a)-(3.1b), For simplicity, we call (3.2) the moment closure model I.

Global basic reproduction number
Firstly, we deduce global basic reproduction number by the approach in van den Driessche and Watmough [28]. Note that model (3.2) admits a unique disease-free equilibrium (DFE) E 0 = ( n 0, . . . , 0), n = k maxk min + 1. In the DFE E 0 , the rate of appearance of new infections F and the rate of transfer of individuals out of the compartments V are given by here i, j ∈ 1, . . . , n. Thus Using the next-generation matrix [28], the global basic reproduction number is R 0 = ρ(FV -1 ), the spectral radius of the matrix FV -1 .

Global stability of disease-free equilibrium
Since i k ∈ [0, N * k ] for k = k min , . . . , k max , we study system (3.2)

Lemma 3.1 The set n is positively invariant for system (3.2).
Theorem 3.1 For system (3.2), if R 0 < 1, DFE E 0 is globally attractive in n .
Proof To complete the proof, it is sufficient to show that lim t→+∞ i k = 0, k = k min , . . . , k max .
For system (3.2) with i k ≤ N * k , we can obtain the inequality group Define an auxiliary linear system The coefficient matrix of the above system is F -V . When R 0 = ρ(FV -1 ) < 1, all eigenvalues of F -V lie in the left half plane. Thus each non-negative solution of (3.3) satisfies This implies that the zero solution of (3.4) is globally asymptotically stable. By comparison, each non-negative solution of (3.2) satisfies Accordingly, DFE E 0 of system (3.2) is globally attractive.

Global stability of endemic equilibrium
Next, we analyze the existence and global stability of endemic equilibrium.
Proof We will use the theory of cooperate system in Corollary 3.2 in [29] to prove the existence and global stability of endemic equilibrium.
In fact, let f : + n → n be defined by the right-hand side of (3.2), f = (f k min , . . . , f k max ). Clearly f is continuously differentiable, f (0) = 0, f i (y) ≥ 0 for all y (= ( i k min , . . . , i k max )) ∈ + n with y i = 0 and ∂f i /∂y j ≥ 0, i = j for y ∈ + n . So f is cooperative. Clearly Dy = (∂f i /∂y j ) k min ≤i,j≤k max is irreducible for every y ∈ + n . Note that, for ∀α ∈ (0, 1) and y i > 0, Thus f is strong sublinear on + n . By Lemma 2 and Corollary 3.2 in [29], we conclude that system (3.2) admits a unique endemic equilibrium E * = ( i * k min , . . . , i * k max ) which is globally asymptotically stable.

Moment closure based on two-dimensional log-normal distribution
Following Keeling [30] we assume that the quasi-stationary distribution is approximately two-dimensional log-normal. The calculation of all cumulants of orders greater than first is performed by considering all possible pairwise combinations of the elements and multiplying by the appropriate moment. In this case, where ξ k is the multiplicative covariance between infectious and susceptible individuals within a subpopulation with degree k. Taking derivatives by t at both ends of the second equality at the same time, transposition, substitution and simplification, leads to This brings in two other parameters,V s k andV i k , the variance in susceptibles and infectious within a subpopulation with degree k. In a similar way, we obtain From the above derivation, we obtain a moment closure model based on the log-normal distribution Accordingly, we have the following limiting system: For simplicity, we denote it as moment closure model II.

Numerical results
In this section, we conduct extensive numerical simulations, stochastic simulations, and compare them. Numerical simulations of models are obtained by the fourth-order Runge-Kutta algorithm on MATLAB. We report stochastic simulation results from Monte Carlo simulations in a variety of different realizations on FORTRAN. Metapopulation networks are generated with the uncorrelated scale-free network model [31,32] with V ranging from 100 to 1000 following the power-law degree distribution P(k) ∼ k -γ , 2 < γ ≤ 3 with minimum degree k min = 2 and maximum degree k max ≤ V 1/2 . According to Ref. [33], w 0 = 1 and θ = 0.5. Simulation results are based on averaging over at least 50 realizations for initial conditions and network structures.

Stability of closed models
First of all, focusing on stability of equilibria, we develop a series of numerical simulations about models (3.2) and (4.3a)-(4.3d) to discuss stability of systems.

Stability of equilibria of moment closure model I
We present numerical integration of model (3.2) under two cases: R 0 < 1 and R 0 > 1 to investigate how the spread of infectious diseases depend on the threshold R 0 and whether the system stable or not. In Fig. 2, when β = 2e -5 , R 0 is 0.9177 < 1 (see Fig. 2(a)), while R 0 becomes 9.1729 > 1 for β = 2e -4 (see Fig. 2(b)). We find that when R 0 < 1 the number of infectious individuals all approach zero under different initial conditions. Besides, the fraction of infectious individuals approaches 0.394 while R 0 > 1.

Stability of equilibrium of moment closure model II
We further study system (4.3a)-(4.3d). As a system of weighted ordinary differential equations, it is unreasonable to perform i k = 0. For simplicity, we assume that every node has infectious seeds at initial moments. In Fig. 3, the density of infectious individuals approaching a constant under different initial conditions, approximately 0.96, means that there exists an endemic equilibrium for system (4.3a)-(4.3d).

Time courses of the density of infectious individuals
According to the forward Kolmogorov differential equation (2.2), we put forward model (2.6) to address the spread of infectious diseases on metapopulation networks. Further, to make the model (2.6) closed, we derive models (3.2) and (4.3a)-(4.3d). In Fig. 4

The effects of total mobility rate on disease transmission
Finally, we discuss the impacts of total mobility rate on disease transmission and perform five magnitudes of total mobility rate. In Fig. 5, when δ = 0.1, diseases globally spread and rapidly achieve steady state. Then with the decrease of total mobility rate, the transmission slows down. When δ reduces to 0.00001, diseases hardly spread in a short time. Instead, in a long time, t = 3000 or longer, infectious diseases will spread to the entire network as in Fig. 6. It is worth noticing that there is no impact of mobility rate on the steady of  In detail, β = 5e -3 and other parameters are the same with Fig. 2. The values are obtained by averaging over 50 stochastic realizations whole metapopulation network. So travel restriction in the early stage is an effective prevention and control of infectious diseases, and restriction must be timely and thoroughly, prohibiting travel over a period of time, for example.

Conclusions and discussion
Following CTMC, we have derived a deterministic unclosed general model to portray disease transmission on metapopulation networks in which the heterogeneity of susceptible and infectious individuals in subpopulations with the same degree was considered. Then we closed the general model under the assumption of a two-dimensional constant distribution and two-dimensional log-normal distribution, respectively. And the existence and global stability of each of feasible equilibria of the system, which is based on twodimensional constant distribution, have been proved mathematically and illustrated by numerical simulations. It has been shown by simulations that the general model we derived generalizes and optimizes the HMF model. In the study of the effects of total mobility rate on infection spread, we have seen that total mobility rate has a huge impact on disease transmission speed, but it has no effect on the steady state of infections. It is worth noting that even for a relatively small total mobility rate, diseases will spread globally as long as time is long enough. Therefore, timely and comprehensive travel restrictions are critical for disease control and prevention. Although we have generalized and optimized the HMF model, there are still some problems in this paper to be further solved.
(1) The deterministic general model extends and optimizes the HMF model in [11], however, numerical simulations did not fit perfectly the stochastic simulations before the transmission was up to steady. Hence, our model needs to be further improved. (2) We conducted numerical simulations for the moment closure II, and we found this system has an endemic equilibrium. It is interesting to verify the existence and global stability of the endemic equilibrium mathematically. (3) On the edge weights selection for metapopulation network, we supposed that weights were only dependent on the degrees of subpopulations. However, a real situation may be more complex, and weights may be the comprehensive outcome of populations, degrees and distances. For the total mobility rate, it is more appropriate to vary with the states of individuals.