Numerical Analysis of the Susceptible Exposed Infected Quarantined and Vaccinated (SEIQV) Reaction-Diffusion Epidemic Model

Abstract

In this paper, two structure-preserving nonstandard finite difference (NSFD) operator splitting schemes are designed for the solution of reaction diffusion epidemic models. The proposed schemes preserve all the essential properties possessed by the continuous systems. These schemes are applied on a diffusive SEIQV epidemic model with a saturated incidence rate to validate the results. Furthermore, the stability of the continuous system is proved, and the bifurcation value is evaluated. A comparison is also made with the existing operator splitting numerical scheme. Simulations are also performed for numerical experiments.

1. Introduction

Mathematical modeling has a prominent role in describing physical phenomena in various disciplines of mathematics, physical sciences, social sciences, engineering, life sciences, and many more [1–6]. The transmission of infectious diseases and the control of their spread can be studied effectively by constructing mathematical models for various strategies like vaccination and quarantine. The word quarantine denotes forced isolation or being cut off from interactions with others. Quarantine is an effective intervention process for restraining the spread of infection by isolating individuals who are affected by the disease. Such isolation has been adopted to decrease the communication of infectious diseases like dengue, measles, smallpox, cholera, leprosy, tuberculosis, and many more.

Epidemic models, that is, mathematical models of infectious diseases, are a simplified way to illustrate the transmission dynamics of the complicated nonlinear processes and complex behavior of an infectious disease in individuals within a population. These are deterministic models that are used to allocate the population to different subclasses or compartments, describing a particular stage of the epidemic. The incidence rate, which is proportional to the number of susceptible and infected persons, is an important parameter of compartment-based epidemic models. The mathematical models of infectious diseases are often based on bilinear incidence rate βSI, but a more concise approach to use the saturated incidence rate rather than the bilinear incidence rate. In the saturated incidence rate , if number of infected individuals I becomes very large, approaches the saturation level. The infection force is measured by βI, which describes the penetration of the disease into a fully susceptible population. is used to measure the inhibition effect of behavioral change of susceptible persons. Liu and Yang [7] proposed the SEIQV epidemic model, which uses the saturated incidence rate. The model is expressed as:

The variables and parameters of the model are defined as:

S(t) = Susceptible persons at time t,

E(t) = Exposed persons at time t,

I(t) = Infected persons at time t,

Q(t) = Quarantined persons at time t,

V(t) = Vaccinated persons at time t,

b = Rate of recruitment,

β = Rate of transmission,

μ = Rate of natural death,

ϵ = Rate of death due to disease in infected compartment,

α = Parameter that measures psychological or inhibitory effects,

γ = Rate at which infected individuals are being vaccinated infected persons,

σ = Rate at which exposed persons become infected,

ω = Rate at which infected individuals are being vaccinated susceptible persons,

ϕ = Rate at which infected individuals are being vaccinated quarantined persons,

q₁, q₂, q₃ = Effective quarantine probabilities.

The above model (1)–(5) assumes a homogeneous population, where the population mixes in such a way that there is no difference between person in one place and person in another place. However, in actual scenarios, the disease may spread faster in one place than in another because of different circumstances like different weather conditions, etc. Hence, it is essential for the variables to depend on space also. Therefore, we extend system (1)–(5) to make it a reaction-diffusion system by adding a diffusion term.

with the initial conditions:

The boundary conditions are no flux,

Epidemic models always demonstrate two equilibrium points: the disease-free equilibrium (DFE) point and the endemic equilibrium (EE) point. The DFE point exists if R₀ < 1, where R₀ is the reproductive number, which basically measures the occurrence of disease. The EE point exists if R₀ > 1. This implies that the SEIQV reaction-diffusion system (6–10) always converges to the DFE point or EE point if R₀ < 1 or R₀ > 1, respectively. Analytical solution of the SEIQV epidemic system is not possible, so we have to use numerical techniques to find its solution. Note that the numerical technique must show the same behavior as is possessed by the continuous SEIQV reaction-diffusion epidemic system.

In this work, we propose two operator-splitting NSFD methods, one explicit and one implicit. These methods are used to solve the SEIQV epidemic model with diffusion. As S, E, I, Q, and V are population sizes and evaluated in absolute scale, we propose NSFD methods because they give a positive solution. Also, the convergence of the proposed NSFD operator splitting methods toward the equilibrium points is the same as the convergence of continuous an SEIQV reaction-diffusion epidemic system. The proposed splitting methods are designed with the aid of rules given by Mickens [8]. In the recent era, positivity preserving FD methods have gained importance, as many physical dynamical systems possess the positivity property [9–11]. The NSFD method presented by Mickens [8, 12, 13] has becomes an effective and important structure-preserving FD method for solving differential equations. In epidemic models, population dynamics and population size cannot be negative, so the numerical technique must be a positivity-preserving technique. Various authors have used different positivity-preserving numerical techniques for the approximate solution of epidemic models: see, for example [14–22].

In this work, we also show the numerical stability of the SEIQV epidemic model with diffusion and evaluate the bifurcation value of the vaccination parameter ω with the aid of the Routh-Hurwitz method.

2. Equilibrium Points

The model (6–10) has two equilibrium points [7], the DFE point and EE point. The DFE point is:

and the EE point is:

where,

Reproductive number R₀ is given as:

R₀ is the reproductive value. Now, if R₀ < 1, the model acquires a DFE point, and if R₀ > 1, the model acquires an EE point.

3. Numerical Stability of the SEIQV Model at Equilibrium Point

We evaluated the small perturbation S₁(x, t), E₁(x, t), I₁(x, t), Q₁(x, t), and V₁(x, t) so that (6)–(10) is linearized at the EE point (S_*, E_*, I_*, Q_*, V_*), as discussed in Chakrabrty et al. [23].

Suppose a Fourier series solution is demonstrated for Equations (23)–(27) of the form:

Here, k = nπ/2, (n = 1, 2, 3,…) exhibits the value of the wave number for the node n. Substituting Equations (28)–(32) in the system (23)–(27), the system is converted into:

The variational matrix for the system (33)–(37) is:

where,

The characteristics equation for matrix is:

The expressions for ξ₁, ξ₂, ξ₃, ξ₄, and ξ₅ with diffusion and without diffusion are mentioned in Islam and Haider [24].

The Routh-Hurwitz stability criterion gives:

and

The Table 2 reflects the numerical stability of the equilibrium point against the various cases, as discussed in Table 1.

4. Bifurcation Value of Vaccination Parameter ω Independent of Diffusion

Considering the vaccination parameter ω, to find its bifurcation value, a₁₁, a₁₂, are used instead of S_*, E_*, I_*, Q_*, and V_*.

The Routh-Hurwitz criterion for stability gives:

where the values of ξ₁, ξ₂, ξ₃, ξ₄, and ξ₅ are obtained from the expression of the characteristic equation (without diffusion) given in paper [24].

f₅(ω) = 0 gives the value of bifurcation for ω. This value transfers the stability of a continuous model from stable to unstable. f₅(ω) = 0 provides the bifurcation value ω = 0.228360507. The EE point is stable for ω less than ω = 0.228360507.

5. Bifurcation Value of Vaccination Parameter ω With Diffusion

For the bifurcation value of vaccination parameter ω, the values of S_*, E_*, I_*, Q_*, and V_* are replaced into a₁₁, a₁₂, to

The Routh-Hurwitz criterion for stability gives:

f₅(ω) = 0 gives ω = 0.24144152. The EE point is stable therefore for any value less than ω = 0.24144152.

It can be seen that the value of bifurcation of ω for the system with diffusion is greater than value of bifurcation of ω for the system without diffusion.

6. Numerical Methods

In this section, we apply two proposed and classical splitting methods to the SEIQV reaction-diffusion epidemic model with diffusion. Operator-splitting techniques very effectively handle the nonlinearity and complexity of reaction-diffusion equations. Therefore, these techniques are used frequently by several researchers for the solution of nonlinear differential equations [23, 25–33]. The SEIQV epidemic model with diffusion is split in two ways. The nonlinear reaction equations are split in the first step as,

and the diffusion equations are split in the second step as:

Now, we apply forward and backward Euler methods with operator splitting on the system (6)–(7).

where and at mτ, m = 0, 1, … and 0 + ih, i = 0, 1, … reflects difference approximations of S, E, I, Q, and V. The values of , , , , and are the values at the half time step. Both forward and backward Euler methods have same process at first step, but, at the second half step of time, they have different procedures. Since the forward Euler operator splitting method is explicit, we use:

For the backward Euler method, we use:

For the proposed NSFD operator splitting methods, we implement the rules constructed by Mickens [8]. The technique for both explicit and implicit schemes is similar at the first half time step:

A positive solution desires that if:

The techniques for the implicit and explicit NSFD schemes are not similar for the second half of the time step. The procedure for the explicit NSFD scheme is as follows:

We use an implicit procedure for the second NSFD method at the second half of the time step:

where,

6.1. Stability and Accuracy of Splitting Schemes

In finite difference operator splitting techniques, the step involving the reaction term is unconditionally stable because it is solved exactly [25, 26]. On the other hand, the step involving the diffusion term has different stability in different techniques. The explicit procedure has conditional stability in the region:

while the implicit procedure has unconditionally stability [25, 26]. The accuracy of both schemes is O(τ) and O(h²) for all the methods under study.

6.2. Positivity of Proposed Schemes

Equations (64)–(65) in the reaction step of both proposed methods preserve the property of positivity depicted by the continuous SEIQV model, as there is no negative term involved in (64)–(65).

As far as the diffusion step is concerned, the proposed explicit scheme (70)–(74) demonstrates the positivity of the solution if:

so,

which is the condition of stability for the explicit operator-splitting NSFD scheme (70)–(74). This verifies that the explicit NSFD scheme retains the positive solution in the region of stability. M matrix theory has been used for the verification of the positivity of the implicit NSFD method (75)–(79). For more details [34] is referred.

6.2.1. Theorem [21, 22]

For any positive τ and h, the system described by (75)–(79) is also positive, i.e., S^m > 0, E^m > 0, Q^m > 0 and V^m > 0, ∀m = 0, 1, 2…

Proof

Equations (75)–(79) can be written as:

In Equations (81)–(85), the letters A, B, C, D, and G represent the square matrices. Where,

The off-diagonal entries of A are a₁ = −2λ₁, a₂ = −λ₁, and the diagonal entries are a₃ = 1 + 2λ₁. The entries of B in the off-diagonal are b₁ = −2λ₁, b₂ = −λ₁, and the diagonal entries are b₃ = 1 + 2λ₂. The entries of C in the off-diagonal are c₁ = −2λ₃, c₂ = −λ₃, and the diagonal entries are c₃ = 1 + 2λ₃. The off-diagonal entries of D are d₁ = −2λ₄, d₂ = −λ₄, and the diagonal entries are d₃ = 1 + 2λ₄. The off-diagonal entries of G are g₁ = −2λ₅, g₂ = −λ₅, and the diagonal entries are g₃ = 1 + 2λ₅. Thus, A, B, C, D, and G are M-matrices, and Equations (81), (82), (83), (84), and (85) are:

If we consider that S^m > 0, E^m > 0, I^m > 0, Q^m > 0, and V^m > 0, then the M-matrix along with (60) implies that the values of all of the state variables, i.e., S^m+1, E^m+1, I^m+1, Q^m+10, and V^m+1 are positive. Hence, the theorem is done by using the principle of mathematical induction.

This theorem is applied for drawing the conclusion that the proposed scheme, which is implicit in nature, guarantees the positive solution unconditionally.

7. Numerical Experiment and Simulations

A numerical test is performed on both the points of equilibrium for all the schemes under consideration. The set of parametric values considered for the test problem at disease-free equilibrium point [7] is given as:

b = 0.7, μ = 0.06, σ = 0.7, β = 0.35, γ = 0.15, ω = 0.3, q₁ = 0.2, q₂ = 0.2, q₃ = 0.1, ϕ = 0.4, ϵ = 0.05, α = 2, d₁ = 0.05, d₂ = 0.01, d₃ = 0.001, d₄ = 0.01, d₅ = 0.01.

For the endemic equilibrium point, the following parametric values are used:

b = 0.7, μ = 0.06, σ = 0.7, β = 0.35, γ = 0.15, ω = 0.06, q₁ = 0.2, q₂ = 0.2,

q₃ = 0.1, ϕ = 0.4, ϵ = 0.05, α = 2, d₁ = 0.05, d₂ = 0.01, d₃ = 0.001, d₄ = 0.01, d₅ = 0.01.

The initial condition for the model (6)–(10) is:

7.1. Disease-Free Equilibrium Point

In this section, graphs of all the state variables against time are presented (for DFE) to illustrate the behavior of the schemes. In Figures 1, 2, we consider h = 0.5, λ₁ = 0.3, λ₂ = 0.06, λ₃ = 0.006, λ₄ = 0.06, and λ₅ = 0.06.

Figure 1

Figure 2

Figures 1, 2 validate the preservation of the positivity property in both of the proposed operator splitting NSFD schemes. Also, all the graphs in Figures 1, 2 show that both proposed NSFD schemes achieve convergence to the DFE point. Next, we examine the behavior of forward Euler and backward Euler splitting schemes at different values of h and τ.

In parts (a) and (b) of Figure 3, we take the same values of h and τ as in Figures 1, 2 for the graphs of forward Euler operator splitting scheme and backward Euler operator splitting scheme. The graphs clearly show the failure of the positivity property of both classic schemes. In parts (c) and (d) of Figure 3, both existing splitting schemes converge to the false DFE equilibrium point for susceptible individuals.

Figure 3

In parts (a) and (b) of Figure 4, we take the same values of h and τ as given in parts (c) and (d) of Figure 3 for the explicit and implicit NSFD operator splitting schemes, respectively. The graphs clearly show that both of the proposed NSFD schemes not only preserve the positivity property but also achieve convergence to the true equilibrium point.

Figure 4

7.2. Endemic Equilibrium Point

In this section, we present simulations of the SEIQV epidemic model at the EE point using all of the operator splitting FD schemes. In Figures 5, 6, we consider h = 0.5, λ₁ = 0.3, λ₂ = 0.06, λ₃ = 0.006, λ₄ = 0.06, and λ₅ = 0.06.

Figure 5

Figure 6

Figures 5, 6 depict the graphs of susceptible, exposed, infected, quarantined, and vaccinated individuals for the EE point using the explicit operator splitting NSFD scheme and implicit operator splitting NSFD scheme, respectively. All the graphs in Figures 5, 6 demonstrate that both of the proposed operator splitting NSFD schemes preserve the property of positivity. These graphs also show that both proposed schemes converge to the EE point.

Again, both the forward and backward Euler FD schemes fail to preserve the positivity property and converge to the false EE point, as shown in Figure 7.

Figure 7

Figure 8 shows that the proposed NSFD operator splitting methods are consistent with the continuous reaction-diffusion system as they not only preserve the positivity property but converge to the EE point.

Figure 8

8. Conclusion

In this work, we consider the SEIQV reaction-diffusion epidemic model. The stability of the SEIQV model is guaranteed numerically by using criteria defined by Routh-Hurwitz. We also find the bifurcation value of the important vaccination parameter ω of SEIQV epidemic systems with diffusion and without diffusion. We design two novel and efficient operator splitting NSFD schemes for the SEIQV reaction-diffusion system. The NSFD schemes put forth, which are technically operator splitting schemes, possess the same behavior as is possessed by the SEIQV epidemic system. To conclude regarding the designed methods, we present two novel numerical schemes, one of which is explicit and the other of which is implicit in nature. The explicit scheme is more computationally efficient than the implicit scheme, but it has conditional stability while the implicit scheme is stable unconditionally. Both schemes employ structural splitting, due to which they deal adroitly with the nonlinearity of the reaction-diffusion system. These schemes avoid the false chaos that is a part of many existing methods. The positive solution of the SEIQV model is sustained by both schemes. Also, the nature of the stability of equilibria is preserved effectively by the proposed NSFD schemes. It is also shown that classical schemes, in parallel to our proposed schemes, produce chaos, leading to inconsistencies and instabilities numerically. The currently designed schemes are a valuable contribution for finding the solutions of nonlinear dynamical systems comprising differential equations. These NSFD schemes will become very efficient for the solution of one- and multi-dimensional reaction-diffusion population models, auto-catalytic chemical reaction models, and many more.

References

• 1.

  AhmedNWeiZBaleanuDRafiqMRehmanMA. Spatio-temporal numerical modeling of reaction-diffusion measles epidemic system. Chaos. (2019) 29:103101. 10.1063/1.5116807

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 2.

  AbroKAKhanINisarKS. Novel technique of Atangana and Baleanu for heat dissipation in transmission line of electrical circuit. Chaos Solitons Fractals. (2019) 129:40–510.1016/j.chaos.2019.08.001

  • CrossRef
  • Google Scholar

• 3.

  AhmedNRafiqMBaleanuDRehmanMA. Spatio-temporal numerical modeling of auto-catalytic Brusselator model. Roman J Phys. (2019) 64:110.

  • Google Scholar

• 4.

  GillVSinghJSinghY. Analytical solution of generalized space-time fractional advection-dispersion equation via coupling of sumudu and fourier transforms. Front Phys. (2019) 6:151. 10.3389/fphy.2018.00151

  • CrossRef
  • Google Scholar

• 5.

  KumarDSinghJAl QurashiMBaleanuD. A new fractional SIRS-SI malaria disease model with application of vaccines, antimalarial drugs, and spraying. Adv Differ Equ. (2019) 278. 10.1186/s13662-019-2199-9

  • CrossRef
  • Google Scholar

• 6.

  SaqibMKhanIShafieS. Application of fractional differential equations to heat transfer in hybrid nanofluid: modeling and solution via integral transforms. Adv Differ Equ. (2019) 52: 10.1186/s13662-019-1988-5

  • CrossRef
  • Google Scholar

• 7.

  LiuXYangL. Stability analysis of an SEIQV epidemic model with saturated incidence rate. Nonlinear Anal Real World Appl. (2012) 13:2671–9. 10.1016/j.nonrwa.2012.03.010

  • CrossRef
  • Google Scholar

• 8.

  MickensRE. Nonstandard finite difference models of differential equations world scientific (1994).

  • Google Scholar

• 9.

  AhmedNRafiqMRehmanMAIqbalMSAliM. Numerical modelling of three dimensional Brusselator reaction diffusion system. AIP Adv. (2019) 9:015205. 10.1063/1.5070093

  • CrossRef
  • Google Scholar

• 10.

  MickensRE. Dynamic consistency: a fundamental principle for constructing nonstandard finite difference schemes for differential equations. J Differ Equ Appl. (2005) 11:645–53. 10.1080/10236190412331334527.

  • CrossRef
  • Google Scholar

• 11.

  Macias-DiazJEPuriA. An explicit positivity-preserving finite-difference scheme for the classical Fisher-Kolmogorov-Petrovsky-Piscounov equation. Appl Math Comput. (2012) 218:5829–37. 10.1016/j.amc.2011.11.064

  • CrossRef
  • Google Scholar

• 12.

  MickensRE. A nonstandard finite difference scheme for a Fisher PDE having nonlinear diffusion. Comput Math Appl. (2003) 45:429–36. 10.1016/S0898-1221(03)80028-7

  • CrossRef
  • Google Scholar

• 13.

  MickensRE. A nonstandard finite difference scheme for an advection-reaction equation J Diff Eq Appl. (2004) 10:307–1312. 10.1080/10236190410001652838

  • CrossRef
  • Google Scholar

• 14.

  AhmedNRafiqMRehmanMAAliMAhmadMO. Numerical modeling of SEIR measles dynamics with diffusion. Commun Math Appl. (2018) 9:315–26.

  • Google Scholar

• 15.

  AhmedNShahidNIqbalZJawazMRafiqMTahiraSSet al. Numerical modeling of SEIQV epidemic model with saturated incidence rate. J Appl Environ Biol Sci. (2018) 8:67–82.

  • Google Scholar

• 16.

  AhmedNJawazMRafiqMRehmanMAAliMAhmadMO. Numerical treatment of an epidemic model with spatial diffusion. J Appl Environ Biol Sci. (2018) 8:17–29.

  • Google Scholar

• 17.

  ChinviriyasitSChinviriyasitW. Numerical modeling of SIR epidemic model with diffusion. Appl Math Comp. (2010) 216:395–409. 10.1016/j.amc.2010.01.028

  • CrossRef
  • Google Scholar

• 18.

  FatimaUAliMAhmedNRafiqM. Numerical modeling of susceptible latent breaking-out quarantine computer virus epidemic dynamics. Heliyon. (2018) 4:e00631. 10.1016/j.heliyon.2018.e00631

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19.

  JansenHTwizellEH. An unconditionally convergent discretization of the SEIR model. Math Comput Simulat. (2002) 58:147–58. 10.1016/S0378-4754(01)00356-1

  • CrossRef
  • Google Scholar

• 20.

  RafiqM. Numerical modeling of infectious diseases dynamics. Ph.D. thesis, Lahore, Pakistan: University of Engineering and Technology (2016).

  • Google Scholar

• 21.

  QinWWangLDingX. A non-standard finite difference method for a hepatitis B virus infection model with spatial diffusion. J Differ Equ Appl. (2014) 20:1641–51. 10.1080/10236198.2014.968565

  • CrossRef
  • Google Scholar

• 22.

  MannaKChakrabartySP. Global stability and a non-standard finite difference scheme for a diffusion driven HBV model with capsids. J Differ Equ Appl. (2015) 21:918–33. 10.1080/10236198.2015.1056524

  • CrossRef
  • Google Scholar

• 23.

  ChakrabrtyASinghMLucyBRidlandP. Predator-prey model with prey- taxis and diffusion. Math Comput Model. (2007) 46:482–98. 10.1016/j.mcm.2006.10.010

  • CrossRef
  • Google Scholar

• 24.

  IslamSHaiderN. Numerical solution of compartmental model by meshless finite difference methods. Appl Math Comput. (2014) 238:408–35. 10.1016/j.amc.2014.04.014

  • CrossRef
  • Google Scholar

• 25.

  HarwoodRC. (2011) Operator splitting method and applications for semilinear parabolic partial differential equations. Ph.D. dissertation. Pullman, WA: Dept. Math., Washington State Univ.

  • Google Scholar

• 26.

  HarwoodRCManoranjanVSEdwardsDB. Lead-acid battery model under discharge with a fast splitting method. IEEE Trans Energy Convers. (2011) 26:1109–17. 10.1109/TEC.2011.2162093

  • CrossRef
  • Google Scholar

• 27.

  YanenkoNN. The method of fractional steps. Berlin; Heidelberg: Springer-Verlag (1971).

  • Google Scholar

• 28.

  ZharnitskyV. Averaging for split-step scheme. Nonlinearity. (2003) 16:1359–66. 10.1088/0951-7715/16/4/310

  • CrossRef
  • Google Scholar

• 29.

  AnsarizadehFSinghMRichardsD. Modelling of tumor cells regression in response to chemotherapeutic treatment. Appl Math Model. (2017) 48:96–112. 10.1016/j.apm.2017.03.045

  • CrossRef
  • Google Scholar

• 30.

  NaheedASinghMLucyD. Numerical study of SARS epidemic model with the inclusion of diffusion in the system. Appl Math Comput. (2014) 229:480–98. 10.1016/j.amc.2013.12.062

  • CrossRef
  • Google Scholar

• 31.

  NaheedASinghMLucyD. Effect of treatment on transmission dynamics of SARS epidemic. J Infect Non Infect Dis. (2016) 2:16. 10.24966/INID-8654/100016

  • CrossRef
  • Google Scholar

• 32.

  SamsuzzohaMD. A Study on numerical solutions of epidemic models. Ph.D. thesis, Australia: Swinburne University of Technology (2012).

  • Google Scholar

• 33.

  WangH. Q., Numerical studies on the split-step finite difference method for nonlinear Schrödinger equations. Appl Math Comput. (2005) 170:17–35. 10.1016/j.amc.2004.10.066

  • CrossRef
  • Google Scholar

• 34.

  FujimotoTRanadeR. Two characterizations of inverse-positive matrices: the Hawkins-Simon condition and the Le Chatelier-Braun principle. Electron J Linear Algebra. (2004) 11:59–65. 10.13001/1081-3810.1122

  • CrossRef
  • Google Scholar