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

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.


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][2][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 susceptibleinfective-recovered (SIR) model, where the total population is divided into three classes: susceptible, infective, and recov-ered [4]. After that, variant of SIR compartmental models were developed, some of them outlined in [5][6][7][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 NðtÞ = SðtÞ + IðtÞ + RðtÞ is the number of total populations, parameter β is the infection transmission rate, and SðtÞ, IðtÞ, and RðtÞ 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 R 0 ; the disease will be eliminated if R 0 < 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 f ðIÞS (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: 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 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 Section 2, 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.

Model Framework
In this section, we describe the mathematical formulation of an SIR epidemic model, where the total population size of individuals is represented by NðtÞ. The total population size is subdivided into three different classes, namely, SðtÞ: susceptible, IðtÞ: infected, and RðtÞ: 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 time t 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 The infected population is decreased by the natural death rate d 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: The recovered population is increased by the recovery rate as in above and decreased by the natural death rate d and the susceptibility of recovered individuals δ. The time rate of change for the population of recovered individuals can be represented by the following equation: All parameters involved in equations (3)-(5) are nonnegative. Equations (3)-(5) can be written in the combined form as follows: Journal of Function Spaces subject to initial conditions For systems (6)-(9), the cone R 3+ is positively invariant. The C 1 smoothness of the right side of systems (6)-(9) implies local existence and uniqueness of solutions with the initial values in R 3+ .
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 which implies that the set is 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 population NðtÞ = SðtÞ + IðtÞ + RðtÞ 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.

Basic Reproduction Number
The basic reproduction number is denoted by R 0 , 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: where We denote the disease-free equilibrium of models (6)-(9) by E 0 , where Now, the Jacobian matrixes of F i and V i at E 0 are given as where The basic reproduction number R 0 of system (14) is defined by the spectral radius of the matrix ρðFV −1 Þ (see Driessche and Watmough [19]) and it is given by R 0 = βA/ ðdðd + γ + μ 1 ÞÞ.

Existence and Types of Equilibria
4.1. Existence of Equilibria. We have already established the existence of the disease-free equilibrium E 0 = ðA/d, 0, 0Þ for systems (6)- (9) for all values of parameters. For any endemic equilibrium E * = ðS, I, RÞ, its coordinates satisfy Substituting this into the first equation in (6)-(9) and solving for S, we obtain 3 Journal of Function Spaces The variable I should be the positive root of the following quadratic equation where and Equation f ðIÞ = 0 may have two roots if Δ 0 > 0, which are , where Δ 0 = C 2 − 4AD.
Proof. The Jacobian matrix JðE 0 Þ for systems (6)-(9) is given by The eigenvalues are obtained from this.
Proof. Consider the Lyapunov function V = I in R 3+ with the Liapunov derivative 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 E 0 . Therefore, E 0 is globally asymptotically stable.☐ Theorem 4. A sufficient condition for the endemic equilibrium E 2 to be locally asymptotically stable is R 0 > 1 Proof. In this case, the Jacobian matrix has the following form: where Σ = ð−βS * /ð1 + I * kÞÞ + ðβkI * S * /ðð1 + I * kÞ 2 ÞÞ and Ξ = ððμ 1 − μ 0 ÞbÞ/ðb + I * ÞððI * /ðb + I * ÞÞ − 1Þ − μ 0 . The characteristic equation for the above Jacobian around its endemic equilibrium E 2 is where C 1 , C 2 , and C 3 are too long to reproduce here. We have observed that C 1 , C 2 , and C 3 are positive and C 1 C 2 > C 3 ; hence, the Routh-Hurwitz criterion is satisfied, so systems (6)-(9) are locally asymptotically stable for R 0 > 1.☐ Theorem 5. The epidemic models (4)- (7) at E 2 are globally asymptotically stable if R 0 > 1:

Journal of Function Spaces
Proof. We consider the Lyapunov function given by Taking the derivative, we have Then, we have Hence, the theorem is proved.☐ Theorem 6. For systems (6)- (9), consider R 0 as the bifurcation parameter. Then, we have, when R 0 = 1, systems (6) Proof. Since R 0 is a function of the parameters β, γ, δ, d, and μ 1 , without loss of generality, we can choose μ 1 as the bifurcation parameter.
Since R 0 is a function of μ 1 , also b and k, we choose μ 1 as the bifurcation parameter. In Theorem 1, there are two endemic equilibria provided that if R 0 < 1, C < 0, and Δ 0 > 0 and if C < 0 and Δ 0 = 0, then, there are two equilibria coalesce.
Based on the above discussion and Theorem 1, for fix k, we define In Figure 1, we see that there is one endemic equilibrium in region D 1 and two equilibria in region D 2 . 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 P 0 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.☐

Further Development of the Model, Numerical Results, and Discussion
We numerically show that equilibrium point E 2 ðS * , I * , R * Þ is locally and globally asymptotically stable. For the parameters in Table 1, R 0 = 6:9307, the endemic equilibrium E 2 exists at E 2 ðS * , I * , R * Þ = ð243:9075,2:455,1:6707Þ. Figure 2 provides that Theorem 1 is satisfied, i.e., E 2 is locally asymptotically stable where initial conditions Sð0Þ = 150, Ið0Þ = 50, and Rð0Þ = 20 are used. In Figure 2, we could see that solutions approach to E 2 ðS * , 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 E 2 ðS * , I * , R * Þ = ð 243:9075,2:455,1:6707Þ is globally asymptotically stable because the solutions of S, I, and R converge to the same E 2 independently from the initial values of S, I, and R. 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 Table 1 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).
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 6 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 dR Figure 1: For fixed k, the bifurcation curves in the ðμ 1 , bÞ plane when Aβ > dμ 0 + γd + d 2 . In D 1 region R 0 > 1, where there exists a unique endemic equilibrium. There are two endemic equilibria in D 2 . Two equilibria coalesce and a saddle-node bifurcation occurs on P − Δ . The forward bifurcation occurs on P − 0 and backward bifurcation occurs on P + 0 . There is no endemic equilibrium in the first region. 8 Journal of Function Spaces subject to initial conditions where v is the vaccination rate of the population, p s represents the vaccine efficacy for the suspected, and p i 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 R 0 = βA/ððd + v p s Þðd + γ + μ 1 + γvp i ÞÞ. We now try to simulate Saudi Arabia cases, where p i = 0:5, p S = 0:95, γ = 0:28, v = 0:00233, and R 0 = 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    Journal of Function Spaces 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 In general, the form of the m i ðtÞ can be considered as in the work of [22,23] as where m i measures the intensity of the control measures, k i has a dimension 1/day and simulates the efficiency of the    Journal of Function Spaces 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 where V is the number of vaccinations per day, the denominator is the total population × the number of new born per daynew death per day-the number of vaccinations per day. The values of λ i , m i , and k i 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

Journal of Function Spaces
Increasing the number of vaccinations per day reduces the reproducing number. Figure 7 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 R 0 ≤ 1:311. We now increase the reproducing number up to R 0 ≤ 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 R 0 ≤ 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.

Conclusions
In this paper, we established a new model with nonstandard nonlinear incidence and recovery rates formulated to con-sider 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 when R 0 < 1 and global asymptotic stability of the endemic equilibria when R 0 > 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 We will compare the prediction of this system with our newly defined system in this paper.

Data Availability
The authors confirm that the data supporting the findings of this study are available within the article and/or its supplementary materials.

Conflicts of Interest
The authors declare that they have no competing interests.

Authors' Contributions
Fehaid Salem Alshammari did the conceptualization, methodology, numerical simulations, and writing of the original