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.

Introduction

In recent years, global transmission of infectious diseases, such as severe acute respiratory syndromes (SARS) [1], influenza A (H1N1) flu [2], avian influenza [3], Middle East respiratory syndrome coronavirus (MERS-CoV) [4], Ebola virus disease [5], and zika [6], has been threatening human beings. Great attention has been paid to the effects of human mobility on disease transmission.

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 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).

Figure 1

A metapopulation network model of SIS infections with individuals’ 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

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 , and denote by and 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 and . Then the sum of susceptible and infectious individuals of all nodes with degree k are respectively. Let 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 at rate , and . When the mobility rate depends on traffic (passengers along the link), where represents the average traffic (passengers along the link node with degree k and node with degree per day) on the link between two subpopulations with degrees k and , and behaves as here and θ are positive constants, 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 represents the conditional probability that any given edge departing from a node of degree k is pointing to a node of degree . Moreover, this work provided a good tool to calculate a global invasion threshold.

Applications of HMF theory have launched a variety of researches: factors including network structures [12-15], human mobility patterns [11, 16–18], human behaviors [19-22], human contact patterns [14, 23], travel restrictions [24], heterogeneous dwelling time in subpopulations [18], pathogen competition [25], and so forth, have been investigated. It has been shown that substrate network structures play an essential role in the spatial spread of infectious diseases [12, 13, 15]. In real-world networks, human mobility patterns vary in a very complicated way, e.g., large-scale air-travel (typical case [11]), small-scale recurrent visits of subpopulations (or commuting flows) [16–18, 26, 27], etc. Safety-driven people’s behavioral responses to the infectious diseases have also been found to be able to accelerate the infectious disease spread [19-22], contrary to willingness. With regard to human contact patterns, heterogeneous mixing [23] and subpopulations with network topology have been investigated.

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 . For , 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 and . And there are susceptible individuals and infectious individuals in a node whose label is 2. From (1.6), we have Substituting (1.1) into (1.7), and multiplying by 2 on the two sides of the above equations lead to yielding To make the equations hold, either or , according to the assumption , so both and are satisfied. For the case of , for simplicity, label these nodes with integers , 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., . Replacing by , and after some algebra, we have that is, Hence , accordingly . 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 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.

Table 1

Abbreviations in the article

Abbreviations	Full names
SARS	Severe Acute Respiratory Syndromes
MERS-CoV	Middle East respiratory syndrome coronavirus
SIS	susceptible-infectious-susceptible
HMF	heterogeneous mean-field
CTMC	continuous-time Markov chain
DFE	disease-free equilibrium

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.

Global spread of infectious diseases is relatively fast timescales such that changes in demographics (e.g., births, aging, deaths) is negligible.

Network is connected and the minimal and maximal degrees are and , respectively.

The memory of individuals for the origin is not taken into account.

For the Markov process, at each time step Δt, the following four events occur simultaneously.

Disease transmission process: each susceptible individual is infected at transmission rate β by infectious individuals.

Recovery process: every infectious individual recovers at rate μ.

Mobility process consists of two processes: emigration and immigration. 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 travels to it at rate .

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.

Table 2

Parameters description

Major parameters	Description
M	The maximum number of individuals among all subpopulations.
(S^k,I^k)	Two-dimensional random variable with the numbers of susceptible and infectious individuals in a subpopulation with degree k.
V_k(s,i)(t)	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.
ps,ik(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, ps,ik(t)=P{(S^k,I^k)(t)=(s,i),S^k,I^k∈0,1,…,M}.
ps+Δs,i+Δik(Δt)	The transition probability from state (s,i) to state (s + Δs,i + Δi) in the interval Δt in a homogeneous Markov process.
〈f(s,i)〉_k	The expectation of the function f(s,i), 〈f(s,i)〉_k=∑_s,if(s,i)ps,ik. Specially, 〈s〉_k and 〈i〉_k are expectations of susceptible individuals and infectious individuals, respectively.

The probability function with respect to susceptible and infectious individuals in a subpopulation with degree k is where , , and . If lies outside of this range, the probability is assumed to be zero.

The stochastic process of two-dimensional random variable

For this stochastic process, firstly, we derive the forward Kolmogorov differential equation about two-dimensional random variable .

The forward Kolmogorov differential equation about

Define the transition probability about from state to state in the time interval Δt as Note that, for a homogeneous Markov process, the transition probability does not depend on the starting time t, thus it is expressed as . In the paper, we assume the Markov process is homogeneous. Therefore, the transition probabilities about satisfy in which the time step Δt must be chosen sufficiently small

Hence, we obtain the forward Kolmogorov equation as follows: Subtracting on both sides, dividing by Δt and letting give rise to the forward Kolmogorov differential equation,

Now we calculate the expectations for susceptible and infectious individuals: and in subpopulations with degree k.

The deterministic equations of expectations

Further, we apply the obtained forward Kolmogorov differential equation (2.2) to derive the deterministic equations of , i.e., the rates of changes of the expectations of , respectively.

Now, we give the deterministic equation of . Multiplying Eq. (2.2) by s and summing over s and i, yield Simplifying Eq. (2.3), we obtain Substituting , and into Eq. (2.4), we get where 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 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 . Multiplying Eq. (2.2) by the product si, summing over s and i and simplifying yield here , and . The derivation above brings new parameters into play, , 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 reads So far, the model is still not closed because of the appearance of the third-order cumulants: and . 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 two-dimensional quasi-stationary distributions of , 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 satisfies Summing (2.7a) and (2.7b), we have In the paper, we overlook the degree correlations, that is to say, . Let Equation (2.8) becomes where N̄, 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 .

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 approaches a two-dimensional constant distribution, we have , so the closed model is

Remark 3.1

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) . In the DFE , the rate of appearance of new infections F and the rate of transfer of individuals out of the compartments V are given by and where here .

Thus Using the next-generation matrix [28], the global basic reproduction number is , the spectral radius of the matrix .

Global stability of disease-free equilibrium

Since for , we study system (3.2) in .

Lemma 3.1

The set is positively invariant for system (3.2).

Theorem 3.1

For system (3.2), if , DFE is globally attractive in .

Proof

To complete the proof, it is sufficient to show that

For system (3.2) with , we can obtain the inequality group

Define an auxiliary linear system

The coefficient matrix of the above system is . When , all eigenvalues of 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 of system (3.2) is globally attractive. □

Global stability of endemic equilibrium

Next, we analyze the existence and global stability of endemic equilibrium.

Theorem 3.2

If , system (3.2) admits a unique endemic equilibrium which is globally asymptotically stable with respect to any initial value :=.

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 be defined by the right-hand side of (3.2), . Clearly f is continuously differentiable, , for all y with and , for . So f is cooperative. Clearly is irreducible for every .

Note that, for and , Thus f is strong sublinear on . By Lemma 2 and Corollary 3.2 in [29], we conclude that system (3.2) admits a unique endemic equilibrium 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 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, and , the variance in susceptibles and infectious within a subpopulation with degree k. In a similar way, we obtain and

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 , with minimum degree and maximum degree . According to Ref. [33], and . 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: and to investigate how the spread of infectious diseases depend on the threshold and whether the system stable or not. In Fig. 2, when , is (see Fig. 2(a)), while becomes for (see Fig. 2(b)). We find that when the number of infectious individuals all approach zero under different initial conditions. Besides, the fraction of infectious individuals approaches 0.394 while .

Figure 2

Evolution of the moment closure model I with the varying initial infectious individuals. Here we take , , , , . Three lines represent different initial infectious seeds in a node chosen randomly. (a) ; (b)

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 . 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).

Figure 3

Evolution of the moment closure model II with the varying initial infectious individuals. Here we take . Three lines represent different initial seeds every node. The other parameters are the same with Fig. 2

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, we compare these models with each other based on numerical simulations, and make a comparison with stochastic simulations at the same time. Red line, green line and blue line are calculated from the general model (2.6) coupling with (2.2), moment closure model I (3.2) and moment closure model II (4.3a)–(4.3d), respectively, while dots are averaging over 100 stochastic realizations. It is showed that the steady state of moment closure model II is higher than the others. Hence the assumption of a two-dimensional log-normal distribution is inappropriate for metapopulation networks. In the steady state, general model I fits perfectly with stochastic simulations, and moment closure model I takes the second place. For an early stage of infectious diseases, they overestimate the spreading of diseases on different levels. Overall, the general model extends and optimizes the HMF model.

Figure 4

Comparison of models and stochastic simulations. In all cases, seed 10 infectious individuals in subpopulation with the maximum degree in the initial time. Here we take , , , , ,

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 , 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, 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 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.

Figure 5

The impacts of total mobility rate on disease transmission. Five magnitudes of total mobility rate are performed. From bottom to top, increase one magnitude in turn. The left panel shows the fraction of individual infectious evolving over time, while the right panel is the proportion of infectious subpopulations. In detail, and other parameters are the same with Fig. 2. The values are obtained by averaging over 50 stochastic realizations

Figure 6

The density of infectious individuals for . The other parameters are the same with Fig. 5

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 two-dimensional 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.

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.

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.

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.