# A SEIR Epidemic Model of Whooping Cough-Like Infections and Its Dynamically Consistent Approximation

Authors: M.M. Alqarni; Arooj Nasir; Maryam Ahmed Alyami; Ali Raza; Jan Awrejcewicz; Muhammad Rafiq; Nauman Ahmed; Tahira Sumbal Shaikh; Emad E. Mahmoud
Journal: Complexity
 (2022)

---

## Abstract

Whooping cough is a highly transmitted disease around the world. According to the World Health Organization (WHO), 0.15 million cases had reported globally in 2018. Most of the Asian and African states are infected regions. Through the study, we investigated the whole population into the four classes susceptible (S), exposed (E), infected (I), and vaccinated or recovered (R). The transmission dynamics of whooping cough disease are studied analytically and numerically. Analytical analyses are positivity, boundedness, reproduction number, equilibria, and local and global stabilities. In numerical analysis, we developed an implicit numerical integration scheme consistent with the biological problem’s properties. The analysis of the implicit method for the said model is dynamically consistent, positive, and bounded. Furthermore, an implicit numerical integration scheme is suitable for studying a particular epidemic model such as the whooping cough disease.


---

## Body

## 1. Introduction

Whooping cough is also known as pertussis. It is a highly contagious disease that is caused by bacteria. The symptoms of the common cold are such as runny nose, low-grade fever, and mild cough. After some time, these symptoms become worse. After a fit of coughing, a high pitch, whoop sound is produced as the patient breathes. This coughing sometimes can be severe. It can cause vomiting, breaking of the ribs, and fatigue in children. Difficulty in breathing may cause whooping cough in infants. After some time, symptoms of whooping cough become apparent. Pertussis is an infectious disease that is caused by bacteria. The name of this bacterium isBordetella pertussis. It is an airborne disease. When an infected person coughs or sneezes, the germs spread in the air through droplets. When a normal person breathes in the infected air, the germ enters their body through air. So, the normal person gets infected. There have been some uncertainties about pertussis being a zoonotic disease. Some scientists believed that some primates are highly vulnerable to Bordetella pertussis. When injected with a low dose of bacterium, primates can develop whooping cough. It is unclear if another wild animal can get infected, although gorillas can get infected by the disease. However, primates in zoo are usually vaccinated to prevent the spread of the disease. The symptoms that appear in the earlier stage of pertussis are almost similar to later-stage symptoms. The difference exists only in the degree of their severity. Symptoms that appear during the early stage are milder than that of the later stage. The initial symptoms resemble the symptoms of the common cold. In babies, the cough may be milder or absent altogether. However, whooping cough can be dangerous for babies, as half of the babies with whooping cough need intensive care. Early symptoms usually last 1 to 2 weeks and usually include runny nose, low fever, milder cough, and difficulty in breathing, especially in children; these are symptoms of whooping cough. Symptoms are not severe at the initial stage, so it is hard to diagnose early. It is diagnosed only when the conditions become much more severe. The symptoms that appear during the early stage of the disease become worse within 1 to 2 weeks. It includes fits of rapid cough, often followed by a high-pitched whoop sound, and vomiting, often caused by excessive coughing and fatigue. Pertussis can cause rapid and violent coughing again and again. When it occurs more frequently, it can decrease the presence of air in the lungs. It can cause vomiting and exhaustion in the patient. Coughing is more frequent at night. The condition of the patient is usually normal between the fits. The condition of the patients becomes worse with time. It can go on for ten weeks, and for that reason, whooping cough is also known as the 100-day cough in China. There have been worldwide periodic occurrences of whooping cough. It has infected a chunk of the population. According to a report in 2014, there were 24.1 million pertussis cases worldwide. About 160,700 children under five years died because of this disease. However, data from the developing world about pertussis are meager. Moreover, this rate is higher among infants. The mortality rate among infants with less than 12 months is 4%. It is 1% in children in the age group of 1–4 years. The risk of whooping cough can be decreased with the help of vaccination. There are two types of vaccines available right now. Whole-cell vaccines contain killed Bordetella pertussis, and the second is the acellular pertussis vaccine. It contains one or several highly purified pertussis antigens. In some countries, adults are given some additional pertussis vaccine to boost their immune systems. Pregnant women are also inoculated with doses to protect the infant who is too young to be vaccinated. The use of antibiotics has been useful in controlling the spread of the disease. It is also useful for the protection of already infected persons. The use of antibiotics such as erythromycin at the early stage can reduce the effects of pertussis. Still, if the conditions have already been worsening, the antibiotic does not change the victim’s condition but can further spread the disease. People who are already infected should not be allowed to contact other people, especially children and pregnant women. People who are not vaccinated until at least five days of infection are more in danger. Contact investigation is the management of the people who have been in contact with the infected people. These include a family member or people in direct contact with infected people. The people who work at medical facilities are also at greater risk because they are constantly contacting infected people. Moreover, for infection to take place, the contact does not have to be very close. Breathing in the air, sneezing, or coughing are enough to cause infection. Infants, pregnant women in their third trimester, and health workers are especially vulnerable to pertussis. People who have been in contact with the infected people and developed the symptoms should be tested. People who are exposed but do not develop the symptoms should not be tested, but they should take necessary precautions. The World Health Organization (WHO), one of the departments of the U.N., has launched several global programs to raise awareness among the masses about the relationship between whooping cough and influenza. This organization is helping the masses by providing them with information and treatment. It also provides the historical record and statistics about control over malaria progress in different parts of the world, where malaria has become an epidemic. The standard reports of the U.N. department have shown that out of 99 countries where whooping cough is outspread, 80 countries are now in the control phase and 19 are in the pre-elimination phase. The following are the uses of mathematical modeling in real-life problems: dynamical models are given as reinforcing mathematics previously learned, determining new mathematical affiliations, joining mathematics to the world, elasticity students and professionals’ problem-solving skills, and linking mathematics to other disciplines such as an epidemic disease. Mathematical modeling is a procedure in which we use mathematical and biological structures, graphs, equations, and diagrams to analyze real-life problems and world situations, or it is the procedure of making mathematical relations by using some conditions. Mathematical modeling is a remarkable ground to construct models of disease spreading. Mathematical modeling of infectious disease is an instrument in which we see how disease transmits individually and control it. In 2019, Ameri and Cooper devised a model to discuss the transmission of whooping cough in Nebraska [1]. In 2019, Ali et al. applied numerical techniques to the SIR epidemic model with a time delay effect [2]. In 2019, Suleman and Riaz observed the whooping cough model with optimally controlled vaccination [3]. In 2018, Althouse and Scarpino studied the epidemiological dynamics of the whooping cough model and investigated its evolutionary and ecological signatures [4]. In 2018, Bhattacharyya et al. explained the periodicity of whooping cough dynamics using species interactions [5]. In 2018, Van et al. studied the cost effectiveness of maternal vaccination for whooping cough in Australia [6]. In 2017, Obaidat et al. solved the second-order whooping cough model by implementing an effective numerical method [7]. In 2013, Obaidat et al. investigated about the two-dimensional diffusion type whooping cough model [8]. In 2010, Black and McKane quantified whooping cough dynamics using the well-known assumption of modeling [9]. In 2010, Arenas et al. applied a predictor-corrector scheme to various epidemic models [10]. Hathout et al. [11] studied the triple age-structured model, including the protection strategy. Bentout et al. [12] studied a bifurcation analysis for an epidemic model by considering two delayed factors. George et al. [13] studied the bifurcation analysis of the discrete-time model. In 2005, Keeling et al. discussed the network of various epidemic models [14]. In 1998, Duncan et al. studied the factors affecting whooping coughs such as malnutrition and dense population [15]. In 1996, Duncan et al. studied the whooping cough epidemics in London and factors such as malnutrition, dynamics of infection, and seasonal forcing [16]. In 1986, Knox et al. presented a model for the control strategies of the whooping cough disease [17]. Two types of the latest epidemic models are studied in references [18, 19]. Mezouaghi et al. [20] studied the global properties of delayed epidemic models. The unreported cases and their modeling regarding COVID-19 are studied in reference [21]. Well-known models studied with stochastic techniques are presented in references [22–25].The main focus of this research is to study the efficiencies of the implicit numerical integration scheme for the whooping cough disease. The results produced by the implicit numerical integration scheme are the cross-check of analytical results. We claim that an implicit numerical integration scheme tool is reliable, efficient, and adopts all the dynamical properties for long-time disease behavior.The flow of the paper is based on the following sections. In section2, the formulation of the deterministic whooping cough epidemic model with fundamental properties of positivity, boundedness, local stability, and implementation of an implicit numerical integration scheme. Finally, in Section 3, the conclusion is presented.


## 2. Model Formulation

The population is divided into four components of subpopulations to study the dynamics of disease. The people will be represented byN:0,∞⟶ℛ with the time t≥0. Furthermore, each of the components of the population will be denoted by a nonnegative differentiable function S,E,I,R:0,∞⟶ℛ. For any time, St represents the susceptible humans, Et represents the exposed humans, I(t) defines the infected humans, and R(t) means the recovered humans. The graphical abstract of whooping cough disease is presented in Figure 1.Figure 1
Flowchart of the whooping cough model.The physical interpretation of the model is presented as follows:μ represents the birth and death rate, β represents the transmission rate of whooping cough from susceptible to exposed humans, V represents the rate at which revealed humans become infected, and α is the rate at which infected humans become recovered after vaccination or natural. The model of the differential equations is as follows:(1)S′t=μ−μSt−βStIt,t≥0,(2)E′t=βStIt−VEt−μEt,t≥0,(3)I′t=VEt−αIt−μIt,t≥0,(4)R′t=αIt−μRt,t≥0,where S0≥0,E0≥0,I0≥0,R0≥0, and St+Et+It+Rt≤N


### 2.1. Model Properties

The feasible regionℂ=S,E,I,Rεℝ+4:St+Et+It+Rt≤N,S≥0,E≥0,I≥0,R≥0, at any time t≥0, and the solution of the model remains positive and bounded.Property 1.
For positivity, considering the system (1)–(4), as follows:(5)dSdt|S=0=μ≥0,dEdt|E=0=βSI≥0,dIdt|I=0=VE≥0,dRdt|R=0=αI≥0,as desired.Property 2.
For boundedness, considering the population function at any timet as follows:(6)Nt=St+Et+It+Rt.
Taking a change with respect to time on both sides, we have as follows.
dN/dt=dS/dt+dE/dt+dI/dt+dR/dt, by using the system (1)–(4), we have as follows:(7)dNdt=μ−μS+E+I+R,
dN/dt≤μ−μN. Since, St+Et+It+Rt≤N
From Gronwall’s inequality, we obtain(8)Nt≤e−μτ+1,t≥0.
limt⟶∞SupNt≤1, as desired.Property 3.
The system (1)–(4) admits two types of equilibria and lies in the region ℂ.Proof.
By assuming the variables as a constant into the system’s (1)–(4), we have as follows:(9)μ−μS−βSI=0,μ−μS−βSI=0,(10)βSI−VE−μE=0,(11)VE−αI−μI=0,(12)αI−μR=0.
From equations (9)–(12), the cough existing equilibrium (CEE-C2=S∗,E∗,I∗,R∗, we have as follows:(13)S∗=αV+μV−μα−μ2βV,E∗=μαV+μV−μα−μ2−βVβVμ−V,I∗=μαV+μV−μα−μ2−βVβμ−Vα+μ,R∗=ααV+μV−μα−μ2−βVTβμ−Vα+μ.
Also, the cough-free equilibrium (CFE-C1 of the system (28)–(31) is as follows:(14)S≠0,E=0,I=0,R=0.
Thus, the cough-free equilibrium (CFE-C1=S1,E1,I1,R1=1,0,0,0.


### 2.2. Local Stability

Consider the system (28)–(31) as a function of F,G,H,andK as follows:(15)F=μ−μs−βSI,G=βSI−VE−μE,H=VE−αI−μI,K=αI−μR.The partial derivates are as follows:(16)∂F∂S=−μ−βI,∂F∂E=0,∂FdI=−βS,∂F∂R=0,∂GdS=βI,∂G∂E=–V−µ,∂G∂I=βS,∂G∂R=0,∂H∂S=0,∂H∂E=V,∂H∂I=−α−μ,∂H∂ω=0,∂K∂S=0,∂K∂E=0,∂K∂I=α,∂K∂R=−μ.Thus, the Jacobean matrix of the model is as follows:(17)J=−μ−βI0−βS0βI−V−μβS00V−α−μ000α−μ.Theorem 1.
The cough-free equilibrium (CFE-C1=S1,E1,I1,R1=1,0,0,0 is locally asymptotically stable if R0<1 and unstable if R0>1.Proof.
The Jacobean matrix evaluated atC1 is as follows:(18)JC1=−μ0−β00−V−μβ00v−α−μ000α−μ,JC1−λI=−μ−λ0−β00−V−μ−λβ00V−α−μ−λ000α−μ−λ=0,λ1=−μ<0,λ2=−μ<0,JC1−λI=−V−μ−λβV−α−μ−λ=0,λ2+λV+2μ+α+αV+μV+αμ+μ2−βV=0,λ2+a1λ+a0=0.
By using the properties of Routh Hurwitz criteria of 2nd order, the coefficients of 2nd order polynomials should be positive. So,wherea1=V+2μ+α>0 and a0=αV+μV+αμ+μ2−βV>0, if. R0<1
IfαV+μV+αμ+μ2−βV>0(19)αV+μV+αμ+μ2>βV,1>βVαV+μV+αμ+μ2,1>βVα+μV+μ,R0=βVα+μV+μ<1.
Hence, using the Routh Hurwitz criterion of 2nd order,C1 is locally asymptotically stable.Theorem 2.
The cough existing equilibriumCEE−C2=S∗,E∗,I∗,R∗ is locally asymptotically stable if R0>1.Proof.
The Jacobean matrix evaluated atC2 is as follows:(20)JC2=−μ−βI∗0−βS∗0βI∗−V−μβS∗00V−α−μ000α−μ,JC2−λI=−μ−βI∗−λ0−βS∗0βI∗−V−μ−λβS∗00V−α−μ−λ000α−μ−λ=0,λ1=−μ<0,λ3+3μ+βI∗+V+αλ2+2Vμ+3μ2+2αμ+VβI∗+2μβI∗+αβI∗+Vα−VβS∗λ+Vαμ+Vμ2+αμ2+μ3−VμβS∗+VαβI∗+VμβI∗+αβI∗μ+μ2βI∗=0,λ3+a2λ2+a1λ+a0=0,where a2=3μ+βI∗+V+α,a1=2Vμ+3μ2+2αμ+VβI∗+2μβI∗+αβI∗+Vα−VβS∗, a0=Vαμ+Vμ2+αμ2+μ3−VμβS∗+VαβI∗+VμβI∗+αβI∗μ+μ2βI∗
Hence, by using the Routh Hurwitz criterion of 3rd order,(21)a2>0,a0>0,a2a1>a0,ifR0>1.
Thus,C2 is locally asymptotically stable.


### 2.3. Global Stability

Theorem 3.
(Global stability atC1 The system (1) at CFE−C1=S1,E1,I1,R1=1,0,0,0 is globally asymptotically stable if R0<1.Proof.
By letting the Lyapunov function,U:Ω⟶R is defined as follows:(22)U=S−S1−S1logSS1+E−E1−E1logEE1+I−I1−I1logII1+R−R1−R1logRR1,dUdt=1−S1SdSdt+dEdt+dIdt+dRdt,dUdt=1−S1Sμ−μS−βSI+βSI−μE−μI−μR,dUEdt=S−S1μS−βI−μ+βSI−μE−μI−μR,dUdt=−μS−S12SS1,dUdt≤0,ifR0<1.
Hence, the system is globally asymptotically stable atC1.Theorem 4.
(Global stability atC2) The system (1) at C2=S∗,E∗,I∗,R∗ is globally asymptotically stable if R0>1.Proof.
By letting the Lyapunov function,Z:Ω⟶R is defined as follows:(23)Z=K1S−S∗−S∗logSS∗+K2E−E∗−E∗logIE∗+K3I−I∗−I∗logII∗+K4R−R∗−R∗logRR∗.where Kii=1,2,3,4 are positive constants to be chosen later. We have as follows:(24)dZdt=K1s−s∗sμ−μS−βSI+K2E−E∗EβSI−VE−μE+K3I−I∗IVE−αI−μI,K4R−R∗RαI−μR,dZdt=K1s−s∗μs−βI−μ+K2E−E∗βSIE−V−μ+K3I−I∗VEI−α−μ+K4R−R∗αIR−μ.
If we chooseKii=1,2,3,4, then(25)dZdt=μS−S∗2ss∗−βS−S∗I−I∗−βSII−I∗2EE∗−VEI−I∗2II∗−αIR−R∗2RR∗where dZ/dt≤0, for R0>1, and dZ/dt=0 only if S=S∗,I=I∗,E=E∗,R=R∗.
Hence, by Lasalle’s invariance principle,C2 is globally asymptotically stable (GAS) in Ω.


### 2.4. Threshold Number

By considering the infectious classes, to find the transmission and transition matrices with the help next-generation matrix method by substituting the cough-free equilibrium as follows:(26)E′I′R′=0β0000000EIR−V+μ00−Vα+μ00−αμEIRwhere G1 is the transmission matrix and G2 is the transition matrix as follows:(27)G1=0β0000000,G2=V+μ00−Vα+μ00−αμ,G2−1=1V+μα+μμαμ+μ200VμμV+μ20αVαV+αμαV+αμ+μ2+Vμ,G1G2−1=1V+μα+μμ0β0000000αμ+μ200VμμV+μ20αVαV+αμαV+αμ+μ2+Vμ,G1G2−1−λI=βVV+μα+μ−λβα+μ00−λ000−λ=0,βVV+μα+μ−λλ2=0.Hence, the dominant eigenvalue of theG1G2−1 is called the threshold number and is denoted as R0=βV/V+μα+μ.Furthermore, the scientific literature is present in Table1, for the simulating behavior of the system (1–4) at both equilibria of the model as follows:Table 1
Values of parameters [26].
Parameters%C1%C2μ0.50.5β22.5V0.80.8α0.90.9


### 2.5. Explicit Methods

Explicit methods have a significant role in handling physical systems in all disciplines. In the mathematical modeling of infectious diseases, explicit methods are used for obtaining numerical approximations to the solutions of time-dependent ordinary and partial differential equations, as required in computer simulations of physical processes. Explicit methods calculate the state of a system at a later time from the state of the system at the current time. The explicit methods are easier to program and can be calculated within a shorter time. However, its stability is so low that the small step size is used to prevent divergence.


### 2.6. Implicit Methods

Forn≥0 and h is represented by the discretization parameter; then, the discretization of the system (1)–(4) under the assumption of the implicit method in the particular case of the SEIR model is as follows:From equation (1),(28)Sn+1=Sn+hμ−μSn+1−βSn+1In,Sn+1=Sn+hμ1+hμ+hβIn.From equation (2),(29)En+1=En+hβSnIn−VEn+1−μEn+1,En+1=En+hβSnIn1+hV+hμ.From equation (3),(30)In+1=In+hVEn−αIn+1−μIn+1,In+1=In+hVEn1+hα+hμ.From equation (4),(31)Rn+1=Rn+hαIn+μRn+1,Rn+1=Rn+hαIn1+hμ.


#### 2.6.1. Convergence Analysis

Theorem 5.
The proposed implicit numerical integration scheme is stable if it satisfies the assumption that the system’s Eigenvalues lie in the unit circle [27, 28].Proof.
Considering the functionsL, M, N, and P from the system (28)–(31) as follows:(32)L=S+hμ1+hμ+hβI,M=E+hβSI1+hV+hμ,P=I+hVE1+hα+hμ,Q=R+hαI1+hμ.
If the all eigenvalues of the Jacobean in formula (33),(33)Ј=∂L∂S∂L∂E∂L∂I∂L∂R∂M∂S∂M∂E∂M∂I∂M∂R∂P∂S∂Q∂S∂P∂E∂Q∂E∂P∂I∂P∂R∂Q∂I∂Q∂R,which satisfies the condition in (34), at the equilibria of the model(34)ρJ<1.
The elements of the Jacobean matrix are as follows:(35)∂L∂S=11+hμ+hβI,∂L∂E=0,∂L∂I=−hβS+hμ1+hμ+hβI2,∂L∂R=0,∂M∂S=hβI1+hV+hμ,∂M∂E=11+hV+hμ,∂M∂I=hβS1+hV+hμ,∂M∂R=0,∂P∂S=0,∂P∂E=hV1+hα+hμ,∂P∂I=11+hα+hμ,∂P∂R=0,∂Q∂S=0,∂Q∂E=0,∂Q∂I=hα1+hμ,∂Q∂R=11+hμ.
The Jacobean cough-free equilibrium (CFE-C1=S1,E1,I1,R1=1,0,0,0.(36)J1,0,0,0=11+hμ0−hβ1+hμ0011+hV+hμhβ1+hV+hμ00hV1+hα+hμ11+hα+hμ000hα1+hμ11+hμ,
The eigenvalues of the Jacobean matrix are as follows.
λ1=1/1+hμ<1,λ2=1/1+hμ≪1 provided that R0<1(37)J−λI=11+hV+hμhβ1+hV+hμhV1+hα+hμ11+hα+hμ,
A= Trace of J, B= det of J(38)A=2+hV+α+2μ1+hV+hμ1+hα+hμ,B=1−h2βV1+hV+hμ1+hα+hμ.
The Jacobean matrix at the cough existing equilibriumC2=S∗,E∗,I∗,R∗ is as follows:(39)JS∗,E∗,I∗,R∗=11+hμ+hβI∗0−hβS∗+hμ1+hμ+hNβI∗20hβI∗1+hV+hμ11+hV+hμhβS∗1+hV+hμ00hV1+hα+hμ11+hα+hμ000hα1+hμ11+hμ.
For equilibriaC2, Figure 2 depicts the behavior of the largest eigenvalues and the values of parameters are presented in Table 1, as desired.Figure 2
Largest Eigenvalue of the system at the equilibriaC2.


#### 2.6.2. Positivity

Positivity is one of the most important physical traits of the state variables involved in the compartmental epidemic models. This property is investigated and ensured by the principle of mathematical induction for an implicit numerical integration scheme.Theorem 6.
Let the state variablesSt,Et,It, and Rt involved in the scheme are positive at t=0; furthermore, if all the parameters are also positive, then Sn+1≥0,En+1≥0,In+1≥0 and Rn+1≥0Proof.
Consider the equations (28)–(31) and apply the principle of mathematical induction, we proceed as follows:(40)Sn+1=Sn+hμ1+hμ+hNβInEn+1=En+hβSnIn1+hV+hμIn+1=In+hVEn1+hα+hμRn+1=Rn+hαIn1+hμ
At first, we putn=0 in the abovementioned system of equations and obtain the following results:(41)S1=S0+hμ1+hμ+hNβI0≥0.
Similarly,(42)E1=E0+hβS0I01+hV+hμ≥0,I1=I0+hVE01+hα+hμ≥0,R1=R0+hαI01+hμ≥0.
Now, forn=1, we arrive at the following expressions:(43)S2=S1+hμ1+hμ+hNβI1≥0.
Similarly,(44)E2=E1+hβS1I11+hV+hμ≥0,I2=I1+hVE11+hα+hμ≥0,R2=R1+hαI11+hμ≥0.
Next, suppose that the given set of equations guarantees the positivity property for the value of.n=2,3,4,…,n−1,i.e.,Sn≥0;En≥0;In≥0; and Rn≥0 for n=2,3,4,…,n−1 and state variables.
Now, the positivity will be investigated for an arbitrary positivity integern.
Forn∈Z, we observe that.(45)Sn+1=Sn+hμ1+hμ+hNβIn≥0,En+1=En+hβSnIn1+hV+hμ≥0,In+1=In+hVEn1+hα+hμ≥0,Rn+1=Rn+hαIn1+hμ≥0.
Hence, the projected scheme ensures the positivity for the state variablesSt,Et,It, and Rt for all the positive integer values of n.


#### 2.6.3. Boundedness

Theorem 7.
LetS0,E0,I0, and R0 are finite, such that S0+E0+I0+R0≤1. Also μ,β,V,α are all positive in the model. Then, the discretized state variables, Sn+1,En+1,In+1, and Rn+1 are bounded by recrossing defined real constant bn+1 such that Sn+1,En+1,In+1 and Rn+1<bn+1∀n∈Z+where bn+1=4bn+hμ+hβbn2+hVbn+hαbn and b1=4+hμ+hβS0I0+hVE0+hαI0.Proof.
Considering the equations of an implicit numerical integration scheme for the state variablesS,E,I, and R as,(46)Sn+11+hμ+hNβIn=Sn+hμEn+11+hV+hμ=En+hβSnInIn+11+hα+hμ=In+hVEnRn+11+hμ=Rn+hαIn
We obtain the following expression by adding all the equations in the system (46).(47)Sn+1+En+1+In+1+Rn+1+Sn+1hμ+hNβIn+En+1hV+hμ+In+1hα+hμ+Rn+1hμ=Sn+En+In+Rn+hμ+hβSnIn+hVEn+hαIn.
We prove the boundedness of the nonstandard numerical scheme by using the theory of mathematical induction and the constraints imposed on the parameters and state variables.
By substitutingn=0, we obtain as follows:(48)S1+E1+I1+R1+S1hμ+hNβI0+E1hV+hμ+I1hα+hμ+R1hμ=S0+E0+I0+R0+hμ+hβS0I0+hVE0+hαI0,S1+E1+I1+R1+S1hμ+hNβI0+E1hV+hμ+I1hα+hμ+R1hμ<4+hμ+hβS0I0+hVE0+hαI0=b1,⇒S1+E1+I1+R1+S1hμ+hNβI0+E1hV+hμ+I1hα+hμ+R1hμ≤b1,⇒S11+hμ+hNβI0+E11+hV+hμ+I11+hα+hμ+R11+hμ≤b1,⇒S1<b1,E1<b1,I1<b1,R1<b1.
Now, by substitutingn=1 in (47) to reach the following step(49)S21+hμ+hNβI1+E21+hV+hμ+I21+hα+hμ+R21+hμ=S1+E1+I1+R1+hμ+hβS1I1+hVE1+hαI1,⇒S21+hμ+hNβI1+E21+hV+hμ+I21+hα+hμ+R21+hμ<4b1+hμ+hβb1+hVb1+hαb1=b2,⇒S2<b2,E2<b2,I2<b2,R2<b2.
Next, consider that the expression (47) is bounded for positive integers. 2,3,4,⋅⋅⋅,n−1i.e., Sn<bn,En<bn,In<bn and Rn<bn∀n∈2,3,4,…,n−1
Now, we investigate the boundedness of a positivity integern.(50)Sn+11+hμ+hNβIn+En+11+hV+hμ+In+11+hα+hμ+Rn+11+hμ=S1+E1+I1+R1+hμ+hβSnIn+hVEn+hαIn,Sn+11+hμ+hNβIn+En+11+hV+hμ+In+11+hα+hμ+Rn+11+hμ<4bn+hμ+hβbn2+hVbn+hαbn=bn+1.
So,Sn+1<bn+1,En+1<bn+1,In+1<bn+1 and. Rn+1<bn+1,where. bn+1=4bn+hμ+hβbn2+hVbn+hαbnb1=4+hμ+hβS0I0+hVE0+hαI0
Hence,Sn+1,En+1,In+1, and Rn+1 are bounded by a real number bn+1 for all positive values of n.


#### 2.6.4. Consistency Analysis

In this section, the consistency analysis of an implicit numerical integration scheme is performed by using Taylor’s series expansion. First, we take equation (28) of the implicit numerical integration scheme and apply Taylor’s series expansion of Sn+1(51)Sn+1=Sn+hdSdt+h22!d2Sdt2+h33!d3Sdt3+⋅⋅⋅In the following expression:(52)Sn+11+hμ+hNβIn=Sn+hμ,Sn+hdSdt+h22!d2Sdt2+h33!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=Sn+hμ,Sn+hμSn+hNβInSn+hdSdt+hμhdSdt+hNβInhdSdt+h22!d2Sdt2+h33!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=Sn+hμ,hμSn+NβInSn+dSdt+hμdSdt+NβInhdSdt+h2!d2Sdt2+h23!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=hμ.We obtain as follows:(53)μSn+NβInSn+dSdt+hμdSdt+NβInhdSdt+h12!d2Sdt2+h3!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=μBy applyingh⟶0, we obtain as follows:(54)μS+NβIS+dSdt=μ,⇒dSdt=μ−μS−NβIS.This result implies that our discretized equation is consistent with equation (1) of the ODE model. Similarly, we take equation (29) and apply Taylor’s series expansion of En+1.(55)En+1=En+hdEdt+h22!d2Edt2+h33!d3Edt3+⋅⋅⋅In the following expression:(56)En+11+hV+hμ=En+hβSnIn.We obtain as follows:(57)VEn+μEn+dEdt+hVdEdt+hdEdtμ+h12!d2Edt2+h3!d3Edt3+⋅⋅⋅1+hV+hμ=βSnIn.Applyingh⟶0, we obtain as follows:(58)dEdt=βSI−VE−μE.Similarly, taking equation (30) and applying Taylor’s series expansion of In+1. After proper simplifications, we obtain as follows:(59)Inα+μIn+dIdt+hαdIdt+hμdIdtμ+h12!d2Idt2+h3!d3Idt3+⋅⋅⋅1+hα+hμ=VEn.Applyingh⟶0, we obtain as follows:(60)dIdt=VE−αI−μI.Similarly, taking equation (31) and applying Taylor’s series expansion of Rn+1, we obtain as follows:(61)μRn+dRdt+hμdRdt+h12!d2Rdt2+h3!d3Rdt3+⋅⋅⋅1+hμ=αIn.Applyingh⟶0, we obtain as follows:(62)dRdt=αI−μR.Hence, our discretized implicit numerical integration scheme is consistent with the ODES system (1–4).


## 2.1. Model Properties

The feasible regionℂ=S,E,I,Rεℝ+4:St+Et+It+Rt≤N,S≥0,E≥0,I≥0,R≥0, at any time t≥0, and the solution of the model remains positive and bounded.Property 1.
For positivity, considering the system (1)–(4), as follows:(5)dSdt|S=0=μ≥0,dEdt|E=0=βSI≥0,dIdt|I=0=VE≥0,dRdt|R=0=αI≥0,as desired.Property 2.
For boundedness, considering the population function at any timet as follows:(6)Nt=St+Et+It+Rt.
Taking a change with respect to time on both sides, we have as follows.
dN/dt=dS/dt+dE/dt+dI/dt+dR/dt, by using the system (1)–(4), we have as follows:(7)dNdt=μ−μS+E+I+R,
dN/dt≤μ−μN. Since, St+Et+It+Rt≤N
From Gronwall’s inequality, we obtain(8)Nt≤e−μτ+1,t≥0.
limt⟶∞SupNt≤1, as desired.Property 3.
The system (1)–(4) admits two types of equilibria and lies in the region ℂ.Proof.
By assuming the variables as a constant into the system’s (1)–(4), we have as follows:(9)μ−μS−βSI=0,μ−μS−βSI=0,(10)βSI−VE−μE=0,(11)VE−αI−μI=0,(12)αI−μR=0.
From equations (9)–(12), the cough existing equilibrium (CEE-C2=S∗,E∗,I∗,R∗, we have as follows:(13)S∗=αV+μV−μα−μ2βV,E∗=μαV+μV−μα−μ2−βVβVμ−V,I∗=μαV+μV−μα−μ2−βVβμ−Vα+μ,R∗=ααV+μV−μα−μ2−βVTβμ−Vα+μ.
Also, the cough-free equilibrium (CFE-C1 of the system (28)–(31) is as follows:(14)S≠0,E=0,I=0,R=0.
Thus, the cough-free equilibrium (CFE-C1=S1,E1,I1,R1=1,0,0,0.


## 2.2. Local Stability

Consider the system (28)–(31) as a function of F,G,H,andK as follows:(15)F=μ−μs−βSI,G=βSI−VE−μE,H=VE−αI−μI,K=αI−μR.The partial derivates are as follows:(16)∂F∂S=−μ−βI,∂F∂E=0,∂FdI=−βS,∂F∂R=0,∂GdS=βI,∂G∂E=–V−µ,∂G∂I=βS,∂G∂R=0,∂H∂S=0,∂H∂E=V,∂H∂I=−α−μ,∂H∂ω=0,∂K∂S=0,∂K∂E=0,∂K∂I=α,∂K∂R=−μ.Thus, the Jacobean matrix of the model is as follows:(17)J=−μ−βI0−βS0βI−V−μβS00V−α−μ000α−μ.Theorem 1.
The cough-free equilibrium (CFE-C1=S1,E1,I1,R1=1,0,0,0 is locally asymptotically stable if R0<1 and unstable if R0>1.Proof.
The Jacobean matrix evaluated atC1 is as follows:(18)JC1=−μ0−β00−V−μβ00v−α−μ000α−μ,JC1−λI=−μ−λ0−β00−V−μ−λβ00V−α−μ−λ000α−μ−λ=0,λ1=−μ<0,λ2=−μ<0,JC1−λI=−V−μ−λβV−α−μ−λ=0,λ2+λV+2μ+α+αV+μV+αμ+μ2−βV=0,λ2+a1λ+a0=0.
By using the properties of Routh Hurwitz criteria of 2nd order, the coefficients of 2nd order polynomials should be positive. So,wherea1=V+2μ+α>0 and a0=αV+μV+αμ+μ2−βV>0, if. R0<1
IfαV+μV+αμ+μ2−βV>0(19)αV+μV+αμ+μ2>βV,1>βVαV+μV+αμ+μ2,1>βVα+μV+μ,R0=βVα+μV+μ<1.
Hence, using the Routh Hurwitz criterion of 2nd order,C1 is locally asymptotically stable.Theorem 2.
The cough existing equilibriumCEE−C2=S∗,E∗,I∗,R∗ is locally asymptotically stable if R0>1.Proof.
The Jacobean matrix evaluated atC2 is as follows:(20)JC2=−μ−βI∗0−βS∗0βI∗−V−μβS∗00V−α−μ000α−μ,JC2−λI=−μ−βI∗−λ0−βS∗0βI∗−V−μ−λβS∗00V−α−μ−λ000α−μ−λ=0,λ1=−μ<0,λ3+3μ+βI∗+V+αλ2+2Vμ+3μ2+2αμ+VβI∗+2μβI∗+αβI∗+Vα−VβS∗λ+Vαμ+Vμ2+αμ2+μ3−VμβS∗+VαβI∗+VμβI∗+αβI∗μ+μ2βI∗=0,λ3+a2λ2+a1λ+a0=0,where a2=3μ+βI∗+V+α,a1=2Vμ+3μ2+2αμ+VβI∗+2μβI∗+αβI∗+Vα−VβS∗, a0=Vαμ+Vμ2+αμ2+μ3−VμβS∗+VαβI∗+VμβI∗+αβI∗μ+μ2βI∗
Hence, by using the Routh Hurwitz criterion of 3rd order,(21)a2>0,a0>0,a2a1>a0,ifR0>1.
Thus,C2 is locally asymptotically stable.


## 2.3. Global Stability

Theorem 3.
(Global stability atC1 The system (1) at CFE−C1=S1,E1,I1,R1=1,0,0,0 is globally asymptotically stable if R0<1.Proof.
By letting the Lyapunov function,U:Ω⟶R is defined as follows:(22)U=S−S1−S1logSS1+E−E1−E1logEE1+I−I1−I1logII1+R−R1−R1logRR1,dUdt=1−S1SdSdt+dEdt+dIdt+dRdt,dUdt=1−S1Sμ−μS−βSI+βSI−μE−μI−μR,dUEdt=S−S1μS−βI−μ+βSI−μE−μI−μR,dUdt=−μS−S12SS1,dUdt≤0,ifR0<1.
Hence, the system is globally asymptotically stable atC1.Theorem 4.
(Global stability atC2) The system (1) at C2=S∗,E∗,I∗,R∗ is globally asymptotically stable if R0>1.Proof.
By letting the Lyapunov function,Z:Ω⟶R is defined as follows:(23)Z=K1S−S∗−S∗logSS∗+K2E−E∗−E∗logIE∗+K3I−I∗−I∗logII∗+K4R−R∗−R∗logRR∗.where Kii=1,2,3,4 are positive constants to be chosen later. We have as follows:(24)dZdt=K1s−s∗sμ−μS−βSI+K2E−E∗EβSI−VE−μE+K3I−I∗IVE−αI−μI,K4R−R∗RαI−μR,dZdt=K1s−s∗μs−βI−μ+K2E−E∗βSIE−V−μ+K3I−I∗VEI−α−μ+K4R−R∗αIR−μ.
If we chooseKii=1,2,3,4, then(25)dZdt=μS−S∗2ss∗−βS−S∗I−I∗−βSII−I∗2EE∗−VEI−I∗2II∗−αIR−R∗2RR∗where dZ/dt≤0, for R0>1, and dZ/dt=0 only if S=S∗,I=I∗,E=E∗,R=R∗.
Hence, by Lasalle’s invariance principle,C2 is globally asymptotically stable (GAS) in Ω.


## 2.4. Threshold Number

By considering the infectious classes, to find the transmission and transition matrices with the help next-generation matrix method by substituting the cough-free equilibrium as follows:(26)E′I′R′=0β0000000EIR−V+μ00−Vα+μ00−αμEIRwhere G1 is the transmission matrix and G2 is the transition matrix as follows:(27)G1=0β0000000,G2=V+μ00−Vα+μ00−αμ,G2−1=1V+μα+μμαμ+μ200VμμV+μ20αVαV+αμαV+αμ+μ2+Vμ,G1G2−1=1V+μα+μμ0β0000000αμ+μ200VμμV+μ20αVαV+αμαV+αμ+μ2+Vμ,G1G2−1−λI=βVV+μα+μ−λβα+μ00−λ000−λ=0,βVV+μα+μ−λλ2=0.Hence, the dominant eigenvalue of theG1G2−1 is called the threshold number and is denoted as R0=βV/V+μα+μ.Furthermore, the scientific literature is present in Table1, for the simulating behavior of the system (1–4) at both equilibria of the model as follows:Table 1
Values of parameters [26].
Parameters%C1%C2μ0.50.5β22.5V0.80.8α0.90.9


## 2.5. Explicit Methods

Explicit methods have a significant role in handling physical systems in all disciplines. In the mathematical modeling of infectious diseases, explicit methods are used for obtaining numerical approximations to the solutions of time-dependent ordinary and partial differential equations, as required in computer simulations of physical processes. Explicit methods calculate the state of a system at a later time from the state of the system at the current time. The explicit methods are easier to program and can be calculated within a shorter time. However, its stability is so low that the small step size is used to prevent divergence.


## 2.6. Implicit Methods

Forn≥0 and h is represented by the discretization parameter; then, the discretization of the system (1)–(4) under the assumption of the implicit method in the particular case of the SEIR model is as follows:From equation (1),(28)Sn+1=Sn+hμ−μSn+1−βSn+1In,Sn+1=Sn+hμ1+hμ+hβIn.From equation (2),(29)En+1=En+hβSnIn−VEn+1−μEn+1,En+1=En+hβSnIn1+hV+hμ.From equation (3),(30)In+1=In+hVEn−αIn+1−μIn+1,In+1=In+hVEn1+hα+hμ.From equation (4),(31)Rn+1=Rn+hαIn+μRn+1,Rn+1=Rn+hαIn1+hμ.


### 2.6.1. Convergence Analysis

Theorem 5.
The proposed implicit numerical integration scheme is stable if it satisfies the assumption that the system’s Eigenvalues lie in the unit circle [27, 28].Proof.
Considering the functionsL, M, N, and P from the system (28)–(31) as follows:(32)L=S+hμ1+hμ+hβI,M=E+hβSI1+hV+hμ,P=I+hVE1+hα+hμ,Q=R+hαI1+hμ.
If the all eigenvalues of the Jacobean in formula (33),(33)Ј=∂L∂S∂L∂E∂L∂I∂L∂R∂M∂S∂M∂E∂M∂I∂M∂R∂P∂S∂Q∂S∂P∂E∂Q∂E∂P∂I∂P∂R∂Q∂I∂Q∂R,which satisfies the condition in (34), at the equilibria of the model(34)ρJ<1.
The elements of the Jacobean matrix are as follows:(35)∂L∂S=11+hμ+hβI,∂L∂E=0,∂L∂I=−hβS+hμ1+hμ+hβI2,∂L∂R=0,∂M∂S=hβI1+hV+hμ,∂M∂E=11+hV+hμ,∂M∂I=hβS1+hV+hμ,∂M∂R=0,∂P∂S=0,∂P∂E=hV1+hα+hμ,∂P∂I=11+hα+hμ,∂P∂R=0,∂Q∂S=0,∂Q∂E=0,∂Q∂I=hα1+hμ,∂Q∂R=11+hμ.
The Jacobean cough-free equilibrium (CFE-C1=S1,E1,I1,R1=1,0,0,0.(36)J1,0,0,0=11+hμ0−hβ1+hμ0011+hV+hμhβ1+hV+hμ00hV1+hα+hμ11+hα+hμ000hα1+hμ11+hμ,
The eigenvalues of the Jacobean matrix are as follows.
λ1=1/1+hμ<1,λ2=1/1+hμ≪1 provided that R0<1(37)J−λI=11+hV+hμhβ1+hV+hμhV1+hα+hμ11+hα+hμ,
A= Trace of J, B= det of J(38)A=2+hV+α+2μ1+hV+hμ1+hα+hμ,B=1−h2βV1+hV+hμ1+hα+hμ.
The Jacobean matrix at the cough existing equilibriumC2=S∗,E∗,I∗,R∗ is as follows:(39)JS∗,E∗,I∗,R∗=11+hμ+hβI∗0−hβS∗+hμ1+hμ+hNβI∗20hβI∗1+hV+hμ11+hV+hμhβS∗1+hV+hμ00hV1+hα+hμ11+hα+hμ000hα1+hμ11+hμ.
For equilibriaC2, Figure 2 depicts the behavior of the largest eigenvalues and the values of parameters are presented in Table 1, as desired.Figure 2
Largest Eigenvalue of the system at the equilibriaC2.


### 2.6.2. Positivity

Positivity is one of the most important physical traits of the state variables involved in the compartmental epidemic models. This property is investigated and ensured by the principle of mathematical induction for an implicit numerical integration scheme.Theorem 6.
Let the state variablesSt,Et,It, and Rt involved in the scheme are positive at t=0; furthermore, if all the parameters are also positive, then Sn+1≥0,En+1≥0,In+1≥0 and Rn+1≥0Proof.
Consider the equations (28)–(31) and apply the principle of mathematical induction, we proceed as follows:(40)Sn+1=Sn+hμ1+hμ+hNβInEn+1=En+hβSnIn1+hV+hμIn+1=In+hVEn1+hα+hμRn+1=Rn+hαIn1+hμ
At first, we putn=0 in the abovementioned system of equations and obtain the following results:(41)S1=S0+hμ1+hμ+hNβI0≥0.
Similarly,(42)E1=E0+hβS0I01+hV+hμ≥0,I1=I0+hVE01+hα+hμ≥0,R1=R0+hαI01+hμ≥0.
Now, forn=1, we arrive at the following expressions:(43)S2=S1+hμ1+hμ+hNβI1≥0.
Similarly,(44)E2=E1+hβS1I11+hV+hμ≥0,I2=I1+hVE11+hα+hμ≥0,R2=R1+hαI11+hμ≥0.
Next, suppose that the given set of equations guarantees the positivity property for the value of.n=2,3,4,…,n−1,i.e.,Sn≥0;En≥0;In≥0; and Rn≥0 for n=2,3,4,…,n−1 and state variables.
Now, the positivity will be investigated for an arbitrary positivity integern.
Forn∈Z, we observe that.(45)Sn+1=Sn+hμ1+hμ+hNβIn≥0,En+1=En+hβSnIn1+hV+hμ≥0,In+1=In+hVEn1+hα+hμ≥0,Rn+1=Rn+hαIn1+hμ≥0.
Hence, the projected scheme ensures the positivity for the state variablesSt,Et,It, and Rt for all the positive integer values of n.


### 2.6.3. Boundedness

Theorem 7.
LetS0,E0,I0, and R0 are finite, such that S0+E0+I0+R0≤1. Also μ,β,V,α are all positive in the model. Then, the discretized state variables, Sn+1,En+1,In+1, and Rn+1 are bounded by recrossing defined real constant bn+1 such that Sn+1,En+1,In+1 and Rn+1<bn+1∀n∈Z+where bn+1=4bn+hμ+hβbn2+hVbn+hαbn and b1=4+hμ+hβS0I0+hVE0+hαI0.Proof.
Considering the equations of an implicit numerical integration scheme for the state variablesS,E,I, and R as,(46)Sn+11+hμ+hNβIn=Sn+hμEn+11+hV+hμ=En+hβSnInIn+11+hα+hμ=In+hVEnRn+11+hμ=Rn+hαIn
We obtain the following expression by adding all the equations in the system (46).(47)Sn+1+En+1+In+1+Rn+1+Sn+1hμ+hNβIn+En+1hV+hμ+In+1hα+hμ+Rn+1hμ=Sn+En+In+Rn+hμ+hβSnIn+hVEn+hαIn.
We prove the boundedness of the nonstandard numerical scheme by using the theory of mathematical induction and the constraints imposed on the parameters and state variables.
By substitutingn=0, we obtain as follows:(48)S1+E1+I1+R1+S1hμ+hNβI0+E1hV+hμ+I1hα+hμ+R1hμ=S0+E0+I0+R0+hμ+hβS0I0+hVE0+hαI0,S1+E1+I1+R1+S1hμ+hNβI0+E1hV+hμ+I1hα+hμ+R1hμ<4+hμ+hβS0I0+hVE0+hαI0=b1,⇒S1+E1+I1+R1+S1hμ+hNβI0+E1hV+hμ+I1hα+hμ+R1hμ≤b1,⇒S11+hμ+hNβI0+E11+hV+hμ+I11+hα+hμ+R11+hμ≤b1,⇒S1<b1,E1<b1,I1<b1,R1<b1.
Now, by substitutingn=1 in (47) to reach the following step(49)S21+hμ+hNβI1+E21+hV+hμ+I21+hα+hμ+R21+hμ=S1+E1+I1+R1+hμ+hβS1I1+hVE1+hαI1,⇒S21+hμ+hNβI1+E21+hV+hμ+I21+hα+hμ+R21+hμ<4b1+hμ+hβb1+hVb1+hαb1=b2,⇒S2<b2,E2<b2,I2<b2,R2<b2.
Next, consider that the expression (47) is bounded for positive integers. 2,3,4,⋅⋅⋅,n−1i.e., Sn<bn,En<bn,In<bn and Rn<bn∀n∈2,3,4,…,n−1
Now, we investigate the boundedness of a positivity integern.(50)Sn+11+hμ+hNβIn+En+11+hV+hμ+In+11+hα+hμ+Rn+11+hμ=S1+E1+I1+R1+hμ+hβSnIn+hVEn+hαIn,Sn+11+hμ+hNβIn+En+11+hV+hμ+In+11+hα+hμ+Rn+11+hμ<4bn+hμ+hβbn2+hVbn+hαbn=bn+1.
So,Sn+1<bn+1,En+1<bn+1,In+1<bn+1 and. Rn+1<bn+1,where. bn+1=4bn+hμ+hβbn2+hVbn+hαbnb1=4+hμ+hβS0I0+hVE0+hαI0
Hence,Sn+1,En+1,In+1, and Rn+1 are bounded by a real number bn+1 for all positive values of n.


### 2.6.4. Consistency Analysis

In this section, the consistency analysis of an implicit numerical integration scheme is performed by using Taylor’s series expansion. First, we take equation (28) of the implicit numerical integration scheme and apply Taylor’s series expansion of Sn+1(51)Sn+1=Sn+hdSdt+h22!d2Sdt2+h33!d3Sdt3+⋅⋅⋅In the following expression:(52)Sn+11+hμ+hNβIn=Sn+hμ,Sn+hdSdt+h22!d2Sdt2+h33!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=Sn+hμ,Sn+hμSn+hNβInSn+hdSdt+hμhdSdt+hNβInhdSdt+h22!d2Sdt2+h33!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=Sn+hμ,hμSn+NβInSn+dSdt+hμdSdt+NβInhdSdt+h2!d2Sdt2+h23!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=hμ.We obtain as follows:(53)μSn+NβInSn+dSdt+hμdSdt+NβInhdSdt+h12!d2Sdt2+h3!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=μBy applyingh⟶0, we obtain as follows:(54)μS+NβIS+dSdt=μ,⇒dSdt=μ−μS−NβIS.This result implies that our discretized equation is consistent with equation (1) of the ODE model. Similarly, we take equation (29) and apply Taylor’s series expansion of En+1.(55)En+1=En+hdEdt+h22!d2Edt2+h33!d3Edt3+⋅⋅⋅In the following expression:(56)En+11+hV+hμ=En+hβSnIn.We obtain as follows:(57)VEn+μEn+dEdt+hVdEdt+hdEdtμ+h12!d2Edt2+h3!d3Edt3+⋅⋅⋅1+hV+hμ=βSnIn.Applyingh⟶0, we obtain as follows:(58)dEdt=βSI−VE−μE.Similarly, taking equation (30) and applying Taylor’s series expansion of In+1. After proper simplifications, we obtain as follows:(59)Inα+μIn+dIdt+hαdIdt+hμdIdtμ+h12!d2Idt2+h3!d3Idt3+⋅⋅⋅1+hα+hμ=VEn.Applyingh⟶0, we obtain as follows:(60)dIdt=VE−αI−μI.Similarly, taking equation (31) and applying Taylor’s series expansion of Rn+1, we obtain as follows:(61)μRn+dRdt+hμdRdt+h12!d2Rdt2+h3!d3Rdt3+⋅⋅⋅1+hμ=αIn.Applyingh⟶0, we obtain as follows:(62)dRdt=αI−μR.Hence, our discretized implicit numerical integration scheme is consistent with the ODES system (1–4).


## 2.6.1. Convergence Analysis

Theorem 5.
The proposed implicit numerical integration scheme is stable if it satisfies the assumption that the system’s Eigenvalues lie in the unit circle [27, 28].Proof.
Considering the functionsL, M, N, and P from the system (28)–(31) as follows:(32)L=S+hμ1+hμ+hβI,M=E+hβSI1+hV+hμ,P=I+hVE1+hα+hμ,Q=R+hαI1+hμ.
If the all eigenvalues of the Jacobean in formula (33),(33)Ј=∂L∂S∂L∂E∂L∂I∂L∂R∂M∂S∂M∂E∂M∂I∂M∂R∂P∂S∂Q∂S∂P∂E∂Q∂E∂P∂I∂P∂R∂Q∂I∂Q∂R,which satisfies the condition in (34), at the equilibria of the model(34)ρJ<1.
The elements of the Jacobean matrix are as follows:(35)∂L∂S=11+hμ+hβI,∂L∂E=0,∂L∂I=−hβS+hμ1+hμ+hβI2,∂L∂R=0,∂M∂S=hβI1+hV+hμ,∂M∂E=11+hV+hμ,∂M∂I=hβS1+hV+hμ,∂M∂R=0,∂P∂S=0,∂P∂E=hV1+hα+hμ,∂P∂I=11+hα+hμ,∂P∂R=0,∂Q∂S=0,∂Q∂E=0,∂Q∂I=hα1+hμ,∂Q∂R=11+hμ.
The Jacobean cough-free equilibrium (CFE-C1=S1,E1,I1,R1=1,0,0,0.(36)J1,0,0,0=11+hμ0−hβ1+hμ0011+hV+hμhβ1+hV+hμ00hV1+hα+hμ11+hα+hμ000hα1+hμ11+hμ,
The eigenvalues of the Jacobean matrix are as follows.
λ1=1/1+hμ<1,λ2=1/1+hμ≪1 provided that R0<1(37)J−λI=11+hV+hμhβ1+hV+hμhV1+hα+hμ11+hα+hμ,
A= Trace of J, B= det of J(38)A=2+hV+α+2μ1+hV+hμ1+hα+hμ,B=1−h2βV1+hV+hμ1+hα+hμ.
The Jacobean matrix at the cough existing equilibriumC2=S∗,E∗,I∗,R∗ is as follows:(39)JS∗,E∗,I∗,R∗=11+hμ+hβI∗0−hβS∗+hμ1+hμ+hNβI∗20hβI∗1+hV+hμ11+hV+hμhβS∗1+hV+hμ00hV1+hα+hμ11+hα+hμ000hα1+hμ11+hμ.
For equilibriaC2, Figure 2 depicts the behavior of the largest eigenvalues and the values of parameters are presented in Table 1, as desired.Figure 2
Largest Eigenvalue of the system at the equilibriaC2.


## 2.6.2. Positivity

Positivity is one of the most important physical traits of the state variables involved in the compartmental epidemic models. This property is investigated and ensured by the principle of mathematical induction for an implicit numerical integration scheme.Theorem 6.
Let the state variablesSt,Et,It, and Rt involved in the scheme are positive at t=0; furthermore, if all the parameters are also positive, then Sn+1≥0,En+1≥0,In+1≥0 and Rn+1≥0Proof.
Consider the equations (28)–(31) and apply the principle of mathematical induction, we proceed as follows:(40)Sn+1=Sn+hμ1+hμ+hNβInEn+1=En+hβSnIn1+hV+hμIn+1=In+hVEn1+hα+hμRn+1=Rn+hαIn1+hμ
At first, we putn=0 in the abovementioned system of equations and obtain the following results:(41)S1=S0+hμ1+hμ+hNβI0≥0.
Similarly,(42)E1=E0+hβS0I01+hV+hμ≥0,I1=I0+hVE01+hα+hμ≥0,R1=R0+hαI01+hμ≥0.
Now, forn=1, we arrive at the following expressions:(43)S2=S1+hμ1+hμ+hNβI1≥0.
Similarly,(44)E2=E1+hβS1I11+hV+hμ≥0,I2=I1+hVE11+hα+hμ≥0,R2=R1+hαI11+hμ≥0.
Next, suppose that the given set of equations guarantees the positivity property for the value of.n=2,3,4,…,n−1,i.e.,Sn≥0;En≥0;In≥0; and Rn≥0 for n=2,3,4,…,n−1 and state variables.
Now, the positivity will be investigated for an arbitrary positivity integern.
Forn∈Z, we observe that.(45)Sn+1=Sn+hμ1+hμ+hNβIn≥0,En+1=En+hβSnIn1+hV+hμ≥0,In+1=In+hVEn1+hα+hμ≥0,Rn+1=Rn+hαIn1+hμ≥0.
Hence, the projected scheme ensures the positivity for the state variablesSt,Et,It, and Rt for all the positive integer values of n.


## 2.6.3. Boundedness

Theorem 7.
LetS0,E0,I0, and R0 are finite, such that S0+E0+I0+R0≤1. Also μ,β,V,α are all positive in the model. Then, the discretized state variables, Sn+1,En+1,In+1, and Rn+1 are bounded by recrossing defined real constant bn+1 such that Sn+1,En+1,In+1 and Rn+1<bn+1∀n∈Z+where bn+1=4bn+hμ+hβbn2+hVbn+hαbn and b1=4+hμ+hβS0I0+hVE0+hαI0.Proof.
Considering the equations of an implicit numerical integration scheme for the state variablesS,E,I, and R as,(46)Sn+11+hμ+hNβIn=Sn+hμEn+11+hV+hμ=En+hβSnInIn+11+hα+hμ=In+hVEnRn+11+hμ=Rn+hαIn
We obtain the following expression by adding all the equations in the system (46).(47)Sn+1+En+1+In+1+Rn+1+Sn+1hμ+hNβIn+En+1hV+hμ+In+1hα+hμ+Rn+1hμ=Sn+En+In+Rn+hμ+hβSnIn+hVEn+hαIn.
We prove the boundedness of the nonstandard numerical scheme by using the theory of mathematical induction and the constraints imposed on the parameters and state variables.
By substitutingn=0, we obtain as follows:(48)S1+E1+I1+R1+S1hμ+hNβI0+E1hV+hμ+I1hα+hμ+R1hμ=S0+E0+I0+R0+hμ+hβS0I0+hVE0+hαI0,S1+E1+I1+R1+S1hμ+hNβI0+E1hV+hμ+I1hα+hμ+R1hμ<4+hμ+hβS0I0+hVE0+hαI0=b1,⇒S1+E1+I1+R1+S1hμ+hNβI0+E1hV+hμ+I1hα+hμ+R1hμ≤b1,⇒S11+hμ+hNβI0+E11+hV+hμ+I11+hα+hμ+R11+hμ≤b1,⇒S1<b1,E1<b1,I1<b1,R1<b1.
Now, by substitutingn=1 in (47) to reach the following step(49)S21+hμ+hNβI1+E21+hV+hμ+I21+hα+hμ+R21+hμ=S1+E1+I1+R1+hμ+hβS1I1+hVE1+hαI1,⇒S21+hμ+hNβI1+E21+hV+hμ+I21+hα+hμ+R21+hμ<4b1+hμ+hβb1+hVb1+hαb1=b2,⇒S2<b2,E2<b2,I2<b2,R2<b2.
Next, consider that the expression (47) is bounded for positive integers. 2,3,4,⋅⋅⋅,n−1i.e., Sn<bn,En<bn,In<bn and Rn<bn∀n∈2,3,4,…,n−1
Now, we investigate the boundedness of a positivity integern.(50)Sn+11+hμ+hNβIn+En+11+hV+hμ+In+11+hα+hμ+Rn+11+hμ=S1+E1+I1+R1+hμ+hβSnIn+hVEn+hαIn,Sn+11+hμ+hNβIn+En+11+hV+hμ+In+11+hα+hμ+Rn+11+hμ<4bn+hμ+hβbn2+hVbn+hαbn=bn+1.
So,Sn+1<bn+1,En+1<bn+1,In+1<bn+1 and. Rn+1<bn+1,where. bn+1=4bn+hμ+hβbn2+hVbn+hαbnb1=4+hμ+hβS0I0+hVE0+hαI0
Hence,Sn+1,En+1,In+1, and Rn+1 are bounded by a real number bn+1 for all positive values of n.


## 2.6.4. Consistency Analysis

In this section, the consistency analysis of an implicit numerical integration scheme is performed by using Taylor’s series expansion. First, we take equation (28) of the implicit numerical integration scheme and apply Taylor’s series expansion of Sn+1(51)Sn+1=Sn+hdSdt+h22!d2Sdt2+h33!d3Sdt3+⋅⋅⋅In the following expression:(52)Sn+11+hμ+hNβIn=Sn+hμ,Sn+hdSdt+h22!d2Sdt2+h33!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=Sn+hμ,Sn+hμSn+hNβInSn+hdSdt+hμhdSdt+hNβInhdSdt+h22!d2Sdt2+h33!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=Sn+hμ,hμSn+NβInSn+dSdt+hμdSdt+NβInhdSdt+h2!d2Sdt2+h23!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=hμ.We obtain as follows:(53)μSn+NβInSn+dSdt+hμdSdt+NβInhdSdt+h12!d2Sdt2+h3!d3Sdt3+⋅⋅⋅1+hμ+hNβIn=μBy applyingh⟶0, we obtain as follows:(54)μS+NβIS+dSdt=μ,⇒dSdt=μ−μS−NβIS.This result implies that our discretized equation is consistent with equation (1) of the ODE model. Similarly, we take equation (29) and apply Taylor’s series expansion of En+1.(55)En+1=En+hdEdt+h22!d2Edt2+h33!d3Edt3+⋅⋅⋅In the following expression:(56)En+11+hV+hμ=En+hβSnIn.We obtain as follows:(57)VEn+μEn+dEdt+hVdEdt+hdEdtμ+h12!d2Edt2+h3!d3Edt3+⋅⋅⋅1+hV+hμ=βSnIn.Applyingh⟶0, we obtain as follows:(58)dEdt=βSI−VE−μE.Similarly, taking equation (30) and applying Taylor’s series expansion of In+1. After proper simplifications, we obtain as follows:(59)Inα+μIn+dIdt+hαdIdt+hμdIdtμ+h12!d2Idt2+h3!d3Idt3+⋅⋅⋅1+hα+hμ=VEn.Applyingh⟶0, we obtain as follows:(60)dIdt=VE−αI−μI.Similarly, taking equation (31) and applying Taylor’s series expansion of Rn+1, we obtain as follows:(61)μRn+dRdt+hμdRdt+h12!d2Rdt2+h3!d3Rdt3+⋅⋅⋅1+hμ=αIn.Applyingh⟶0, we obtain as follows:(62)dRdt=αI−μR.Hence, our discretized implicit numerical integration scheme is consistent with the ODES system (1–4).


## 3. Conclusion

The implicit numerical integration scheme is designed for the consistent dynamical approximation of the SEIR whooping cough epidemic model. The results proved that the particular model restores the dynamical properties of a continuous model such as positivity, boundedness, consistency, and stability. Also, the implicit numerical integration scheme predicts the long-term perspectives for the dynamics of whooping cough disease. The implicit numerical integration scheme is unconditionally stable and consistent with the same behavior as the continuous model.



---

*Source: 3642444-2022-09-29.xml*