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

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.


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][2][3][4][5][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 βSI 1+αI , if number of infected individuals I becomes very large, βSI 1+αI approaches the saturation level. The infection force is measured by βI, which describes the penetration of the disease into a fully susceptible population. 1 1+αI 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 1 , q 2 , q 3 = 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 0 < 1, where R 0 is the reproductive number, which basically measures the occurrence of disease. The EE point exists if R 0 > 1. This implies that the SEIQV reaction-diffusion system (6-10) always converges to the DFE point or EE point if R 0 < 1 or R 0 > 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 reactiondiffusion 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][10][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 positivitypreserving numerical techniques for the approximate solution of epidemic models: see, for example [14][15][16][17][18][19][20][21][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.

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 0 is given as: R 0 is the reproductive value. Now, if R 0 < 1, the model acquires a DFE point, and if R 0 > 1, the model acquires an EE point.
The Routh-Hurwitz stability criterion gives: The Table 2 reflects the numerical stability of the equilibrium point against the various cases, as discussed in Table 1.

BIFURCATION VALUE OF VACCINATION PARAMETER ω INDEPENDENT OF DIFFUSION
Considering the vaccination parameter ω, to find its bifurcation value, a 11 , a 12 , are used instead of S * , E * , I * , Q * , and V * .  The Routh-Hurwitz criterion for stability gives: where the values of ξ 1 , ξ 2 , ξ 3 , ξ 4 , and ξ 5 are obtained from the expression of the characteristic equation (without diffusion) given in paper [24]. f 5 (ω) = 0 gives the value of bifurcation for ω. This value transfers the stability of a continuous model from stable to unstable. f 5 (ω) = 0 provides the bifurcation value ω = 0.228360507. The EE point is stable for ω less than ω = 0.228360507.

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 11 , a 12 , to a 11 = −0.6893640968ω − 0.2309369724, The Routh-Hurwitz criterion for stability gives: 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.

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][26][27][28][29][30][31][32][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).
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,

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 2 ) for all the methods under study.

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: which is the condition of stability for the explicit operatorsplitting 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.

NUMERICAL EXPERIMENT AND SIMULATIONS
A numerical test is performed on both the points of equilibrium for all the schemes under consideration.

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

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, λ 1 = 0.3, λ 2 = 0.06, λ 3 = 0.006, λ 4 = 0.06, and λ 5 = 0.06. 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 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.

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 reactiondiffusion population models, auto-catalytic chemical reaction models, and many more.

AUTHOR CONTRIBUTIONS
NA, MR, and MAR designed the study. NA, MF, and MR developed the methodology. NA and MF collected the data. NA, DB, KN, and IK performed the analysis. NA, MF, DB, MA, and IK wrote the manuscript.