GLOBAL STABILITY OF A MULTI-GROUP MODEL WITH VACCINATION AGE, DISTRIBUTED DELAY AND RANDOM PERTURBATION

. A multi-group epidemic model with distributed delay and vaccination age has been formulated and studied. Mathematical analysis shows that the global dynamics of the model is determined by the basic reproduction number R 0 : the disease-free equilibrium is globally asymptotically stable if R 0 ≤ 1, and the endemic equilibrium is globally asymptotically stable if R 0 > 1. Lyapunov functionals are constructed by the non-negative matrix theory and a novel grouping technique to establish the global stability. The stochastic perturbation of the model is studied and it is proved that the endemic equilibrium of the stochastic model is stochastically asymptotically stable in the large under certain conditions.


1.
Introduction. An epidemic model was proposed by Cooke [1] to describe the disease spread via a vector (such as a mosquito). It is assumed that when a susceptible vector is infected by a person, there is a time delay, τ > 0, during which the infectious agents develop in the vector, and the infected vector becomes itself infectious after the delay. It is also assumed that the vector population size is large enough such that at any time t the infectious vector population is simply proportional to the infectious human population at time t − τ . Let S(t) and I(t) denote the numbers of the human susceptible and infective individuals, respectively. The force of infection at time t is assumed to be given by βS(t)I(t − τ ). Beretta and Takeuchi [2] studied the model with distributed delay, where β is the contact rate; µ is the birth and death rate; λ is the recovery/removal rate. f (τ ) represents the proportion of the infectious vector population, and function f (τ ) is assumed to be non-negative, square integrable on R + = [0, +∞) with +∞ 0 f (τ )dτ = 1 and +∞ 0 τ f (τ )dτ < +∞. One essential assumption in classical compartmental epidemic models is that the individuals are homogeneously mixed, and each individual has the same chance to get infected. More realistic models divide the host population into groups to consider the disease transmission in heterogeneous cases. The groups can be classified according to education levels, ethnic backgrounds, gender, age, professions, communities, or geographic distributions for their diversities in disease transmission. The vital epidemic parameters vary among different groups. Based on these factors, Shu et al. [3] investigated the following general multi-group model with distributed delay: The global stability of the unique endemic equilibrium of model (2) was proved by using a graph-theoretical approach and Lyapunov functionals. Global stability results were also obtained for other multi-group epidemic models [4,5,6,7,8,9,10,11,12]. Vaccination is one of the commonly used control measures to prevent and reduce the transmission of infectious diseases. The eradication of smallpox has been considered as the most spectacular success of vaccination. Some vaccines can offer lifelong immunity with one dose, while others require boosters to maintain immunity since the acquired immunity may wane with time. It is natural to consider the vaccination and the waning immunity in modeling disease dynamics. Li et al. [13] has investigated the global dynamics of an epidemic model with vaccination for newborns and susceptibles. Blower and McLean [14] have argued that a mass vaccination campaign may increase the severity of disease, if the vaccination is applied to only 50% of the population and the vaccine efficacy is 60%. Xiao et al. [15] assumed that the vaccinated individuals can be infected at a reduced rate compared to the susceptibles. Other mathematical models on vaccination have been studied in [5,16,17].
Although waning immunity has been included in several models [5,13,15], it was assumed that the rate of the immunity loss is a constant. A better assumption on the waning immunity is that the protection immunity depends on the vaccination age of an individual (the time from the vaccination). The epidemiological models with vaccination age structure can be a suitable choice to describe the dynamics of an infectious diseases with waning immunity after vaccination.
The mathematical models with the chronological age, disease age, and vaccination age have been widely used to describe the impact of the age on the disease evolution [18,19,20,21,22]. Iannelli et al. [18] have studied an epidemic model with vaccination age by assuming the immunity decreases with the time after vaccination. Li et al. [19] have proposed an epidemic model with vaccination age and treatment to show that backward bifurcation occurs due to a piecewise treatment function. Duan et al. [20] have simplified the model [19] by assuming no treatment and have obtained the global stabilities of the disease-free equilibrium and the endemic equilibrium.
Motivated by [2,3,20], we formulate and study the following multi-group epidemic models.
Here S k , E k , I k , and R k (k = 1, 2, · · · , n) denote the numbers of susceptible, latent, infectious, and recovered individuals at time t in the k-th group, respectively. Function v k (θ, t) is the age density of vaccinated individuals at time t in the k-th group. The kernel function f j (τ ) satisfies the conditions f j (τ ) ≥ 0 and ∞ τ =0 f j (τ )dτ = a j > 0. The non-negative constant β kj is the transmission rate due to the contact of susceptible individuals in the k-th group with infectious individuals in the j-th group. The new infection occurred in the k-th group with distributed delays and the nonlinear transmission is given by n j=1 β kj S k ∞ 0 f j (τ )g j (I j (t−τ ))dτ , where g j (I j ) denotes the force of the infection. The vaccinated compartment is structured by the vaccination age θ, and it is assumed that the newly vaccinated individuals enter the vaccinated class v k (θ, t) with vaccination age zero. Function α k (θ) is the immunity wane rate, and it is a nonnegative, bounded and continuous function of θ. For two given vaccination ages θ 1 , θ 2 , (0 ≤ θ 1 ≤ θ 2 ≤ +∞), the number of the vaccinated individuals with the vaccination age θ between θ 1 and θ 2 at time t is θ2 θ1 v k (θ, t)dθ. The immunity lose rate (the number of individuals moving from the vaccinated class into the susceptible class due to the waning immunity) at time t is +∞ 0 are the natural death rates of S k , v k , E k , I k and R k , respectively. Λ k , ξ k , δ k , γ k , and k are the recruitment rate of the susceptible class, the rate of vaccination of the susceptible individuals, the rate at which exposed individuals become infectious, the recovery rate of infectious individuals, and the disease induced mortality in the k-th group, respectively. We also assume that the function g k (I k ) is sufficiently smooth and satisfies following properties [3] where g k (I k ) describes the infectivity of the individuals in I k compartment, and it is natural to assume that g k (0) = 0, g k (I k ) > 0 for I k > 0 due to the fact that the disease can not spread if there is no infection. Note that g k (I k ) I k is the per capita infectivity of the infected individuals in compartment I k . The assumption sup I k >0 g k (I k ) I k = b k says that the per capita infectivity is bounded. The limit 1086 JINHU XU AND YICANG ZHOU lim I k →0 + g k (I k ) I k = b k indicates that the per capital infectivity is the largest at the beginning of the disease outbreak.
There exist functions g j (I j ) satisfying assumptions given in (4). Four examples are listed here: Preliminary results on the dynamics of model (3) are presented in Section 2. The global stability of the equilibrium of model (3) is proved by Lyapunov functional and graph theory in Section 3. The stochastic version of the model is derived and its asymptotic behavior is studied in Section 4. A brief summary is given in the concluding section.
The comparison principle implies Integrating the second equation in (3) along the characteristic line t−θ = constant yields where Γ 0k (θ) = e − θ 0 (d V k +α k (s))ds . Substituting (6) into the first equation of (3) gives where After replacing the first equation in (3) by (7) and dropping the equation for R k (t), we have the limiting model (8). The qualitative behavior of the limiting GLOBAL STABILITY OF A MULTI-GROUP MODEL 1087 model is equivalent to that of model (3) [23].
From above analysis, we know that is the positive invariant set of model (8). The global stability of model (8) will be discussed for the solutions with the initial values in Ω. Once the solution of model (8) is determined, we can obtain v k (θ, t) from (6). The stability of the equilibrium of model (3) is the same as that of model (8). We will focus on the dynamical analysis of the reduced model (8).
The basic reproduction number is defined to as the spectrum radius of matrix F 0 , i.e., It is the expected number of infected individuals produced by a typical infected individual during its entire infectious period [26].
3. Stability analysis of equilibria. R 0 may serve as the threshold to describe disease spread. Usually, the disease will extinct in the host population if R 0 ≤ 1, and the disease will become endemic if R 0 > 1.
Proof. The irreducibility of B implies that matrix F 0 is also irreducible. F 0 has a positive left eigenvector ω = (ω 1 , · · · , ω n ) corresponding to the spectral radius we get the derivative of U 1 along the solution of (8) It is clear U 1 = 0 if and only if I = 0 or R 0 = 1 and S k = S 0 k . It can be verified that the largest compact invariant set in Ω = {(S 1 , E 1 , I 1 (.), · · · , S n , E n , I n (.)) | U 1 = 0 } is the singleton {P 0 }. By the LaSalle invariance principle for delay systems [25,27,28], we obtain that the equilibrium P 0 of system (8) is globally asymptotically stable. This completes the proof of Theorem 3.2.
Next, we can prove that the endemic equilibrium P * is globally asymptotically stable when it exists. The method is based on the graph approach and Lyapunov functionals [9,10,11].
then the endemic equilibrium P * of (8) is globally asymptotically stable when it exists.
Proof. For convenience of notations, define B is also irreducible. By Lemma 2.1 in [9], the solution space of the linear system Bv = 0 has dimension 1 with a base where c kk > 0 is the co-factor of the k-th diagonal entry of B. We construct the Lyapunov functional dsdτ .

JINHU XU AND YICANG ZHOU
Computing the derivative of U 2 along the solution of model (8), we obtain Using the equilibrium equations (10), we have The inequality in dU 2 dt is obtained by using conditions (11) and following two inequalities We first show that H 1 ≡ 0 for all I 1 , I 2 , · · · , I n > 0. It follows from the equality Let G denote the directed graph associated with matrix (β kj ), which has vertices {1, 2, · · · , n} with a directed arc (k, j) from k to j if and only if β kj = 0. E(G) denotes the set of all directed arcs of G. Using Kirchhoff's Matrix-tree Theorem (see Lemma 2.1 [9]), we know that v k = c kk can be interpreted as a sum of weights of all directed spanning subtrees T of G that are rooted at vertex k. Consequently, each term in v k β kj is the weight w(Q) of a unicyclic subgraph Q of G, obtained from such a tree T by adding a directed arc (k, j) from vertex k to vertex j. Note that the arc (k, j) is part of the unique cycle CQ of Q, and that the same unicyclic graph Q can be formed when each arc of CQ is added to a corresponding rooted tree T . It is not difficult to see that the double sum in H 2 can be interpreted as a sum over all the arcs in the cycles of all the unicyclic subgraphs Q containing of G. Therefore, H 2 can be rewritten as Since E(CQ) is the set of arcs of a cycle CQ, we have This implies that H n,Q = 0 for each Q, and H 2 ≡ 0 for all I 1 , I 2 , · · · , I n > 0. From assumption (11) we have U 2 ≤ 0. It can be verified that the largest compact invariant subset of set (S 1 , E 1 , I 1 (.), · · · , S n , E n , I n (.)) dU 2 dt = 0 is the singleton {P * }. By the LaSalle invariance principle and an argument similar to that in the proof of Theorem 3.2 we know that the equilibrium P * of model (8) is globally asymptotically stable.

Remark 1. (i)
It is easy to see that the functions listed in (5) satisfy condition (11). But the following three functions may be excluded by condition (11) I j e −αj Ij , I j 1 + α j I 2 j , I j 1 + α j I j + β j I 2 j .
(ii) Condition (11) holds if g j (0) = 0, g j (I j ) > 0, and g j (I j ) ≤ 0 (I j > 0). 4. Stochastic model. The nature of epidemic growth and spread is inherently random due to the unpredictability of person-to-person contacts [29], and the population is subject to a continuous spectrum of disturbances [30]. The deterministic approach has some limitations in modelling the transmission of an infectious disease due to environmental noises. Stochastic differential equation (SDE) models have been applied to different infectious diseases in many circumstances [31,32,33,34,35,36]. Motivated by [3], we take the randomly fluctuating environment into consideration by stochastic perturbations of white noise type. We study the following stochastic model, where B 1k (t), B 2k (t) and B 3k (t) (k = 1, 2, · · · , n) are independent standard Brownian motions defined on a complete probability space (Ω, F, {F t } t≥0 , P) with a filtration {F t } t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F 0 contains all P-Null set). And σ 2 ik > 0 (i = 1, 2, 3) represent the intensities of B ik (i = 1, 2, 3), respectively. If there are no noises, i.e. σ ik = 0 (i = 1, 2, 3), then model (12) is The basic reproduction number of model (13) is The similar arguments in [3] can lead to the following result.
Theorem 4.1. If R 1 ≤ 1, then the disease-free equilibrium of model (13) is globally asymptotically stable. If R 1 > 1, then model (13) has an endemic equilibriumP * , which is globally asymptotically stable. P * is also an equilibrium of the stochastic model (12). We will focus on the stability of the equilibriumP * by using Lyapunov functionals. We firstly give some definitions and auxiliary statements.
Consider the n-dimensional stochastic functional differential equation where X t = X (t+s), s ≤ 0, BC((−∞, 0], R n ) is the space of bounded and continuous functions from (−∞, 0] to R n with the norm φ = sup s≤0 |φ(s)|. Suppose that the existence and uniqueness theorem holds, and (14) has a zero solution. Let C 2,1 (R n ×R + ; R + ) be the family of all nonnegative functions V (X , t) defined on R n × R + such that they are continuously differentiable twice in X and once in t. For a function V ∈ C 2,1 (R n × R + ; R + ), define the operator L by where T means the transposition. (14) is said to be stochastically stable or stable in probability if for every pair of ε ∈ (0, 1) and r > 0, there exists δ > 0 such that

Definition 4.2. [37] (1) The trivial solution of
(2) The trivial solution is said to be stochastically asymptotically stable if it is stochastically stable, and for every ε ∈ (0, 1), there exists δ > 0 such that If R 1 > 1, then the stochastic system (12) can be centered at its endemic equi-libriumP * (S * 1 , E * 1 , I * 1 , · · · , S * n , E * n , I * n ). The change of variables It is easy to see that the stability ofP * of (12) is equivalent to the stability of zero solution of system (18).
Proof. From the stability theory of stochastic functional differential equations, it is sufficient to find a Lyapunov functional V (X ) such that LV (X ) ≤ 0 for sufficiently small δ > 0 and the identity holds if and only if X = 0 (see [37]). The endemic equilibriumP * satisfies following equations We define where m k (k = 1, 2, · · · , n) are constants to be determined. Itô's formula leads to Define Similarly, from Itô's formula, we obtain and From (21), (22) and (23) we have Finally, we define which gives Let where c kk has the similar definition as that given in Section 3. Then we have c kk = m k E * k , and which gives Let V = V 1 + V 2 + V 3 + V 4 , from above analysis we obtain where In (28) we choose n k such that Then Let Substituting (30) and (31) into (29) yields where From (19), we have A k > 0, B k > 0 and D k > 0. Consequently Assume P{|z j (s)| ≤ ρ} = 1 (ρ > 0, j = 1, 2, · · · , n). Then Therefore, Consequently, for sufficiently small ρ > 0, we have LV < 0 except the zero point.  (12). The main idea in simulations is the Milstein method [38] by discretizing the differential equations. We do simulations for k = 1, 2 and f j (τ ) = e −τ , τ ≥ 0. The initial condition is I j (θ) = v j e θ , θ ≤ 0, where v j > 0. For k = 1, 2 and f j (τ ) = e −τ , τ ≥ 0, model (12) is The discretization of model (33) is 3k (k = 1, 2.) are independent random variables with normal distribution N (0, 1), which can be generated numerically by pseudo-random number generators. For simplicity, we always assume that d S k = d E k = d I k = d k and parameters take following values: The endemic equilibriumP * of model (33) In the absence of noise, i.e. σ ik = 0, i = 1, 2, 3; k = 1, 2, the global stability of the endemic equilibrium corresponding deterministic model (33) is shown in Figure 1.
In Figure 2 (34), are given in Figure 3. The endemic equilibrium of model (33) is also asymptotically stable (see Figure 3). The comparison of Figures 2 and 3 shows that the solutions of stochastic model (33) fluctuate at the beginning, and converge to the equilibrium position finally. The fluctuations of the stochastic model (33) may enhance with the increasing noises.

5.
Conclusions. We have studied a multi-group SVEIR epidemic model with delays and vaccination age. The global stability of model (8) is established by Lyapunov functionals. The disease-free equilibrium is globally asymptotically stable if R 0 ≤ 1, and the unique endemic equilibrium is globally asymptotically stable if R 0 > 1. We can define R 0 = ρ β kj S 0 k aj d I k +γ k + k 1≤k,j≤n = lim δ k →∞ R 0 to investigate the influence of the latent period on R 0 . It is obvious that R 0 > R 0 and ∂R0 ∂δ k > 0, which implies that latent period has a positive role in disease control: a long latent period may lead to the extinction of the disease. Similarly, the fact that implies that the vaccination is very helpful to eradicate the disease. Though the immunity of a vaccine may not be permanent, a long immunity period of vaccines is still expected for diseases prevention.  We have investigated a multi-group stochastic model (12) with white noise perturbations. We obtained sufficient conditions for stochastic stability of model (12). The results reveal that the stochastic stability of the endemic equilibrium depends on the magnitude of the noise. Numerical simulations show that the stochastic  Other factors, such as population migration, can be integrated into the model to make it more realistic. The efficacy of some vaccines may not be 100% though it is assumed that the vaccinated individuals can not be infected. It is more reasonable to assume that the vaccinated individuals can be infected at a reduced rate [15]. Furthermore, we just suppose that the stochastic perturbations are proportional to S k −S * k , E k −E * k and I k −I * k . It is also interesting to study other types of stochastic perturbations [33] and consider a stochastic system with nonlinear incidence.