Analysis and Numerical Simulations of a Stochastic SEIQR Epidemic System with Quarantine-Adjusted Incidence and Imperfect Vaccination.

This paper considers a high-dimensional stochastic SEIQR (susceptible-exposed-infected-quarantined-recovered) epidemic model with quarantine-adjusted incidence and the imperfect vaccination. The main aim of this study is to investigate stochastic effects on the SEIQR epidemic model and obtain its thresholds. We first obtain the sufficient condition for extinction of the disease of the stochastic system. Then, by using the theory of Hasminskii and the Lyapunov analysis methods, we show there is a unique stationary distribution of the stochastic system and it has an ergodic property, which means the infectious disease is prevalent. This implies that the stochastic disturbance is conducive to epidemic diseases control. At last, computer numerical simulations are carried out to illustrate our theoretical results.

1. Introduction

Mathematical models for differential equations have been widely applied in various fields [1-7]. Specifically, they have had a realistic significance to analyze the dynamical behaviors in the field of mathematical biology [8-17], which obtained some novel results.

In fact, the main meaning of the research of infectious disease dynamics is to make people more comprehensively and deeply understand the epidemic regularity of infectious disease; then more effective control strategies are adopted to provide better theoretical support for the prevention and control of epidemics. To this end, many mathematical biology workers considered more realistic factors in the course of the study, such as population size change, migration, cross infection, and other practical factors. In the course of epidemics and outbreaks of infectious diseases, people always take various measures to control the epidemic in order to minimize the harm of epidemic diseases. Quarantine is one of the important means to prevent and control epidemic diseases; it has been used to control contagious diseases with some success. Specifically, during the severe acute respiratory syndrome (SARS) outbreak in 2002, remarkable results were also achieved. Among them, mathematical models have been used to investigate their impact on the dynamics of infectious diseases under quarantine [18-22], which attracts deep research interest of many mathematicians and biologists. Recently, Hethcote et al. [21] considered an SIQR (susceptible-infected-quarantined-recovered) model with quarantine-adjusted incidence. The system can be expressed as follows:where the total population size is given by , Λ is the inflow rate corresponding to birth and immigration, and μ is the outflow rate corresponding to natural death and emigration. Since the quarantine process, using the standard incidence , the contact rate with the quarantined fraction does not occur. Hence the standard incidence is replaced by (it is called quarantine-adjusted incidence); here β is the transmission coefficient between susceptible individuals and infected individuals. λ is the quarantine rate of infected individuals, φ is the recovery rate of quarantined individuals, and ε1 and ε2 stand for the rate of disease-related death of infected and quarantined individuals, respectively. γ is the recovery rate of infected individuals. Furthermore, all the parameters are positive and the region is a positively invariant set of system (1). In the region , they established the basic reproduction number R0, which determines disease extinction or permanence, whereMeanwhile, they analyzed the global dynamics of system (1) and derived the equilibria (including the disease-free equilibrium and the endemic equilibrium) and their global stability. In addition, the parameter conditions for the existence of a Hopf bifurcation are obtained.

In the real world, with the development of modern medicine, vaccination has become an important strategy for disease prevention and control in addition to quarantine, and numerous scholars have investigated the effect of vaccination on disease [23-30]. For another, many infectious diseases incubate inside the hosts for a period of time before becoming infectious, so it is very meaningful to consider the effect of the incubation period. Motivated by the aforementioned work, this paper considers an SEIQR (susceptible-exposed-infected-quarantined-recovered) epidemic model with imperfect vaccination, which is described by the following system:where the total population size is given by N(t) = S(t) + E(t) + I(t) + Q(t) + R(t), δ  (0 ≤ δ < 1) is the vaccine coverage rate, p  (0 ≤ p ≤ 1) is the vaccine efficacy, and α is the rate at which the exposed individuals become infected individuals. Other parameters are the same as in system (1). Now we assume that all the parameters are positive constants here except that δ, p are nonnegative constants. Clearly, the region D = {(S, E, I, Q, R)∣S ≥ 0, E ≥ 0, I ≥ 0, Q ≥ 0, R ≥ 0, S + E + I + Q + R ≤ Λ/μ} is a positively invariant set of system (3). For system (3), the basic reproduction number isand it has the following properties:

When R1 ≤ 1 holds, system (3) has a unique disease-free equilibrium E0 = (S0, 0,0, 0, R0) = (Λ/(δp + μ), 0,0, 0, δpΛ/μ(δp + μ)) which is globally asymptotically stable in the region D. That means the epidemic diseases will die out and the total individuals will become the susceptible and recovered individuals.

When R1 > 1 holds, system (3) has a unique globally asymptotically stable positive equilibrium E = (S, E, I, Q, R) in the region D, which means the epidemic diseases will persist.

In the natural world, deterministic model is not enough to describe the species activities. Sometimes, the species activities may be disturbed by uncertain environmental noises. Consequently, some parameters should be stochastic [31-40]. There is no denying that this phenomenon is ubiquitous in the ecosystem. Therefore numerous scholars have introduced the effect of stochastic perturbation on diseases [41-50]. To the best of our knowledge, the research on global dynamics of the stochastic SEIQR epidemic model with imperfect vaccination is not too much yet. In this paper, to make system (3) more reasonable and realistic, we assume the environmental noise is directly proportional to S(t), E(t), I(t), Q(t), and R(t). Then corresponding to system (3), a stochastic version can be reached bywhere B(t)  (i = 1,2, 3,4, 5) is the mutually independent standard Wiener process with B(0) = 0 a.s. σ(t)  (i = 1,2, 3,4, 5) is a continuous and bounded function for any t ≥ 0 and σ2(t)  (i = 1,2, 3,4, 5) are the intensities of Wiener processes.

In this paper, we are mainly concerned with two interesting problems as follows:

(P1) Under what parameter conditions, will the disease die out?

(P2) Under what conditions, will system (5) have a unique ergodic stationary distribution?

Throughout this paper, let (Ω, ℱ, {ℱ}, ℙ) be a complete probability space with a filtration {ℱ} satisfying the usual conditions (i.e., it is increasing and right continuous while ℱ0 contains all ℙ-null sets). Further B(t)  (i = 1,2, 3,4, 5) is defined on the complete probability space.

For simplicity and convenience, we introduce the following notations:

ℝ+ = [0, +∞), ℝ+5 = {x = (x1, x2, x3, x4, x5) ∈ ℝ5 : x > 0  (i = 1,2, 3,4, 5)}.

For an integrable function x(t)∈[0, +∞), 〈x(t)〉 = (1/t)∫0x(r)dr.

a∧b = min⁡{a, b}, a∨b = max⁡{a, b}.

2. Global Positive Solution

To investigate the dynamical behaviors of a population system, we first concern the global existence and positivity of the solutions of system (5).

Lemma 1 .

For any given initial value (S(0), E(0), I(0), Q(0), R(0)) ∈ ℝ+5, system (5) has a unique positive local solution (S(t), E(t), I(t), Q(t), R(t)) for t ∈ [−ω, τ), where τ is the explosion time [51].

Theorem 2 .

For any given initial value (S(0), E(0), I(0), Q(0), R(0)) ∈ ℝ+5, system (5) has a unique positive solution (S(t), E(t), I(t), Q(t), R(t)) ∈ ℝ+5 on t ≥ 0 a.s.

Proof

The following proof is divided into two parts.

Part I. By Lemma 1, it is easy to see that system (5) has a unique positive local solution (S(t), E(t), I(t), Q(t), R(t)) for any given initial value (S(0), E(0), I(0), Q(0), R(0)) ∈ ℝ+5.

Part II. Now we prove that the positive solution is global, that is, τ = ∞ a.s. Let k0 ≥ 0 be sufficiently large such that S(0), E(0), I(0), Q(0), and R(0) all lie in [1/k0, k0]. For each integer k ≥ k0, let us define the stopping time as follows:where we define inf⁡∅ = ∞ (∅ represents the empty set). Evidently, τ is strictly increasing when k → ∞. Let τ = limτ; thus τ ≤ τ a.s. So we just need to show that τ = ∞ a.s. If τ = ∞ is untrue, then there exist two constants T > 0 and ς ∈ (0,1) such that ℙ{τ ≤ T} > ς. Thus there exists k1 ≥ k0  (k1 ∈ N+) such thatDefine a C2-function byApplying Itô's formula and system (5), we havewhereand here M0 is a positive constant. HenceIntegrating both sides of (11) from 0 to τ∧T and then taking the expectation, we haveSet Ω = {τ ≤ T}, k ≥ k1 and by (7) we can get that P(Ω) ≥ ς. Notice that, for every ω ∈ Ω, there exists S(τ, ω), E(τ, ω), I(τ, ω), Q(τ, ω), or R(τ, ω) which equals either 1/k or k. ThusBy virtue of (12) and (13), one has and here 1 is the indicator function of Ω(ω). Let k → ∞, which implies is a contradiction. Obviously, we get that τ = ∞. This completes the proof of Theorem 2.

3. Extinction

In this section, we mainly explore the parameter conditions for extinction of the disease in system (5). Before proving the main results, we first give a useful lemma as follows.

Lemma 3 .

For any given initial value (S(0), E(0), I(0), Q(0), R(0)) ∈ ℝ+5, the solution (S(t), E(t), I(t), Q(t), R(t)) of the system (5) has the following properties:Furthermore, when μ > (1/2)(σ12∨σ22∨σ32∨σ42∨σ52) holds, then

The proof of Lemma 3 is similar to [25, 41]; thus we omit it here.

Theorem 4 .

Let μ > (1/2)(σ12∨σ22∨σ32∨σ42∨σ52). For any given initial value (S(0), E(0), I(0), Q(0), R(0)) ∈ ℝ+5, ifholds, thenMoreover,

Define a differentiable function V0 byFrom Itô's formula and system (5), we haveIntegrating from 0 to t and dividing by t on both sides of (22), we have Making use of Lemma 3, we havewhich shows thatFrom the fourth equation of system (5), it is easy to get thatMoreover, integrating from 0 to t and dividing by t on both sides of the first equation of system (5) yieldand considering (25), (26), and Lemma 3, it then follows thatSimilarly, we also getThe proof of Theorem 4 is complete.

4. Stationary Distribution and Ergodicity

Ergodicity is a significant property in a population system. Recently, it attracts deep research interest of numerous scholars [52, 53]. In this section, based on the theory of Hasminskii et al. [54] and the Lyapunov analysis methods, we study the conditions of the existence of an ergodic stationary distribution.

Assume X(t) as a time-homogeneous Markov process in 𝔼 ⊂ ℝ, which is described by the stochastic differential equationand here 𝔼 stands for a n-dimensional Euclidean space. The diffusion matrix takes the following form:

Assumption 5 .

Assume that there is a bounded domain U ⊂ 𝔼 with regular boundary Γ such that , satisfying the following properties:

In the domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix is bounded away from zero.

If x ∈ 𝔼∖U, the mean time τ at which a path issuing from x reaches the set U is finite, and sup𝔼τ < ∞ for every compact subset Θ ⊂ 𝔼.

Lemma 6 (see [54]).

When Assumption 5 holds, then the Markov process X(t) has a stationary distribution π(·). Furthermore, when f(·) is a function integrable with respect to the measure π, then for all x ∈ 𝔼.

Remark 7 .

To demonstrate Assumption 5(i) [55], it suffices to demonstrate that F is uniformly elliptical in any bounded domain H; herenamely, there exists a positive number Z such thatTo verify Assumption 5(ii) [56], it suffices to demonstrate that there exist some neighborhood U and a nonnegative C2-function V such that ∀x ∈ 𝔼∖U, LV(x) < 0.

Using Lemma 6, we can get the following main results.

Theorem 8 .

For any given initial value (S(0), E(0), I(0), Q(0), R(0)) ∈ ℝ+5. Ifholds, then system (5) has a unique stationary distribution π(·) and it has ergodic property.

Define a C2-function byand here ϱ and a  (i = 1,2, 3) are positive constants satisfying the following conditions: and we take Υ > 0 large enough such that hereObviously, and here U = (1/w, w) × (1/w, w) × (1/w, w) × (1/w, w) × (1/w, w). Since is a continuous function, then there exists a unique point in ℝ+5 which is the minimum point of . Therefore let us construct a positive-definite C2-function V: ℝ+5 → ℝ+ by From Itô's formula, we haveSimilarly,and hereWe also haveTherefore,where M = Λ + Γ + δp + α + φ + ε2 + 4μ + σ12/2 + σ22/2 + σ42/2 + σ52/2.

Next let us consider the following compact subset D: and here ϵ is a sufficiently small constant satisfying the following conditions:ThenwithNow we analyze the negativity of ℒV for any (S, E, I, Q, R) ∈ ℝ+5∖D.

Case I. If (S, I, Q, R) ∈ D1, (46) and (48) derive that

Case II. If (S, I, Q, R) ∈ D2, (46) and (49) imply that

Case III. If (S, I, Q, R) ∈ D3, it follows from (46) and (50) that

Case IV. If (S, I, Q, R) ∈ D4, (46) and (51) yield that

Case V. If (S, I, Q, R) ∈ D5, it follows from (46) and (52) that

Case VI. If (S, I, Q, R) ∈ D6, (46) and (53) lead to

Case VII. If (S, I, Q, R) ∈ D7, (46) and (54) derive that

Case VIII. If (S, I, Q, R) ∈ D8, it follows from (46) and (55) that

Case IX. If (S, I, Q, R) ∈ D9, (46) and (56) derive that

Case X. If (S, I, Q, R) ∈ D10, it follows from (46) and (54) thatClearly, from the discussion of the above ten cases, one sees that, for a sufficiently small ϵ, which shows that Assumption 5(ii) is satisfied. In addition, the diffusion matrix of system (5) takes the following form:There exists a positive numbersuch thatwhich shows that Assumption 5(i) is satisfied. Consequently, system (5) has a unique stationary distribution π(·) and it has ergodic property. The proof of Theorem 8 is complete.

Remark 9 .

From Theorem 8, we see that if R > 1 holds, then system (5) has a unique ergodic stationary distribution. It is worthwhile noting that if σ = 0  (i = 1,2, 3,4, 5), the expression of R coincides with the basic reproduction number R1 of system (3). This shows that we generalize the results of system (3). For another, this theorem also shows that the disease can resist a small environmental noise to maintain the original persistence.

5. Conclusions and Simulations

This paper studies the stochastic SEIQR epidemic model with quarantine-adjusted incidence and imperfect vaccination and obtains two thresholds which govern the extinction and the spread of the epidemic disease. Firstly, the existence of a unique positive solution of system (5) with any positive initial value is proved. Then, from Theorems 4 and 8, the sufficient conditions for extinction of the disease and existence of ergodic stationary distribution of the stochastic system are derived by using the theory of Hasminskii and the Lyapunov analysis methods, which means the infectious disease is prevalent. This implies that the stochastic disturbance is conducive to epidemic diseases control. Now we summarize the main conclusions as follows:

When μ > (1/2)(σ12∨σ22∨σ32∨σ42∨σ52) and hold, then the infected individuals go to extinction almost surely.

When R = (μ(1 − δp)βα)/((δp + μ + σ12/2)(α + μ + σ22/2)(λ + ε1 + γ + μ + σ32/2)) > 1 holds, then system (5) has a unique stationary distribution π(·) and it has ergodic property.

To illustrate the results of the above theorems, we next carry out some numerical simulations by the Matlab software. Let us consider the following discretization equations of system (5):and here ζ  (i = 1,2, 3,4, 5; k = 1,2,…, n) stands for N(0,1) distributed independent random variables and time increment Δt > 0.

In Figure 1, take S(0) = E(0) = I(0) = Q(0) = R(0) = 0.3, Λ = 2, μ = 0.55, β = 0.25, γ = 0.2, α = 2.5, δ = 0.5, p = 0.5, ε1 = 0.1, ε2 = 0.1, λ = 0.18, φ = 0.15, σ1 = 0.15, σ2 = 1, σ3 = 1, σ4 = 0.5, σ5 = 0.25, and Δt = 0.01. Thensatisfy the conditions in Theorem 4; we can obtain that the exposed, infected, and quarantined individuals go to extinction almost surely. Moreover,Obviously, Figure 1 supports our results of Theorem 4.

Figure 1

Time sequence diagram of system (5) for extinction of the disease.

In Figure 2, take S(0) = E(0) = I(0) = Q(0) = R(0) = 0.5, Λ = 0.3, μ = 0.1, β = 1.5, γ = 0.2, α = 0.3, δ = 0.1, p = 0.2, ε1 = 0.05, ε2 = 0.05, λ = 0.18, φ = 0.15, σ1 = σ2 = σ3 = σ4 = σ5 = 0.03, and Δt = 0.01. Then satisfies the condition in Theorem 8; we can obtain that system (5) has a unique stationary distribution π(·) and it has ergodic property. Figure 2 shows that the solution of system (5) swings up and down in a small neighborhood. According to the density functions in Figures 2(b)–2(f), we can see that there exists a stationary distribution. As expected, Figure 2 confirms our results of Theorem 8.

Figure 2

(a) represents the solutions of system (5); (b)–(f) stand for the density functions of S(t), E(t), I(t), Q(t), and R(t), respectively.

In Figure 3, take S(0) = E(0) = I(0) = Q(0) = R(0) = 0.1, Λ = 0.2, μ = 0.55, β = 2.65, γ = 0.2, α = 2.5, δ = 0.1, p = 0.1, ε1 = 0.1, ε2 = 0.1, λ = 0.18, φ = 0.15, and Δt = 0.01.

Figure 3

Time sequence diagram of system (5) for persistence and extinction of the disease.

In Figure 3(a), take σ = 0  (i = 1,2, 3,4, 5); then which means the disease will persist. As expected, Figure 3(a) shows the disease persists in real life.

Synchronously, in Figure 3(b), take σ1 = 0.15, σ2 = 1, σ3 = 1, σ4 = 0.5, σ5 = 0.25. Obviously, Figure 3(b) shows the exposed, infected, and quarantined individuals go to extinction and we can get that the permanent disease of system (3) can die out under stochastic effects. This implies that the stochastic disturbance is conducive to epidemic diseases control.