# Global Stability for Novel Complicated SIR Epidemic Models with the Nonlinear Recovery Rate and Transfer from Being Infectious to Being Susceptible to Analyze the Transmission of COVID-19

Authors: Fehaid Salem Alshammari; F. Talay Akyildiz
Journal: Journal of Function Spaces
 (2021)

---

## Abstract

Epidemiological models play pivotal roles in predicting, anticipating, understanding, and controlling present and future epidemics. The dynamics of infectious diseases is complex, and therefore, researchers need to consider more complicated mathematical models. In this paper, we first describe the dynamics of a complex SIR epidemic model with nonstandard nonlinear incidence and recovery rates. In this model, we consider the rate at which individuals lose immunity. Rigorous mathematical results have been established from the point of view of stability and bifurcation. The basic reproduction number (R0) is determined. We then apply LaSalle’s invariance principle and Lyapunov’s direct method to prove that the disease-free equilibrium is globally asymptotically stable when R0<1. The model has a unique endemic equilibrium when R0>1. A nonlinear Lyapunov function is used together with LaSalle’s invariance principle to show that the endemic equilibrium is globally asymptotically stable under some conditions. Further, for the case when R0=1, we analyze the model and show a backward bifurcation under certain conditions. In the second part of this paper, we analyze a modified SIR model with a vaccination term, which must be a function of time. We show that the modified model agrees well with COVID-19 data in Saudi Arabia. We then investigate different future scenarios. Simulation results suggest that a two-pronged strategy is crucial to control the COVID-19 pandemic in Saudi Arabia.


---

## Body

## 1. Introduction

There is a long and rich history of mathematical modeling of epidemiology. Most often, compartmental deterministic models are used for modeling the spread of infectious diseases [1–3]. In these models, a population of susceptible individuals evolves into other categories representing different stages of infection. Among many epidemiological models which had been used for infectious disease, SIR types of models have received more attention. In 1927, Kermack and McKendrick were the first to develop the susceptible-infective-recovered (SIR) model, where the total population is divided into three classes: susceptible, infective, and recovered [4]. After that, variant of SIR compartmental models were developed, some of them outlined in [5–8]. Here, we consider a complex SIR epidemic model with nonstandard nonlinear incidence and recovery rates. We also consider in our model the rate of losing immunity, which has not been considered before.In classic SIR epidemic models, the bilinear incident rateβIS/N (where Nt=St+It+Rt is the number of total populations, parameter β is the infection transmission rate, and St,It,andRt represent the number of susceptible and infected and recovered individuals at time t. Also, a linear recovery rate μI (μ is the per capita recovery rate) is often used. These classic models do not have bistability and periodicity in their solution which is not realistic, especially for COVID-19. Their dynamics basically depend on the basic reproducing number R0; the disease will be eliminated if R0<1; otherwise, the disease persists [9]. However, in reality, many infectious diseases show multiple peaks and/or periodic oscillations during the outbreak. The monotone incident rate term fIS (which describes the mechanism of disease transmission) does not capture the “psychological” or behavioral change and crowding effect of infected individuals. Therefore, we use the following general incidence rate:
(1)fIS=βIqS1+kIp,where k>0 measures the psychological or inhibitory effect and q and p are constants. This type of incidence rate was first considered by Liu et al. [10]. This incidence rate was used by many other scholars [11, 12]. In the rest of the paper, we use the incidence rate term (1) with q=1.Determining the treatment rate is not an easy task and many factors are involved in this process. The main factor is the number of health workforce that includes physicians, nurses, pharmacists, and other health care workers (HCW). The facilities of the hospital such as medical equipment and apparatus, the availability of the intensive care unit, and the number of the hospital beds and medicines are the other significant factors which are necessary and essential for safe and effective avoiding, diagnosis, and treatment of illness [13, 14]. In this work, we follow the work of Shan and Zhu [15] and use the following nonlinear treatment function
(2)μb,I=μ0+μ1−μ0bb+I,where μ0,μ1μ1>μ0 are the minimum and maximum per capita recovery rates, respectively. Parameter b is considered as a measure of available hospital resources.We need herd immunity to eradicate the COVID-19 pandemic from the human population. This could be achieved either by previous infection or by vaccination. Several pharmacological companies declared high efficacy rates of their vaccine products [16, 17]. In Saudi Arabia, the Ministry of Health launched a vaccine campaign through a mobile application, which is called “Sehaty,” that provides registration for COVID-19 vaccination, and vaccination centers were established in different cities around the country. The campaign was launched offering both the Pfizer-BioNTech and AstraZeneca’s COVID-19 vaccines. They aim to provide free vaccination to all citizens and residents until getting herd immunity [18]. Hence, it is extremely important to create public awareness on the importance of vaccination. In this paper, we develop a mathematical model to show the effectiveness of vaccination in reducing the infection rates in Saudi Arabia.The organization of this paper is as follows: the model framework is given in Section2, the existence of equilibria and global stability of disease-free and endemic equilibria is studied in Section 3. In Section 4, we study the backward bifurcation. Section 5 is devoted to the modification of the model that was previously defined in Section 2. In Section 6, we summarize our results and provide a short discussion on possible extensions of our model with a possibility of additional vaccinated component.


## 2. Model Framework

In this section, we describe the mathematical formulation of an SIR epidemic model, where the total population size of individuals is represented byNt. The total population size is subdivided into three different classes, namely, St: susceptible, It: infected, and Rt: removed or recovered individuals. We consider that the (viruses) parasites of the diseases are transmitted to the susceptible populations by the direct contact with the infected populations.We assume that the total recruitment at any timet is "A" and all the new recruited populations go to the susceptible class. We consider that β≥0, which is the disease transmission rate and the incidence rate to be βIS/1+kI, where k is the half-saturation constant for which the susceptible population decreases. The population of susceptible people is decreased by the natural death rate d, and the recovery population goes to the susceptible class at a rate δ (the recovered individuals could become susceptible again due to lose of immunity). Hence, the governing equation can be modeled as
(3)dSdt=A−βIS1+kI−dS+δR.The infected population is decreased by the natural death rated and the disease death rate γ. The medical treatments, determining how well the diseases are controlled, are normally expressed as constant recovery rates. There are several suggested functions to describe the treatment term. In this paper, we follow the work of Shan and Zhu [15], where they defined the recovery rate as a function of b, the number of hospital beds, and the number of infectives I. Thus, the time rate of change for this can be represented by the following equation:
(4)dIdt=βIS1+kI−μ0+μ1−μ0bb+II−d+γI.The recovered population is increased by the recovery rate as in above and decreased by the natural death rated and the susceptibility of recovered individuals δ. The time rate of change for the population of recovered individuals can be represented by the following equation:
(5)dRdt=μ0+μ1−μ0bb+II−d+δR.All parameters involved in equations (3)–(5) are nonnegative. Equations (3)–(5) can be written in the combined form as follows:
(6)dSdt=A−βIS1+kI−dS+δR,(7)dIdt=βIS1+kI−μ0+μ1−μ0bb+II−d+γI,(8)dRdt=μ0+μ1−μ0bb+II−d+δR,subject to initial conditions
(9)S0=S0≥0,(10)I0=I0≥0,(11)R0=R0≥0.For systems (6)–(9), the cone R3+ is positively invariant. The C1 smoothness of the right side of systems (6)–(9) implies local existence and uniqueness of solutions with the initial values in R3+.It is not difficult to show that every solution of (6)–(9) with nonnegative initial conditions remains nonnegative. If we add up the three equations of systems (6)–(8), we get
(12)dS+I+Rdt=A−dS+I+R−γI≤A−dS+I+R,which implies that the set
(13)Ω=S,I,R∈R3+:S≥0,I≥,R≥0,S+I+R≤Adis positively invariant and an attractive set for (6)–(9); hence, all solutions in the first octant approach enter or stay inside the set defined above and also bounded and, hence, globally exist. Thus, the initial value problem of systems (6)–(9) is mathematically well posed and epidemiologically reasonable.Because of total populationNt=St+It+Rt is regulated by the disease, the system cannot be reduced to the lower dimension; hence, we need to analyze systems (6)–(9) in the three-dimensional phase space.


## 3. Basic Reproduction Number

The basic reproduction number is denoted byR0, and it is defined as the number of newly infected individuals caused by a single infection. We now find the basic reproduction number of systems (6)–(9) using the next-generation matrix method developed by Driessche and Watmough [19]. Let us write systems (6)–(9) as follows:
(14)x.i=Fix−Vix,i=1,⋯,3,where
(15)xi=I,S,RT,Fi=βS1+kI00,Vi=μ0+μ1−μ0bb+II+d+γI−A+dS−δR+βS1+kI−μ0+μ1−μ0bb+II+d+δR.We denote the disease-free equilibrium of models (6)–(9) by E0, where
(16)E0∗=0,Ad,0.Now, the Jacobian matrixes ofFi and Vi at E0 are given as
(17)F=JFi=βAd000,JVi=V0M1M2,where
(18)V=d+γ+μ10βAdd.The basic reproduction numberR0 of system (14) is defined by the spectral radius of the matrix ρFV−1 (see Driessche and Watmough [19]) and it is given by R0=βA/dd+γ+μ1.


## 4. Existence and Types of Equilibria

### 4.1. Existence of Equilibria

We have already established the existence of the disease-free equilibriumE0=A/d,0,0 for systems (6)–(9) for all values of parameters. For any endemic equilibrium E∗=S,I,R, its coordinates satisfy
(19)R=1d+δμ0+μ1−μ0bb+II.Substituting this into the first equation in (6)–(9) and solving for S, we obtain
(20)S=A+δ/d+δμ0+μ1−μ0b/b+IβIS/1+kI+d.The variableI should be the positive root of the following quadratic equation
(21)fI=AI2+CI+D,where
(22)A=d3k+γ+β+δk+kμ0d2+γβ+γδk+βδ+βμ0+δkμ0d+γβδdγ+d+μ1,C=bk+1d3+γbk+bβ+bkδ+bkμ1+δ+μ0+γd2+γbδβ+kbδμ+γbβ+bβδ+bβδ+bβμ1+kμ0+γδ−Aβd+γbβδ−δβAdγ+d+μ1,and(23)D=bd+δdγ+d+μ11−R0.EquationfI=0 may have two roots if Δ0>0, which are
(24)I1=−C−Δ02A,I1=−C+Δ02A,where Δ0=C2−4AD.IfI>0, we see from (19) and (20) that S>0 and R>0, respectively.We study the existence of equilibria in the following three cases:(25)R0>1.In this case,D>0; since A>0, we have I1<0 and I2>0, so systems (6)–(9) have a unique endemic equilibrium E2=SI2,I2,RI2.
(26)R0=1.In this case,I1=−C/A and I2=0; if C<0, then, systems (6)–(9) have a unique endemic equilibrium E1=SI1,I1,RI1.
(27)R0<1.In this case, ifC>0, there is no endemic equilibrium; if C<0 and Δ0>0, we have two endemic equilibria E1=SI1,I1,RI1 and E2=SI2,I2,RI2. If Δ0=0, then, we have one root multiplicity 2. We summarize the result in the following theorem.Theorem 1.
For systems (6)–(9),
(1)
The disease-free equilibrium exists(2)
IfR0>1, there exists a unique endemic equilibrium E2(3)
IfR0=1, there exists a unique endemic equilibrium E1 provided that C<0; otherwise, there is no endemic equilibrium(4)
IfR0<1 and if C>0, there is no endemic equilibrium; if C<0 and Δ0>0, the system has two endemic equilibria E1 and E2; if C<0 and Δ0=0, the system has two equilibria coalesce into E∗Theorem 2.
For systems (6)–(9), the disease-free equilibrium E0=A/d,0,0 is
(i)
R0<1: an attracting node(ii)
R0>1:a hyperbolic saddle(iii)
R0=1: and b>d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d: a saddle-node of codimension 1
b<d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d is a saddle-node of codimension 1; b=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d is an attracting semihyperbolic node of codimension 2.Proof.
The Jacobian matrixJE0 for systems (6)–(9) is given by
(28)JE0=−d−βAdδ0βAd−μ1−γ00μ1−d−δ.
The eigenvalues are obtained from this.(29)λ1=−d,λ2=−d−δ,λ3=Aβ−γd−d2−dμ1d.
SinceR0<1, λ3 is negative, because λ3=γd+d2+dμ1/dR0−1; hence, all eigenvalues are negative. Thus, if R0<1, E0 is an attracting node, and if R0>1, then, E0 is a hyperbolic saddle. If R0=1, the third eigenvalue is zero. In order to determine the type of E0, we first transform the disease-free equilibrium point E0 to the origin. We use S~=S−A/d; then, using the Taylor expansion as in [15], we get
(30)dS~dt=δRt−dS~t+βAItd+OS~,I,R3,(31)dIdt=−βAkd+−μ1+μ0bI2+OS~,I,R3,(32)dRdt=−d+δR+μ1It+OS~,I,R3,
Using eigenvectors, we find following transform:(33)It=Xt,S~t=−δμ1/d+δ+βA/dIt+Yt−δd/d+δ2d+δZtd,Rt=Zt−μ1Itd+δ.
Thus, using these variables in (30), we obtain
(34)dXdt=−βAkd+−μ1+μ0b+β−δμ1/d+δ+βA/ddX2+XOY,Z+OY,Z2,X,Y,Z3,(35)dYdt=−dY+OX,Y,Z2,(36)dZdt=−d+δZ+OX,Y,Z2.
It is unnecessary to calculate the center manifold ifb≠d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d, and E0 is a saddle node. If b=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d, then, from the central manifold theorem for systems (6)–(9), we obtain
(37)dXdt=−βAkd+β−δμ1/d+δ+βA/dd+δd−μ1+μ0b2d+δ+δμ1d+δ−βAkdX3+OX4,(38)dYdt=−dY+OX,Y,Z2,(39)dZdt=−d+δZ+OX,Y,Z2.
Hence,E0 is a semihyperbolic attracting node.☐From Theorem1, one can see that for the existence of endemic equilibrium, βA<dd+γ+μ1 is a necessary condition. If <dd+γ+μ1, E2 is the unique equilibrium point and we have the following theorem.Theorem 3.
For the system (4)-(7), if <dd+γ+μ1, the disease-free equilibrium E0=A/d,0,0 is globally asymptotically stable.Proof.
Consider the Lyapunov functionV=I in R3+ with the Liapunov derivative
(40)V′=βIS1+kI−μ0+μ1−μ0bb+II−d+γI≤βS−μ0+μ1−μ0bb+I−d+γI≤limS⟶Ad,I⟶0βS−μ0+μ1−μ0bb+I−d+γI≤βAd−μ1+d+γlimI⟶0I≤μ1+d+γR0−1limI⟶0I.
The Lyapunov–Lasalle theorem implies that solutions inΩ approach the largest positively invariant subset of set V′=0, i.e., plane I=0. In this plane, S⟶A/d and R⟶0 as t⟶∞. Thus, all solutions in plane I=0 go to the disease-free equilibrium E0. Therefore, E0 is globally asymptotically stable.☐Theorem 4.
A sufficient condition for the endemic equilibriumE2 to be locally asymptotically stable is R0>1Proof.
In this case, the Jacobian matrix has the following form:(41)J=−βI∗1+I∗k−dΣδβI∗1+I∗k−Σ+Ξ−d−γ00−Ξ−d−δ,where Σ=−βS∗/1+I∗k+βkI∗S∗/1+I∗k2 and Ξ=μ1−μ0b/b+I∗I∗/b+I∗−1−μ0. The characteristic equation for the above Jacobian around its endemic equilibrium E2 is
(42)x3+C1x2+C2x+C3=0,where C1,C2, and C3 are too long to reproduce here. We have observed that C1,C2, and C3 are positive and C1C2>C3; hence, the Routh–Hurwitz criterion is satisfied, so systems (6)–(9) are locally asymptotically stable for R0>1.☐Theorem 5.
The epidemic models (4)–(7) at E2 are globally asymptotically stable if R0>1.Proof.
We consider the Lyapunov function given by(43)L=∫S∗S1−S∗zdz+∫I∗I1−I∗zdz+∫R∗S1−R∗zdz.
Taking the derivative, we have(44)L.=1−S∗SdSdt+1−I∗IdIdt+1−R∗RdRdt=1−S∗SA−βIS1+kI−dS+δR+1−I∗IβIS1+kI−μ0+μ1−μ0bb+II−d+γI+1−R∗Rμ0+μ1−μ0bb+II−d+δR=1−S∗SβI∗S∗1+kI∗+dS∗+δR∗−βIS1+kI−dS−δR+1−I∗IβIS1+kI−μ0+μ1−μ0bb+II−d+γI−βI∗S∗1+kI∗+μ0+μ1−μ0bb+I∗I∗+d+γI∗+1−R∗Rμ0+μ1−μ0bb+II−d+δR−μ0+μ1−μ0bb+I∗I∗+d+δR∗,L.=1−S∗S11+kI∗1+kIβI∗SS∗S−II∗+kβI∗IS∗1−SS∗+δR∗1−RR∗+dS∗1−SS∗+1−I∗I−11+kI∗1+kIβI∗SS∗S−II∗+kβI∗IS∗1−SS∗+b2I∗b+I∗b+I1−II∗+μ0+d+γI∗1−II∗+1−R∗Rd+δR∗1−RR∗−b2I∗b+I∗b+I1−II∗−μ0+d+γI∗1−II∗,L.=1−S∗SkβI∗IS∗1−SS∗+1−I∗Ib2I∗b+I∗b+I1−II∗+μ0+d+γI∗1−II∗+1−R∗Rd+δR∗1−RR∗+1−R∗RδR∗1−S∗S−b2I∗b+I∗b+I1−II∗−μ0+d+γI∗1−II∗+1−I∗I−11+kI∗1+kIβI∗SS∗S−II∗+kβI∗IS∗1−SS∗+11+kI∗1+kIβI∗SS∗S−II∗1−S∗S,if
(45)1−R∗R1−S∗S<0,1−R∗R1−II∗>0,1−I∗I−11+kI∗1+kIβI∗SS∗S−II∗+kβI∗IS∗1−SS∗+11+kI∗1+kIβI∗SS∗S−II∗1−S∗S<0.
Then, we have(46)L.<0,
Hence, the theorem is proved.☐Theorem 6.
For systems (6)–(9), consider R0 as the bifurcation parameter. Then, we have, when R0=1, systems (6)–(9) which undergo forward bifurcation if b>d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d; systems (6)–(9) undergo backward bifurcation if b<d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d; systems (6)–(9) undergo pitchfork bifurcation if b=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d.Proof.
SinceR0 is a function of the parameters β,γ,δ,d,andμ1, without loss of generality, we can choose μ1 as the bifurcation parameter.
Letμ1=βA/d−d−γ+ε, we substitute this into (6)–(9); we note that if ε=0, then, it reduces to R0=1. We then use Taylor expansion at E0; diagonalizing the linear part, we then apply the center manifold theorem for the parameter ε. We found that
(47)dXdt=−ε+Oε2X−βAkd+−Π+μ0b+β−δΠ/d+δ+βA/dd+OεX2+OX3,where Π=βA/d−d−γ. Denoting the right-hand side of (47) as Γε,X, we have
(48)Γ0,0=0,∂∂εΓ0,0=0,∂∂XΓ0,0=0,∂2∂X∂εΓ0,0=−1,∂2∂X2Γ0,0=−βAkd+−Π+μ0b+β−δΠ/d+δ+βA/dd=−βAkd+−μ1+μ0b+β−δμ1/d+δ+βA/dd.
Therefore, the system has transcritical bifurcation if(49)b≠d2μ1−μ0d+μ1−μ0δβAd+kdk+β−δμ1d.
Ifb=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d, using the center manifold theory, as the above, from (6)–(9), we obtain
(50)dXdt=−ε+Oε2X+OεX2−βAkd+β−δΠ/d+δ+βA/dd+δd−Π+μ0b2d+δ+δΠd+δ−βAkd+OεX3+OX4.
Again, denoting the right-hand side of (50) as Γε,X, we have
(51)Γ0,0=0,∂∂εΓ0,0=0,∂∂XΓ0,0=0,∂2∂X∂εΓ0,0=−1,∂2∂X2Γ0,0=0,∂3∂X3Γ0,0=−βAkd+β−δΠ/d+δ+βA/dd+δd−Π+μ0b2d+δ+δΠd+δ−βAkd=−βAkd+β−δΠ/d+δ+βA/dd+δd−Π+μ0b2d+δ+δΠd+δ−βAkd.
It is clear that systems (6)–(9) have pitchfork bifurcation if b=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d when R0=1 [20]. We note that If ε=0, the system related with equation (47) and the system related with equation (50) reduce to the systems in equations (34) and (37), respectively.
SinceR0 is a function of μ1, also bandk, we choose μ1 as the bifurcation parameter. In Theorem 1, there are two endemic equilibria provided that if R0<1, C<0, and Δ0>0 and if C<0 and Δ0=0, then, there are two equilibria coalesce.
The basic reproducing numberR0=1 defines a straight line P0 in the μ1,b plane
(52)P0:μ1=Aβ−γd−d2d,
C=0 defines one branch of hyperbola PC (see Figure 1):
(53)PC:b=HCμ1=Aβ−γd−d2−dμ0d+kd3k+γ+δ+μ1k+βd2+δγ+μ1k+βγ+δ+μ1d+γβδ.
The branch ofPC has an intersection with P0 at point K, where K=Aβ−γd−d2/d,d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d. We also have dHCμ1/dμ1<0andd2HCμ1/dμ12>0, so HCμ1 is decreasing and a convex function of μ1.
Now, let(54)P0+=μ1,b:μ1=Aβ−γd−d2d,b>d2μ1−μ0d+μ1−μ0δβAd+kdk+β−δμ1d,P0−=μ1,b:μ1=Aβ−γd−d2d,b<d2μ1−μ0d+μ1−μ0δβAd+kdk+β−δμ1d.
Then,P0=P0+∪P0−∪K, where P0+andP0− are the parts of P0 which are separated by K. Now, we consider the curve defined by Δ0μ1,b=0 and we indicate these curves as PΔ±; solving for b from Δ0μ1,b=0, we got
(55)PΔ±:b=fΔ±μ1.
We found thatfΔ+Aβ−γd−d2/d=fΔ−Aβ−γd−d2/d and straightforward calculations lead to
(56)fΔ−μ1<HCμ1<fΔ+μ1,μ1∈Aβ−γd−d2d,∞.
Furthermore, we have found that(57)dfΔ±μ1dμ1<0,μ1∈Aβ−γd−d2d,∞andlimμ1⟶∞fΔ−μ1⟶0.d2fΔ±μ1dμ12>0,μ1∈Aβ−γd−d2d,∞andlimμ1⟶∞fΔ−μ1⟶0.
Hence,fΔ±μ1 and HCμ1 are decreasing and they are convex functions of μ1.
Based on the above discussion and Theorem1, for fix k, we define
(58)D0=b,μ1:b>fΔ−μ1,μ1>Aβ−γd−d2d,D1=b,μ1:b>0,μ0<μ1<Aβ−γd−d2d,D2=b,μ1:0<b<fΔ−μ1,μ1>Aβ−γd−d2d.
In Figure1, we see that there is one endemic equilibrium in region D1 and two equilibria in region D2. The system of equations (6)–(9) undergoes saddle-node bifurcation on the curve PΔ−. The forward bifurcation occurs on PΔ− and we have backward bifurcation, which occurs on PΔ−. The pitchfork bifurcation occurs when transversally passing through curve P0 at point K. The mathematical system of equations (6)–(9) has a semihyperbolic node of codimension 2 at point K. The similar discussion can be done for fixed b and varying value of k to show the backward bifurcation. This part is omitted since the analysis is the same as before.☐Figure 1
For fixedk, the bifurcation curves in the μ1,b plane when Aβ>dμ0+γd+d2. In D1 region R0>1, where there exists a unique endemic equilibrium. There are two endemic equilibria in D2. Two equilibria coalesce and a saddle-node bifurcation occurs on PΔ− . The forward bifurcation occurs on P0− and backward bifurcation occurs on P0+. There is no endemic equilibrium in the first region.


## 4.1. Existence of Equilibria

We have already established the existence of the disease-free equilibriumE0=A/d,0,0 for systems (6)–(9) for all values of parameters. For any endemic equilibrium E∗=S,I,R, its coordinates satisfy
(19)R=1d+δμ0+μ1−μ0bb+II.Substituting this into the first equation in (6)–(9) and solving for S, we obtain
(20)S=A+δ/d+δμ0+μ1−μ0b/b+IβIS/1+kI+d.The variableI should be the positive root of the following quadratic equation
(21)fI=AI2+CI+D,where
(22)A=d3k+γ+β+δk+kμ0d2+γβ+γδk+βδ+βμ0+δkμ0d+γβδdγ+d+μ1,C=bk+1d3+γbk+bβ+bkδ+bkμ1+δ+μ0+γd2+γbδβ+kbδμ+γbβ+bβδ+bβδ+bβμ1+kμ0+γδ−Aβd+γbβδ−δβAdγ+d+μ1,and(23)D=bd+δdγ+d+μ11−R0.EquationfI=0 may have two roots if Δ0>0, which are
(24)I1=−C−Δ02A,I1=−C+Δ02A,where Δ0=C2−4AD.IfI>0, we see from (19) and (20) that S>0 and R>0, respectively.We study the existence of equilibria in the following three cases:(25)R0>1.In this case,D>0; since A>0, we have I1<0 and I2>0, so systems (6)–(9) have a unique endemic equilibrium E2=SI2,I2,RI2.
(26)R0=1.In this case,I1=−C/A and I2=0; if C<0, then, systems (6)–(9) have a unique endemic equilibrium E1=SI1,I1,RI1.
(27)R0<1.In this case, ifC>0, there is no endemic equilibrium; if C<0 and Δ0>0, we have two endemic equilibria E1=SI1,I1,RI1 and E2=SI2,I2,RI2. If Δ0=0, then, we have one root multiplicity 2. We summarize the result in the following theorem.Theorem 1.
For systems (6)–(9),
(1)
The disease-free equilibrium exists(2)
IfR0>1, there exists a unique endemic equilibrium E2(3)
IfR0=1, there exists a unique endemic equilibrium E1 provided that C<0; otherwise, there is no endemic equilibrium(4)
IfR0<1 and if C>0, there is no endemic equilibrium; if C<0 and Δ0>0, the system has two endemic equilibria E1 and E2; if C<0 and Δ0=0, the system has two equilibria coalesce into E∗Theorem 2.
For systems (6)–(9), the disease-free equilibrium E0=A/d,0,0 is
(i)
R0<1: an attracting node(ii)
R0>1:a hyperbolic saddle(iii)
R0=1: and b>d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d: a saddle-node of codimension 1
b<d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d is a saddle-node of codimension 1; b=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d is an attracting semihyperbolic node of codimension 2.Proof.
The Jacobian matrixJE0 for systems (6)–(9) is given by
(28)JE0=−d−βAdδ0βAd−μ1−γ00μ1−d−δ.
The eigenvalues are obtained from this.(29)λ1=−d,λ2=−d−δ,λ3=Aβ−γd−d2−dμ1d.
SinceR0<1, λ3 is negative, because λ3=γd+d2+dμ1/dR0−1; hence, all eigenvalues are negative. Thus, if R0<1, E0 is an attracting node, and if R0>1, then, E0 is a hyperbolic saddle. If R0=1, the third eigenvalue is zero. In order to determine the type of E0, we first transform the disease-free equilibrium point E0 to the origin. We use S~=S−A/d; then, using the Taylor expansion as in [15], we get
(30)dS~dt=δRt−dS~t+βAItd+OS~,I,R3,(31)dIdt=−βAkd+−μ1+μ0bI2+OS~,I,R3,(32)dRdt=−d+δR+μ1It+OS~,I,R3,
Using eigenvectors, we find following transform:(33)It=Xt,S~t=−δμ1/d+δ+βA/dIt+Yt−δd/d+δ2d+δZtd,Rt=Zt−μ1Itd+δ.
Thus, using these variables in (30), we obtain
(34)dXdt=−βAkd+−μ1+μ0b+β−δμ1/d+δ+βA/ddX2+XOY,Z+OY,Z2,X,Y,Z3,(35)dYdt=−dY+OX,Y,Z2,(36)dZdt=−d+δZ+OX,Y,Z2.
It is unnecessary to calculate the center manifold ifb≠d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d, and E0 is a saddle node. If b=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d, then, from the central manifold theorem for systems (6)–(9), we obtain
(37)dXdt=−βAkd+β−δμ1/d+δ+βA/dd+δd−μ1+μ0b2d+δ+δμ1d+δ−βAkdX3+OX4,(38)dYdt=−dY+OX,Y,Z2,(39)dZdt=−d+δZ+OX,Y,Z2.
Hence,E0 is a semihyperbolic attracting node.☐From Theorem1, one can see that for the existence of endemic equilibrium, βA<dd+γ+μ1 is a necessary condition. If <dd+γ+μ1, E2 is the unique equilibrium point and we have the following theorem.Theorem 3.
For the system (4)-(7), if <dd+γ+μ1, the disease-free equilibrium E0=A/d,0,0 is globally asymptotically stable.Proof.
Consider the Lyapunov functionV=I in R3+ with the Liapunov derivative
(40)V′=βIS1+kI−μ0+μ1−μ0bb+II−d+γI≤βS−μ0+μ1−μ0bb+I−d+γI≤limS⟶Ad,I⟶0βS−μ0+μ1−μ0bb+I−d+γI≤βAd−μ1+d+γlimI⟶0I≤μ1+d+γR0−1limI⟶0I.
The Lyapunov–Lasalle theorem implies that solutions inΩ approach the largest positively invariant subset of set V′=0, i.e., plane I=0. In this plane, S⟶A/d and R⟶0 as t⟶∞. Thus, all solutions in plane I=0 go to the disease-free equilibrium E0. Therefore, E0 is globally asymptotically stable.☐Theorem 4.
A sufficient condition for the endemic equilibriumE2 to be locally asymptotically stable is R0>1Proof.
In this case, the Jacobian matrix has the following form:(41)J=−βI∗1+I∗k−dΣδβI∗1+I∗k−Σ+Ξ−d−γ00−Ξ−d−δ,where Σ=−βS∗/1+I∗k+βkI∗S∗/1+I∗k2 and Ξ=μ1−μ0b/b+I∗I∗/b+I∗−1−μ0. The characteristic equation for the above Jacobian around its endemic equilibrium E2 is
(42)x3+C1x2+C2x+C3=0,where C1,C2, and C3 are too long to reproduce here. We have observed that C1,C2, and C3 are positive and C1C2>C3; hence, the Routh–Hurwitz criterion is satisfied, so systems (6)–(9) are locally asymptotically stable for R0>1.☐Theorem 5.
The epidemic models (4)–(7) at E2 are globally asymptotically stable if R0>1.Proof.
We consider the Lyapunov function given by(43)L=∫S∗S1−S∗zdz+∫I∗I1−I∗zdz+∫R∗S1−R∗zdz.
Taking the derivative, we have(44)L.=1−S∗SdSdt+1−I∗IdIdt+1−R∗RdRdt=1−S∗SA−βIS1+kI−dS+δR+1−I∗IβIS1+kI−μ0+μ1−μ0bb+II−d+γI+1−R∗Rμ0+μ1−μ0bb+II−d+δR=1−S∗SβI∗S∗1+kI∗+dS∗+δR∗−βIS1+kI−dS−δR+1−I∗IβIS1+kI−μ0+μ1−μ0bb+II−d+γI−βI∗S∗1+kI∗+μ0+μ1−μ0bb+I∗I∗+d+γI∗+1−R∗Rμ0+μ1−μ0bb+II−d+δR−μ0+μ1−μ0bb+I∗I∗+d+δR∗,L.=1−S∗S11+kI∗1+kIβI∗SS∗S−II∗+kβI∗IS∗1−SS∗+δR∗1−RR∗+dS∗1−SS∗+1−I∗I−11+kI∗1+kIβI∗SS∗S−II∗+kβI∗IS∗1−SS∗+b2I∗b+I∗b+I1−II∗+μ0+d+γI∗1−II∗+1−R∗Rd+δR∗1−RR∗−b2I∗b+I∗b+I1−II∗−μ0+d+γI∗1−II∗,L.=1−S∗SkβI∗IS∗1−SS∗+1−I∗Ib2I∗b+I∗b+I1−II∗+μ0+d+γI∗1−II∗+1−R∗Rd+δR∗1−RR∗+1−R∗RδR∗1−S∗S−b2I∗b+I∗b+I1−II∗−μ0+d+γI∗1−II∗+1−I∗I−11+kI∗1+kIβI∗SS∗S−II∗+kβI∗IS∗1−SS∗+11+kI∗1+kIβI∗SS∗S−II∗1−S∗S,if
(45)1−R∗R1−S∗S<0,1−R∗R1−II∗>0,1−I∗I−11+kI∗1+kIβI∗SS∗S−II∗+kβI∗IS∗1−SS∗+11+kI∗1+kIβI∗SS∗S−II∗1−S∗S<0.
Then, we have(46)L.<0,
Hence, the theorem is proved.☐Theorem 6.
For systems (6)–(9), consider R0 as the bifurcation parameter. Then, we have, when R0=1, systems (6)–(9) which undergo forward bifurcation if b>d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d; systems (6)–(9) undergo backward bifurcation if b<d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d; systems (6)–(9) undergo pitchfork bifurcation if b=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d.Proof.
SinceR0 is a function of the parameters β,γ,δ,d,andμ1, without loss of generality, we can choose μ1 as the bifurcation parameter.
Letμ1=βA/d−d−γ+ε, we substitute this into (6)–(9); we note that if ε=0, then, it reduces to R0=1. We then use Taylor expansion at E0; diagonalizing the linear part, we then apply the center manifold theorem for the parameter ε. We found that
(47)dXdt=−ε+Oε2X−βAkd+−Π+μ0b+β−δΠ/d+δ+βA/dd+OεX2+OX3,where Π=βA/d−d−γ. Denoting the right-hand side of (47) as Γε,X, we have
(48)Γ0,0=0,∂∂εΓ0,0=0,∂∂XΓ0,0=0,∂2∂X∂εΓ0,0=−1,∂2∂X2Γ0,0=−βAkd+−Π+μ0b+β−δΠ/d+δ+βA/dd=−βAkd+−μ1+μ0b+β−δμ1/d+δ+βA/dd.
Therefore, the system has transcritical bifurcation if(49)b≠d2μ1−μ0d+μ1−μ0δβAd+kdk+β−δμ1d.
Ifb=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d, using the center manifold theory, as the above, from (6)–(9), we obtain
(50)dXdt=−ε+Oε2X+OεX2−βAkd+β−δΠ/d+δ+βA/dd+δd−Π+μ0b2d+δ+δΠd+δ−βAkd+OεX3+OX4.
Again, denoting the right-hand side of (50) as Γε,X, we have
(51)Γ0,0=0,∂∂εΓ0,0=0,∂∂XΓ0,0=0,∂2∂X∂εΓ0,0=−1,∂2∂X2Γ0,0=0,∂3∂X3Γ0,0=−βAkd+β−δΠ/d+δ+βA/dd+δd−Π+μ0b2d+δ+δΠd+δ−βAkd=−βAkd+β−δΠ/d+δ+βA/dd+δd−Π+μ0b2d+δ+δΠd+δ−βAkd.
It is clear that systems (6)–(9) have pitchfork bifurcation if b=d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d when R0=1 [20]. We note that If ε=0, the system related with equation (47) and the system related with equation (50) reduce to the systems in equations (34) and (37), respectively.
SinceR0 is a function of μ1, also bandk, we choose μ1 as the bifurcation parameter. In Theorem 1, there are two endemic equilibria provided that if R0<1, C<0, and Δ0>0 and if C<0 and Δ0=0, then, there are two equilibria coalesce.
The basic reproducing numberR0=1 defines a straight line P0 in the μ1,b plane
(52)P0:μ1=Aβ−γd−d2d,
C=0 defines one branch of hyperbola PC (see Figure 1):
(53)PC:b=HCμ1=Aβ−γd−d2−dμ0d+kd3k+γ+δ+μ1k+βd2+δγ+μ1k+βγ+δ+μ1d+γβδ.
The branch ofPC has an intersection with P0 at point K, where K=Aβ−γd−d2/d,d2μ1−μ0d+μ1−μ0δ/βAd+kdk+β−δμ1d. We also have dHCμ1/dμ1<0andd2HCμ1/dμ12>0, so HCμ1 is decreasing and a convex function of μ1.
Now, let(54)P0+=μ1,b:μ1=Aβ−γd−d2d,b>d2μ1−μ0d+μ1−μ0δβAd+kdk+β−δμ1d,P0−=μ1,b:μ1=Aβ−γd−d2d,b<d2μ1−μ0d+μ1−μ0δβAd+kdk+β−δμ1d.
Then,P0=P0+∪P0−∪K, where P0+andP0− are the parts of P0 which are separated by K. Now, we consider the curve defined by Δ0μ1,b=0 and we indicate these curves as PΔ±; solving for b from Δ0μ1,b=0, we got
(55)PΔ±:b=fΔ±μ1.
We found thatfΔ+Aβ−γd−d2/d=fΔ−Aβ−γd−d2/d and straightforward calculations lead to
(56)fΔ−μ1<HCμ1<fΔ+μ1,μ1∈Aβ−γd−d2d,∞.
Furthermore, we have found that(57)dfΔ±μ1dμ1<0,μ1∈Aβ−γd−d2d,∞andlimμ1⟶∞fΔ−μ1⟶0.d2fΔ±μ1dμ12>0,μ1∈Aβ−γd−d2d,∞andlimμ1⟶∞fΔ−μ1⟶0.
Hence,fΔ±μ1 and HCμ1 are decreasing and they are convex functions of μ1.
Based on the above discussion and Theorem1, for fix k, we define
(58)D0=b,μ1:b>fΔ−μ1,μ1>Aβ−γd−d2d,D1=b,μ1:b>0,μ0<μ1<Aβ−γd−d2d,D2=b,μ1:0<b<fΔ−μ1,μ1>Aβ−γd−d2d.
In Figure1, we see that there is one endemic equilibrium in region D1 and two equilibria in region D2. The system of equations (6)–(9) undergoes saddle-node bifurcation on the curve PΔ−. The forward bifurcation occurs on PΔ− and we have backward bifurcation, which occurs on PΔ−. The pitchfork bifurcation occurs when transversally passing through curve P0 at point K. The mathematical system of equations (6)–(9) has a semihyperbolic node of codimension 2 at point K. The similar discussion can be done for fixed b and varying value of k to show the backward bifurcation. This part is omitted since the analysis is the same as before.☐Figure 1
For fixedk, the bifurcation curves in the μ1,b plane when Aβ>dμ0+γd+d2. In D1 region R0>1, where there exists a unique endemic equilibrium. There are two endemic equilibria in D2. Two equilibria coalesce and a saddle-node bifurcation occurs on PΔ− . The forward bifurcation occurs on P0− and backward bifurcation occurs on P0+. There is no endemic equilibrium in the first region.


## 5. Further Development of the Model, Numerical Results, and Discussion

We numerically show that equilibrium pointE2S∗,I∗,R∗ is locally and globally asymptotically stable. For the parameters in Table 1, R0=6.9307, the endemic equilibrium E2 exists at E2S∗,I∗,R∗=243.9075,2.455,1.6707. Figure 2 provides that Theorem 1 is satisfied, i.e., E2 is locally asymptotically stable where initial conditions S0=150, I0=50, and R0=20 are used. In Figure 2, we could see that solutions approach to E2S∗,I∗,R∗=243.9075,2.455,1.6707. Furthermore, in Figure 3, we use the parameters in Table 1 to show that the endemic equilibrium E2S∗,I∗,R∗=243.9075,2.455,1.6707 is globally asymptotically stable because the solutions of S, I, and R converge to the same E2 independently from the initial values of S, I, and R.Table 1
Random parameters used in the numerical simulations in Figures2 and 3.
ParameterValueDimensionA1.75Individual/timeβ0.011/(individual× time)k21/individuald0.0051/timeμ00.21/timeμ10.31/timeb0.2Individualγ0.21/timeα0.31/timeTable 2
Random parameters used in the numerical simulations for the real data of Saudi Arabia.
ParameterValueDimensionPopulation of SaudiArabia=N35575027IndividualA=nN/3651252Individual/timek0,21/individuald3.6593.10−51/timeμ00.00021/timeμ10.151/timeb100Individualγ0.00281/timeβ00.00000041/(individual× time)Figure 2
Susceptible (S), infected (I), and recovered (R) individuals as functions of time.Figure 3
Global stability of endemic equilibriumE2 for different initial values.We further simulate the inhibition effect due to the behavioral change of susceptible population when the number of infected individuals increases. We use parameter values given in Table1 with k=0,0.2,1,2, and 4. Figure 4 shows different levels of the inhibition effect: low, moderate, and significant. The results are depicted in Figure 5, which shows that moderate and significant inhibition considerably reduces the number of infected individuals (I).Figure 4
Prediction of model (46–49) for a long time,β=2.79.10−7 and R0=1.108.Figure 5
Comparing the prediction of model (46–49 (black line) with real data of Saudi Arabia (red asterisks) forβ=2.79.10−7 and R0=1.108.We now consider a modification of the model defined in equations (6)–(9). We define a new term, which is the vaccination ratio v, and incorporate this term into our existing model. The modified model is shown in Figure 6Figure 6
Schematic diagram of the modified model. The arrows represent a progression from one compartment to the next and death rates.We assume that the vaccination rate is constant, but in general,v is a function of susceptible, infected, and recovered individuals. Using the idea from [21], we can further modify the defined model in (6)–(9) to be
(59)dSdt=A−βIS1+kI−d+vpsS+δR,(60)dIdt=βIS1+kI−μ0+μ1−μ0bb+II−d+γ+γvpiI,(61)dRdt=vpsS+μ0+μ1−μ0bb+II+γvpiI−d+δR,subject to initial conditions(62)S0=S0≥0,(63)I0=I0≥0,(64)R0=R0≥0,where v is the vaccination rate of the population, ps represents the vaccine efficacy for the suspected, and pi is the effectiveness of the vaccination in infected individuals. We could analyze systems (59)–(62) as we did for the original system; therefore, there is no need to duplicate the same analysis here. We apply the modified model to describe COVID-19 scenarios in Saudi Arabia. The total population of Saudi Arabia is 35575027. The total reported cases in Saudi Arabia on April 1, 2021, are 390007, and the active cases are 5452. The Saudi government has provided 4571478 doses of vaccine until April 1, 2021, and 3818608 people of the population are considered to be immune. It is clear from the model that we can derive the reproducing number easily as we explained above using R0=βA/d+vpsd+γ+μ1+γvpi. We now try to simulate Saudi Arabia cases, where pi=0.5,pS=0.95, γ=0.28,v=0.00233, and R0=1.108. We use the least square method to fit the other parameters. Figure 5 shows the model predictions in comparison with real COVID-19 data of Saudi Arabia for infected individuals. The model predicts current cases nicely, but the model prediction in the future cannot be accurate (see Figure 4 after 1000 days). Furthermore, we know from other transmitted diseases that vaccination together with some simple restrictions on the society could control the spread of transmitted diseases. This fact tells us that the current model must be modified as we do below.It is well known that the coefficientsβ, μ0, and μ1 are functions of time in general. But these functions can be written using empirical assumptions by using constant parameters, where empirical relations depend on disease control measures, such as social distancing, partial closure, and tracing of suspected people. In this report, we model these functions as
(65)βt=β0m1t,μ0t=μ0m2t,μ1t=μ1m3t.In general, the form of themit can be considered as in the work of [22, 23] as
(66)mit=m0−m1expk1t−t0+m1t∈t0,λ1,mλ1−m2expk2t−λ1+m2t∈λ1,λ2,−−−−mλp−1−mpexpkpt−λp−1+mpt∈λp−1,∞,where mi measures the intensity of the control measures, ki has a dimension 1/day and simulates the efficiency of the control measures, and λi represents the number of days for each implemented control strategy. The vaccination rate is also a function of time and can be assumed to equal the following equation
(67)vt=VN+numberofNBXday−NDXday−VXday,where V is the number of vaccinations per day, the denominator is the total population× the number of new born per day – new death per day– the number of vaccinations per day. The values of λi, mi, and ki are not fixed; they are different for each country. In this study, we use data from Saudi Arabia. Under the above discussion, our model can be written as
(68)dSdt=A−β0m1tIS1+kI−d+vtpsS+δR,dIdt=β0m1tIS1+kI+−μ0m2t+μ1m3t−μ0m2tbb+II−d+γ+αtvtI,dRdt=μ0m2t+μ1m3t−μ0m2tbb+II−d+δR+αtvtI.Increasing the number of vaccinations per day reduces the reproducing number. Figure7 shows the total infected individuals, where the continuous black line is the model prediction and symbol red asterisk represents the real data of Saudi Arabia (WHO) for the number of vaccinations v=75000 individual/day and R0≤1.311. We now increase the reproducing number up to R0≤1.64, where we increase the transmission rate of infection from 4.10−7 to 5.10−7. Model prediction of the current infected individuals against the real data of Saudi Arabia is shown in Figure 8. We see that the prediction agrees with real data even for the first month unlike the results in Figure 7. The prediction of the model is not a feasible region of real data; from here, we can conclude that the reproducing number of disease for Saudi Arabia is around 1.311. We now consider 50000 vaccinations per day and then consider 1000000 vaccinations per day to show the effect of vaccination as in Figure 9 for β0=4.10−7, where the dashed line represents 100000 vaccinations per day, while the continuous line represents 75000 vaccinations per day and the dashed-dotted line represents 50000 vaccinations per day. We see that increasing the number of vaccinations per day decreases the number of infected individuals as we expected. We also observe from Figure 10 that we could certainly control the spread of COVID-19 after 200 days for more than 75000 vaccinations per day. Next, we show the effect of the reproducing number on the number of infected individuals. Figure 10 shows the effect of the reproducing number on infected individuals. The dashed line, black line, and dashed-dotted line represent the reproducing number R0≤0.98, 1.311, and 1.96, respectively, for fixed vaccinations v=75000. Now, we check the effect of vaccination and reproducing number on the suspected and recovered individuals. Figure 11 shows the effect of vaccination on the suspected population for β0=4.10−7. It is clear that increasing the number of vaccinations reduces the number of suspected individuals. Finally, Figure 12 shows the effect of the vaccination on the recovered population. We find that increasing the number of vaccination increases the number of recovered individuals.Figure 7
Compared the prediction of model (53–55) (black line) with real data of Saudi Arabia (red asterisks) forv=75000 individual/per day and β0=4.10−7, i.e., R0≤1.311.Figure 8
Compared the prediction of model (53–55) (black line) with real data of Saudi Arabia (red asterisks) forv=75000 individual/per day and β0=5.10−7, i.e., R0≤1.64.Figure 9
The effect of the number of vaccination on the infected individuals (dashed line,v=100000, continuous line v=75000, dashed-dotted line, v=50000).Figure 10
The effect of the reproducing number on the infected individuals (dashed line,R0≤0.98, continuous line R0≤1.311, and dashed-dotted line R0≤1.96.Figure 11
The effect of vaccination on the suspected individuals (dashed line,v=100000, continuous line v=75000, and dashed-dotted line, v=50000) for β0=4.10−7.Figure 12
The effect of vaccination on recovered individuals (dashed line,v=100000, continuous line v=75000, dashed-dotted line, v=50000) for β0=4.10−7.


## 6. Conclusions

In this paper, we established a new model with nonstandard nonlinear incidence and recovery rates formulated to consider the impact of available resources of public health, in particular the number of hospital beds and rate of losing immunity. We used Lyapunov’s direct method to show global asymptotic stability of the disease-free equilibrium whenR0<1 and global asymptotic stability of the endemic equilibria when R0>1. We also solved the system numerically, which confirms our theoretical results.Vaccination is an effective method to prevent individuals from contracting transmitted diseases like flu and cholera. In the second part of this paper, we included a vaccination term and found that the modified model agrees well with the real data of Saudi Arabia.In a forthcoming study, we will study a system with an additional vaccinated component such as(69)dVdt=fV,S,I,R,dSdt=gV,S,I,R,dIdt=hV,S,I,R,dRdt=zV,S,I,R.We will compare the prediction of this system with our newly defined system in this paper.



---

*Source: 5207152-2021-10-21.xml*