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 () 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 . The model has a unique endemic equilibrium when . 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 , 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.
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 (where is the number of total populations, parameter is the infection transmission rate, and represent the number of susceptible and infected and recovered individuals at time . Also, a linear recovery rate ( 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 ; the disease will be eliminated if ; 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 (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 measures the psychological or inhibitory effect and and 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 .
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 are the minimum and maximum per capita recovery rates, respectively. Parameter 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.
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 by . The total population size is subdivided into three different classes, namely, : susceptible, : infected, and : 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 is and all the new recruited populations go to the susceptible class. We consider that , which is the disease transmission rate and the incidence rate to be , where is the half-saturation constant for which the susceptible population decreases. The population of susceptible people is decreased by the natural death rate , 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 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 , the number of hospital beds, and the number of infectives . 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 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: subject to initial conditions
For systems (6)–(9), the cone is positively invariant. The smoothness of the right side of systems (6)–(9) implies local existence and uniqueness of solutions with the initial values in .
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 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 by , 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 , where
Now, the Jacobian matrixes of and at are given as where
The basic reproduction number of system (14) is defined by the spectral radius of the matrix (see Driessche and Watmough [19]) and it is given by .
4. Existence and Types of Equilibria
4.1. Existence of Equilibria
We have already established the existence of the disease-free equilibrium for systems (6)–(9) for all values of parameters. For any endemic equilibrium , its coordinates satisfy
Substituting this into the first equation in (6)–(9) and solving for , we obtain
The variable should be the positive root of the following quadratic equation where
and
Equation may have two roots if , which are where .
If , we see from (19) and (20) that and , respectively.
We study the existence of equilibria in the following three cases:
In this case, ; since we have and , so systems (6)–(9) have a unique endemic equilibrium .
In this case, and ; if , then, systems (6)–(9) have a unique endemic equilibrium .
In this case, if , there is no endemic equilibrium; if and , we have two endemic equilibria and . If , 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)If , there exists a unique endemic equilibrium (3)If there exists a unique endemic equilibrium provided that ; otherwise, there is no endemic equilibrium(4)If and if , there is no endemic equilibrium; if and , the system has two endemic equilibria and ; if and , the system has two equilibria coalesce into
Theorem 2. For systems (6)–(9), the disease-free equilibrium is (i) an attracting node(ii)a hyperbolic saddle(iii) and : a saddle-node of codimension 1 is a saddle-node of codimension 1; is an attracting semihyperbolic node of codimension 2.
Proof. The Jacobian matrix for systems (6)–(9) is given by The eigenvalues are obtained from this. Since , is negative, because ; hence, all eigenvalues are negative. Thus, if , is an attracting node, and if , then, is a hyperbolic saddle. If , the third eigenvalue is zero. In order to determine the type of , we first transform the disease-free equilibrium point to the origin. We use ; then, using the Taylor expansion as in [15], we get Using eigenvectors, we find following transform: Thus, using these variables in (30), we obtain It is unnecessary to calculate the center manifold if , and is a saddle node. If , then, from the central manifold theorem for systems (6)–(9), we obtain Hence, is a semihyperbolic attracting node.☐
From Theorem 1, one can see that for the existence of endemic equilibrium, is a necessary condition. If , is the unique equilibrium point and we have the following theorem.
Theorem 3. For the system (4)-(7), if , the disease-free equilibrium is globally asymptotically stable.
Proof. Consider the Lyapunov function in with the Liapunov derivative The Lyapunov–Lasalle theorem implies that solutions in approach the largest positively invariant subset of set , i.e., plane . In this plane, and as . Thus, all solutions in plane go to the disease-free equilibrium . Therefore, is globally asymptotically stable.☐
Theorem 4. A sufficient condition for the endemic equilibrium to be locally asymptotically stable is
Proof. In this case, the Jacobian matrix has the following form: where and . The characteristic equation for the above Jacobian around its endemic equilibrium is where , and are too long to reproduce here. We have observed that , and are positive and ; hence, the Routh–Hurwitz criterion is satisfied, so systems (6)–(9) are locally asymptotically stable for .☐
Theorem 5. The epidemic models (4)–(7) at are globally asymptotically stable if
Proof. We consider the Lyapunov function given by Taking the derivative, we have if Then, we have Hence, the theorem is proved.☐
Theorem 6. For systems (6)–(9), consider as the bifurcation parameter. Then, we have, when , systems (6)–(9) which undergo forward bifurcation if ; systems (6)–(9) undergo backward bifurcation if ; systems (6)–(9) undergo pitchfork bifurcation if .
Proof. Since is a function of the parameters , without loss of generality, we can choose as the bifurcation parameter.
Let , we substitute this into (6)–(9); we note that if , then, it reduces to . We then use Taylor expansion at ; diagonalizing the linear part, we then apply the center manifold theorem for the parameter . We found that
where . Denoting the right-hand side of (47) as , we have
Therefore, the system has transcritical bifurcation if
If , using the center manifold theory, as the above, from (6)–(9), we obtain
Again, denoting the right-hand side of (50) as , we have
It is clear that systems (6)–(9) have pitchfork bifurcation if when [20]. We note that If , the system related with equation (47) and the system related with equation (50) reduce to the systems in equations (34) and (37), respectively.
Since is a function of , also , we choose as the bifurcation parameter. In Theorem 1, there are two endemic equilibria provided that if , , and and if and , then, there are two equilibria coalesce.
The basic reproducing number defines a straight line in the plane
=0 defines one branch of hyperbola (see Figure 1):
The branch of has an intersection with at point K, where . We also have , so is decreasing and a convex function of .
Now, let
Then, , where are the parts of which are separated by K. Now, we consider the curve defined by and we indicate these curves as ; solving for from , we got
We found that and straightforward calculations lead to
Furthermore, we have found that
Hence, and are decreasing and they are convex functions of .
Based on the above discussion and Theorem 1, for fix , we define
In Figure 1, we see that there is one endemic equilibrium in region and two equilibria in region . The system of equations (6)–(9) undergoes saddle-node bifurcation on the curve . The forward bifurcation occurs on and we have backward bifurcation, which occurs on . The pitchfork bifurcation occurs when transversally passing through curve at point . The mathematical system of equations (6)–(9) has a semihyperbolic node of codimension 2 at point . 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.☐
5. Further Development of the Model, Numerical Results, and Discussion
We numerically show that equilibrium point is locally and globally asymptotically stable. For the parameters in Table 1, , the endemic equilibrium exists at . Figure 2 provides that Theorem 1 is satisfied, i.e., is locally asymptotically stable where initial conditions , , and are used. In Figure 2, we could see that solutions approach to . Furthermore, in Figure 3, we use the parameters in Table 1 to show that the endemic equilibrium is globally asymptotically stable because the solutions of , , and converge to the same independently from the initial values of , , and .
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 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 ().
We now consider a modification of the model defined in equations (6)–(9). We define a new term, which is the vaccination ratio , 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, 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
subject to initial conditions where is the vaccination rate of the population, represents the vaccine efficacy for the suspected, and 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 . We now try to simulate Saudi Arabia cases, where , and . 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 , and 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 can be considered as in the work of [22, 23] as where measures the intensity of the control measures, has a dimension and simulates the efficiency of the control measures, and 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 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 , and 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
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 individual/day and . We now increase the reproducing number up to , 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 , 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 0.98, 1.311, and 1.96, respectively, for fixed vaccinations . 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 . 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.
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 when and global asymptotic stability of the endemic equilibria when . 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 draft. F. Talay Akyildiz did the conceptualization, methodology, numerical simulations, and writing of the original draft.
Acknowledgments
This research was supported by the Deanship of Scientific Research, Imam Mohammad Ibn Saud Islamic University (IMSIU), Saudi Arabia, Grant no. 21-13-18-004.