IMPACT OF HETEROGENEITY ON THE DYNAMICS OF AN SEIR EPIDEMIC MODEL

An SEIR epidemic model with an arbitrarily distributed exposed stage is revisited to study the impact of heterogeneity on the spread of infectious diseases. The heterogeneity may come from age or behavior and disease stages, resulting in multi-group and multi-stage models, respectively. For each model, Lyapunov functionals are used to show that the basic reproduction number R0 gives a sharp threshold. If R0 ≤ 1, then the disease-free equilibrium is globally asymptotically stable and the disease dies out from all groups or stages. If R0 > 1, then the disease persists in all groups or stages, and the endemic equilibrium is globally asymptotically stable.


1.
Introduction. Heterogeneity exists in many aspects of disease transmission processes [1,2,12,33], for example, heterogeneous spatial distribution of host populations, heterogeneous susceptibility among age groups, heterogeneous social behavior among groups for sexually transmitted diseases, and multi-hosts for many diseases such as West Nile virus and Avian influenza. Different types of heterogeneous epidemic models and attempts to understand the impact of heterogeneity on the disease transmission have appeared in the literature; see, for example, [19,22,33,37]. It has been found that heterogeneity sometimes does not alter the dynamical structure of epidemic models, while sometimes it produces more complicated dynamical behavior than found in homogeneous models. For example, if the basic reproduction number R 0 completely determines the dynamics of homogeneous epidemic models in which the disease eventually either dies out or persists at a positive level, then the associated heterogeneous models might have the same dichotomy: there exists either a globally attracting disease-free equilibrium or a globally attracting endemic equilibrium. On the other hand, if the original homogeneous model has rich dynamical behavior such as oscillations, then the associated heterogeneous models have much richer dynamics such as a switch of oscillations from one group to the other [41].
In Section 2, we formulate a basic susceptible-exposed-infectious-removed (SEIR) epidemic model in which an arbitrarily distributed exposed stage is assumed. Using this basic model as a building block, we consider two different heterogeneous models: multi-group models for heterogeneity in a host population (Section 3) and multi-stage models for heterogeneity of infectious individuals (Section 4). For both 394 ZHISHENG SHUAI AND P. VAN DEN DRIESSCHE models, we establish their global dynamics and demonstrate that the heterogeneity does not alter the dynamical structure of the SEIR model with an arbitrarily distributed exposed stage.
Our proof for the global stability of endemic equilibria utilizes global Lyapunov functionals that are motivated by the work in [20,21,34,35] and the graph-theoretic approach for large-scale systems developed in [10,11,29]. This graph-theoretic approach, based on Kirchhoff's Matrix Tree Theorem, has previously only been applied to one study of delay epidemic models, namely, the multi-group model in [30]. We demonstrate that, for the two types of heterogeneous delay epidemic models considered, this graph-theoretic approach can be successfully applied by choosing an appropriate weight matrix. 2. Basic model. Since our main goal is to investigate the impact of heterogeneity on dynamics of epidemic models, we start with a relatively simple model and then build on this to formulate and analyze different heterogeneous models. We now formulate a simple epidemic model for a homogeneous population that includes an arbitrarily distributed exposed stage, while two different heterogeneous models will be formulated in Sections 3 and 4. Let S(t), E(t), I(t), and R(t) be the numbers of individuals in the susceptible, exposed, infectious, and removed compartments, respectively, with the total population N (t) = S(t) + E(t) + I(t) + R(t). Suppose that A > 0 represents the constant recruitment, m > 0 represents the natural mortality rate, and α ≥ 0 represents the mortality rate due to the disease. The rate of change of N (t) is Assuming mass action for the disease transmission and letting β > 0 denote the effective contact rate, the rate of change of S(t) is Let P (t) denote the fraction of exposed individuals remaining in the exposed class t units after entering the exposed class. Throughout we assume the following properties of P (t), which are biologically reasonable.
(H) P (t) is nonincreasing, piecewise continuous with possibly finitely many jumps, satisfies P (0) = 1 and lim t→∞ P (t) = 0, and the mean latent period ω = ∞ 0 P (u)du is finite. The number of exposed individuals can be expressed by the integral where E(0) is the number exposed at t = 0. Differentiating (2.3) gives Here the integral is in the Riemann-Stieltjes sense, and d t P (t − u) = dP (t−u) dt whenever the derivative exists. It can be verified that E(t) given in (2.3) is the unique solution of (2.4) with the initial condition E(0). Assuming that the recovery rate is γ, γ ≥ 0, the rate of change of I(t) is (2.5) If γ > 0, then 1/γ is the average infectious period; whereas if γ = 0, then there is no removed compartment. Substituting (2.1), (2.2), (2.4), and (2.
Therefore, the epidemic model can be written as the system with nonnegative initial conditions. Model (2.6) includes as special cases several earlier models, such as the standard SEIR ordinary differential equation (ODE) model [21,24] and the SEIR model with a discrete delay [17,39]. Epidemic models such as (2.6) with an arbitrarily distributed exposed stage have been studied in the literature; see, for example, [8,14,15]. Recently, a model of this type, but including the possibility of disease relapse, has been proposed in [40,43] to study the transmission and spread of some infectious disease such as herpes, and its global dynamics have been completely investigated in [31,40]. The model in [31,40,43] can be regarded as a generalization of our model (2.6), and thus the stability results there can be immediately applied to our model (2.6) by setting the relapse rate to zero. We outline some of their results in the following; their proofs and more detailed study can be found in [31,40]. The existence, uniqueness, and continuity of solutions of system (2.6) follow from the standard theory of Volterra integro-differential equations, see, for example, [36, p.338]. The feasible region D = (S, E, I, R) ∈ R 4 S, E, I, R ≥ 0, S + E + I + R ≤ A m is positively invariant with respect to (2.6). Let • D denote the interior of D. System (2.6) always has a disease-free equilibrium (DFE) and note that with a positive mean latent period, Define the basic reproduction number as which completely determines the stability of the DFE. We refer the reader to [3,7,40,42,43] for biological interpretation of R 0 .
(1) If R 0 ≤ 1, then the DFE is globally asymptotically stable in D.
Proposition 2.2 (Liu et al. [31]). If R 0 > 1, then the EE of system (2.8) is globally asymptotically stable in When P (t) = e − t with > 0, system (2.6) becomes an ODE model, for which the global dynamics have been established by Li and Muldowney [24] using the theory of compound matrices and by Korobeinikov and Maini [21] using the method of Lyapunov functions. When P (t) is a step function, system (2.6) becomes a delay differential equation model [39]. Propositions 2.1 and 2.2 establish the global dynamics of the model (2.6) with an arbitrarily exposed stage.
3. Multi-group model. System (2.6) assumes homogeneous mixing of individuals in the host population, and this is certainly unrealistically simple. Heterogeneity in the host population can result from different contact patterns such as those among children and adults for childhood diseases (e.g., measles and mumps), or different behavior (e.g., numbers of sexual partners for some sexually transmitted infections). In this section, we extend model (2.6) to the situation in which the population is divided into n groups according to different contact patterns or differential infectivity. Each group is further partitioned into four compartments: S i , E i , I i , and R i , denoting the number of susceptible, exposed, infectious, and removed individuals in group i, respectively. For 1 ≤ i, j ≤ n, the disease effective contact rate between compartments S i and I j is denoted by β ij , so that the new infection occurring in the i-th group is given by Hence, the n-group model associated with (2.6) can be written as the following system with the stated assumptions: Assume that A i > 0, m i > 0, γ i ≥ 0, α i ≥ 0, and each P i satisfies assumption (H) in Section 2. The contact matrix B = (β ij ) is assumed to be nonnegative and irreducible; thus any two groups i and j have a direct or indirect route of transmission. Model (3.1) has been previously seen in [44]. In the special case when P i , 1 ≤ i ≤ n, are gamma distributions, the global dynamics have been established in [44] by using a linear chain trick. However, the linear chain trick fails to apply to model (3.1) with general distributions P i . In this section we demonstrate that the graph-theoretic approach developed in [10,11,29] can be successfully applied to construct suitable Lyapunov functionals, and thus prove the global stability of the endemic equilibrium for model (3.1) with general distributions P i . This model differs from the multi-group delay epidemic model in [30] in the delay terms: arbitrarily distributed exposed stages are modeled in (3.1) while age structure is modeled in [30]. Model (3.1) can also be applied to study heterogeneous spatial distribution of the host population with implicit migration among groups; see, for example, [33].
For model (3.1) and the other heterogeneous model in Section 4, the existence, uniqueness, and continuity of solutions follow as for the basic model (2.6) from [36, p.338]. It can be easily verified that every solution of (3.1) with nonnegative initial conditions remain nonnegative. From the first equation of (3.1), it follows that , and thus lim sup t→∞ S i (t) ≤ Ai mi for all i. Adding the four equations of (3.1) together yields It can be verified that 0 < Q i < 1 for all i. With ρ denoting the spectral radius, define the basic reproduction number as the spectral radius of the n × n matrix Y = βij QiS 0 i mi+αi+γi , that is, The following result establishes that R 0 is a threshold value for the DFE. (1) If R 0 ≤ 1, then the DFE is globally asymptotically stable in D g .
Proof. Since the variables E i and R i do not appear in the equations of S i and I i of (3.1), we can first study the following reduced system consisting of variables S i and I i We prove that the equilibrium (S 0 1 , 0, S 0 2 , 0, . . . , S 0 n , 0) of system (3.4) is globally asymptotically stable by the method of Lyapunov functions. Then our stability results for the reduced system (3.4) can be extended to system (3.1) by the theory of asymptotically autonomous systems. Since is also irreducible, and Y has a positive left eigenvector (w 1 , w 2 , . . . , w n ) corresponding to the spectral radius ρ(Y ) = R 0 > 0. Motivated by [10], let Observe that L ≥ 0, and L = 0 if and only if S i = S 0 i , I i (r) = 0 for all i and 0 ≤ r ≤ t. Before we differentiate L along the solution of system (3.4), we follow the idea in [34,35] and using integration by parts obtain Now differentiating L along the solution of system (3.4), using (3.6) and the nonincreasing property of P i (see assumption (H)), and setting I = (I 1 , I 2 , . . . , I n ) T give ZHISHENG SHUAI AND P. VAN  such that (A.2) holds. Therefore, for the reduced system (3.4), the equilibrium (S 0 1 , 0, . . . , S 0 n , 0) is globally asymptotically stable when R 0 ≤ 1. Using the fact that S i (t) → S 0 i and I i (t) → 0 as t → ∞ along with the theory of asymptotically autonomous systems, E i (t) → 0 and R i (t) → 0 as t → ∞. Therefore, the DFE P 0 for system (3.1) is globally asymptotically stable in D g when R 0 ≤ 1.
To analyze the dynamical behavior of (3.1) when R 0 > 1, we consider the limiting system of (3.1) For any κ i ∈ (0, m i ), the Banach space of fading memory type C κi can be defined similarly as given in (2.9). We study system (3.8) in the phase space It can be easily verified that is positively invariant with respect to system (3.8). Let • ∆ g be the interior of ∆ g . The following result establishes the existence and uniqueness of an endemic equilibrium for system (3.8).
Using the graph-theoretic approach recently developed in [10,11,29], we prove the global stability of the EE for (3.8).
Theorem 3.3. If R 0 > 1, then the unique EE of system (3.8) is globally asymptotically stable in • ∆ g . As a consequence, all solutions of (3.1) starting in • D g approach the EE of (3.8) if R 0 > 1.

ZHISHENG SHUAI AND P. VAN DEN DRIESSCHE
Proof. As in the proof of Theorem 3.1, we first study the reduced system which is the limiting system of (3.4). In the following, we prove that the equilibrium (S * 1 , I * 1 , S * 2 , I * 2 , . . . , S * n , I * n ) of (3.12) is globally asymptotically stable by constructing a suitable Lyapunov functional. Set Here Q i and Q i (r) are defined in (3.2) and (3.5), respectively, and Q i (0) = Q i . Lyapunov functions as U i have been successfully applied to epidemic models since the work in [20,21], while Lyapunov functionals as W i have been shown to be powerful for delay epidemic models [34,35]. Functional W i is different from those in [30,34,35] because of different delay terms in model (3.1), but is motivated by those in [17,18,27,28]. Differentiating U i along the solution of (3.12) and using the equilibrium equations (3.9) and (3.10) yield For the integral in W i , integration by parts gives

It follows that
(3.14) Adding equations (3.13) and (3.14) yields where the following inequalities have been used: Considering (3.15), take the weight matrix W = (w ij ) with constant entry w ij = β ij Q i S * i I * j ≥ 0 and denote the corresponding weighted digraph as (G, W ). Let c i = T ∈Ti w(T ) ≥ 0 be as given in (B.1) in Appendix B with (G, W ). Then, by (B.2), the following identity holds n i,j=1

ZHISHENG SHUAI AND P. VAN DEN DRIESSCHE
Using (3.15) and (3.16) gives Therefore, V is a Lyapunov functional for system (3.12). This rules out the possibility that solutions of (3.12) approach the boundary S i = 0 or , and thus V | (3.12) = 0 implies that S i = S * i for all i. It can be verified that as in the proof of Theorem 3.1 the largest invariant set where V | (3.12) = 0 is the singleton {(S * 1 , I * 1 , . . . , S * n , I * n )}. Therefore, by the LaSalle-Lyapunov Theorem and a similar argument as in the proof of Theorem 3.1, it follows that (S * 1 , I * 1 , . . . , S * n , I * n ) is globally asymptotically stable for system (3.12). Thus for the limiting system (3.8) if R 0 > 1, then the EE P * is globally asymptotically stable in • ∆ g . An immediate consequence of Theorem 7.2 in [36] is that P * attracts all solutions of (3.1) Theorems 3.1 and 3.3 completely determine the dynamical behavior of the multigroup model (3.1) under the stated assumptions: if the basic reproduction number R 0 ≤ 1, then the disease dies out from all groups; if R 0 > 1, the disease persists at a positive level in all groups. It follows that the multi-group structure does not alter the qualitative behavior of the SEIR model (2.6) with an arbitrarily distributed exposed stage. 4. Multi-stage model. Multi-stage models has been proposed in the literature to describe the progression of infectious diseases with a long infectious period, for example, HIV/AIDS. Individuals infected with HIV are highly infectious in the first few weeks after infection, then remain in an asymptotic stage of low infectiousness for many years, and become gradually more infectious as their immune systems become compromised and they progress to AIDS [9,19,38]. So infectious individuals can be categorized into n different stages according to the age of infection. We include such stage structure into our model (2.6) and arrive at the following system in which I i (t) denotes the number of individuals in infection stage i with effective contact rate β i ≥ 0, γ i > 0 is the rate of transferring from stage I i to I i+1 (or to compartment R when i = n), δ i ≥ 0 is the rate of transferring from stage I i to I i−1 when i ≥ 2, and α n+1 ≥ 0 is the mortality due to disease in the removed (AIDS) compartment: I n (t) = γ n−1 I n−1 (t) − (m + α n + δ n + γ n )I n (t), (4.1) Without loss of generality, assume that there exists k ≥ 1 such that β k > 0, β j = 0, and δ j > 0 for all k < j ≤ n. In fact, if β j = 0 and δ j = 0 for all k < j ≤ n, then the variables I j , k < j ≤ n, do not appear in the first k + 2 equations and play the same role as the removed compartment R. In the case when k = n − 1, that is, β n−1 > 0, β n = 0, and δ n > 0, the compartment I n can be treated as a temporally removed compartment and δ n I n represents the relapse of diseases, so the models with disease relapse in [40,43] become special cases of (4.1). To the best of our knowledge, model (4.1) is the first multi-stage model in the literature with an arbitrarily distributed exposed stage. With assumptions as in the previous sections, it can be verified that D s = (S, E, I 1 , . . . , I n , R) ∈ R n+3 S, E, R, I i ≥ 0, 1 ≤ i ≤ n, is positively invariant with respect to system (4.1). Let • D s denote the interior of D s . For system (4.1), there always exists a disease-free equilibrium P 0 = (S 0 , 0, 0, . . . , 0, 0) where Q is defined as in (2.7) and tridiagonal The following result establishes that R 0 is a threshold value for the DFE.
Theorem 4.1. The following results hold for system (4.1).
(1) If R 0 ≤ 1, then the DFE is globally asymptotically stable in D s .
To study the dynamical behavior of (4.1) when R 0 > 1, we consider the limiting system of (4.1) (4.3) We consider system (4.3) in the phase space where κ ∈ (0, m) and C κ is given in (2.9). It can be verified that  Proof. Choose the weight matrix W = (w ij ) as given by otherwise.
Let c i = T ∈Ti w(T ) be as given in (B.1) in Appendix B with the weighted digraph (G, W ). Since β k > 0, β j = 0, δ j > 0, and γ i > 0 for some k ≥ 1, all 1 ≤ i ≤ n, and all k < j ≤ n, it follows that W is irreducible. As a consequence, c i > 0 for all i (see Appendix B). A Lyapunov functional can be constructed for system (4.3). Using similar derivations as in the proof of Theorem 3.3 and identity (B.2), differentiating V along the limiting system (4.3) gives + ln S(t − r)I j (t − r)I * a consequence, the stage structure in the long infection period does not alter the dynamics of the SEIR model (2.6).
5. Summary and discussion. We have investigated the impact of heterogeneity, which comes from age, behavior, or disease stages, on the dynamics of the SEIR epidemic model (2.6) with an arbitrarily distributed exposed stage. The resulting models are the multi-group model (3.1) and the multi-stage model (4.1), respectively. Both models are analyzed by determining the basic reproduction number R 0 (i.e., (3.3) or (4.2)) and proving that it determines a sharp threshold as for the basic model (2.6). That is, if R 0 ≤ 1, then the disease-free equilibrium is globally asymptotically stable and the disease dies out from all groups or stages; if R 0 > 1, then the endemic equilibrium is globally asymptotically stable and the disease persists in all groups or stages. Thus heterogeneity that arises from age or behavior and disease stages does not alter the qualitative behavior of the basic SEIR model (2.6), although the value of R 0 depends on the heterogeneity and the distribution of the exposed stage; see (3.3) and (4.2).
In the special case of the multi-group model (3.1) when P i (t), 1 ≤ i ≤ n, are gamma distributions, the global stability of the DFE and EE has been established in [44]. When P i (t), 1 ≤ i ≤ n, are negatively exponentially distributions, system (3.1) becomes the multi-group SEIR ODE model for which the global dynamics have been established in [10,11]. Theorems 3.1 and 3.3 extend these stability results to the case with general distributions P i (t).
Our models in this paper are of SEIR type as individuals recovering from disease are assumed to obtain permanent immunity. Epidemic models of SEIRS type have been used in the literature modeling infectious diseases with only temporary immunity, for example, see [6,25,26,32]. Our SEIR models can be easily modified to include temporary immunity and become SEIRS epidemic models. For the SEIRS model, the basic reproduction number R 0 has the same expression as for the corresponding SEIR model, and the global stability of the disease-free equilibrium can be proved in the same way as the proof of Theorems 3.1 and 4.1. However, the uniqueness and global stability of the endemic equilibrium for the SEIRS models remain open, even for the ODE models. Disease control strategies (e.g., isolation, quarantine and vaccination) can also be incorporated into our models. Our methods can sometimes be applied to analyze global dynamics of these resulting models and thus may assist public health planning to control disease.