Dynamics of an SIR-Based COVID-19 Model With Linear Incidence Rate, Nonlinear Removal Rate, and Public Awareness

The paper studies the dynamics of the classical susceptible-infectious-removed (SIR) model when applied to the transmission of COVID-19 disease. The model includes the classical linear incidence rate but considers a nonlinear removal rate that depends on the hospital-bed population ratio. The model also includes the effects of media on public awareness. We prove that when the basic reproduction number is less than unity the model can exhibit a number of nonlinear phenomena including saddle-node, backward, and Hopf bifurcations. The model is fitted to COVID-19 data pertinent to Saudi Arabia. Numerical simulations are provided to supplement the theoretical analysis and delineate the effects of hospital-bed population ratio and public awareness on the control of the disease.


INTRODUCTION
The severity and the global scale of the COVID-19 pandemic have pushed research in many areas including the modeling of the disease dynamics with the goal of using such models to better understand the effects of intervention strategies on the disease control [1][2][3][4][5]. Several factors are known to affect the disease dynamics including the incidence rate, recovery rate, quarantine strategy, and awareness effects. Classical susceptible-infectious-removed (SIR) epidemic models were used extensively in the literature for modeling infectious diseases. These models are structured on the numbers of susceptible (S), infectious (I), and removed (recovered or deceased) individuals (R). When these classical models are used with the standard linear incidence rate and linear recovery rate they are known to predict at most one endemic static equilibrium and cannot exhibit bistability and periodic behavior. In fact, the model dynamic behavior depends mainly on the basic reproduction number. The disease is controlled and subsequently eliminated if the basic reproduction number is less than one and keeps spreading otherwise [6][7][8]. However, in practice infectious diseases can exhibit periodic oscillations and other nonlinear phenomena during the outbreak. Therefore, several studies attempted to model such nonlinear behavior. In this regard, nonlinear incidence rates were shown to play a key role in producing rich dynamics in epidemic models including periodic behavior [9][10][11][12][13].
The recovery rate is also an important parameter in the SIR model. Many of the studied SIR models in the literature assumed a constant recovery rate with the implicit assumption that health resources are sufficient to tackle the disease. However, the COVID-19 pandemic has shown that health services even in the most developed countries were overwhelmed during the first wave of the pandemic. It is therefore more realistic to consider a nonlinear recovery rate that incorporates factors such as hospitals capacity and efficiency of treatment. In this regard, various studies have analyzed epidemic models that incorporated nonlinear expressions of the recovery rate [14,15]. Shan and Zhu [14], in particular, showed that the SIR model with linear incidence rate and nonlinear recovery rate can exhibit rich nonlinear phenomena such as saddle-node, backward, Hopf, and Bogdanov-Takens bifurcations.
Another parameter that can also affect the dynamics of the disease is public awareness and the role played by the media in raising such awareness. The media continuously alert the susceptible individuals regarding the severity of the disease and the need to take preventive measures which can lead to a decrease in the disease transmission [16,17].
In this paper we propose to re-investigate the simple SIR model which is based on the standard incidence rate. We propose two modifications to the model: we assume a nonlinear removal rate that depends on hospital beds availability, and we model the media effects on public awareness in the form of a nonlinear rate with a saturation constant. We show that the proposed model can predict backward and Hopf bifurcations. Numerical simulations are provided to supplement the theoretical proofs. These simulations make use of parameter values fitted to COVID-19 data in Saudi Arabia.
A note should be made about the important issue of model validation. One objective of this work is to identify the parameters of the proposed model using real data pertinent to Saudi Arabia and then use the computed parameter values to analyze the static and dynamic behavior of the model. This parameter estimation is known as the inverse modeling problem. This inverse problem is important for model calibration and to allow public health officials to make reliable predictions about the pandemic. Various approaches were used in the literature for this purpose. These include least squares method that minimizes the sum of squared residuals [2] and variational method applied to time-depended model parameters [18]. Comunian et al. [19] provided a critical analysis of the inverse problem. The inverse problem is known to lack uniqueness of solutions [20] and also to lack the continuous dependence of the parameters to be identified on the data, so that any small errors in the data can lead to large discrepancies in the parameters to be identified. Moreover, the fitting can provide a tool for extracting some key quantities of interest, such as the basic reproduction number, that are not explicit in the epidemiological data streams. Unfortunately, as pointed out rightfully in [21] most epidemiological data streams are not designed for modeling. The challenge for data collection and fitting becomes more severe when the proposed model is highly structured. The balance between model structure and quality of fit was made clear in a number of studies concerning Saudi Arabia, where simple SIR models were shown to provide better fits to COVID-19 data than more complex network-based models [22,23].
The rest of the paper is structured as follows. The model is presented in section 2, then in section 3 we study the model positivity, boundedness of its solutions and derive the basic reproduction number. Section 4 discusses the existence of equilibria while section 5 covers the stability of the disease-free equilibrium. In section 6, we study the bifurcations of the model. Numerical simulations are presented in section 7 followed by concluding remarks.

THE DIMENSIONAL MODEL
The basic SIR model under study is described by the following set of equations: In these equations, S and I denote the numbers of susceptible and infectious individuals. R, on the other hand, represents the removed (recovered or deceased) individuals. In order to simplify the model we did not include a separate equation for the deceased individuals. We assume a linear incidence rate β 1 I with β 1 representing the transmission rate. The term β 2 I M+I represents the decrease in the disease transmission as result of increased public awareness due to media. This awareness can lead the public to take preventive measures. The parameter M represents the half-saturation constant. We assume that β 1 is always larger than β 2 .
A rigorous functional description of public awareness is quite challenging. The expression used to quantify public awareness due to media was adopted from [17], but it is not unique. Other expressions were proposed in literature. Liu et al. [12] modeled the psychological effects of media by assuming the transmission rate β 1 to be an exponential decreasing function that depends on numbers of reported infectious individuals. Feng et al. [16], on the other hand, used a simple constant parameter to quantify the effects of media on public awareness. The parameter γ represents the removal rate. Since the model considers R to represent both the recovered or deceased cases, γ also accounts for the fatality rate i.e., the rate at which infected people decease. Unlike many epidemic models where γ was assumed constant, we assume in this work that γ depends explicitly on the hospital-bed to population ratio (B) and the number of infectious individuals (I). The parameter (B) is generally considered by health authorities to be a key indicator of the adequacy of public health services [24]. The following relation which was first used in [14] is adopted in this work: γ 1 is the value taken by γ (B, I) when B → ∞ and/or the number of infectious individuals is very small i.e., I → 0. γ 1 therefore represents the maximum removal rate when the number of hospital beds is more than sufficient and/or when there are few infectious individuals. γ 0 is, on the other hand, the value taken by γ (B, I) when B → 0 and/or the number of infectious individuals is very large i.e., I → ∞. γ 0 therefore represents the minimum per capita removal rate that can be sustained in the face of a large number of infections. Both values of γ 0 and γ 1 also depend on the inherent characteristics of the specific disease (COVID-19).
In fact, when B is large, then the health system can provide a good recovery rate, so that the removed individuals R includes mostly recovered people. On the other hand, when B is small, then the health system cannot provide care to infected persons, and removed persons are mostly deceased persons. In other words, γ 0 should be close to the "natural" fatality rate.

MODEL ANALYSIS
The model is rendered partially dimensionless by defining: The initial conditions of the system (5-8) are given by: In the following we prove that the proposed model is wellposed by showing that it satisfies the positivity condition.
It can be noted that by adding Equations (5-7) we can deduce that the time derivative of n(t): = s(t) + i(t) + r(t) is always zero and therefore n(t) = 1 for every t and the solutions are therefore always bounded.

Positivity of Solutions
We start by showing that the model solutions are always positive for non-negative initial conditions.
Proof: The first equation (Equation 5) can be written as where Integrating Equation (9) we obtain the following expression: This implies that s(t) is positive for all t. From Equation (6) we have: Integrating Equation (10) we obtain: It can also be shown in the same way that which gives We conclude therefore that the model solutions are always positive.

EXISTENCE AND CLASSIFICATION OF EQUILIBRIA
In this section we study the number of real and positive equilibrium solutions of the system (Equations 5-8). The model has always a disease-free equilibrium given by E 0 (s, i, r) = (1, 0, 0). The non-trivial equilibria are obtained by setting the right sides of equations (Equations 5-7) to zero. The equilibrium condition of (7) is given by r = γ µ i and Equations (5) and (6) are independent of r, so that they can be treated as an independent sub-system. From Equation (6) we have: Substituting Equation (11) into Equation (5) leads to the following cubic polynomial: where a 0 , a 1 , a 2 and a 3 are defined by is the basic reproduction number derived following the techniques presented in [25].
We can see from Equation (12) that the coefficient a 3 is always negative. The coefficient a 0 , on the other hand, is positive when R 0 > 1 and negative when R 0 < 1. The possible number of positive roots can be obtained by the Descartes rule of signs as shown in Table 1.
Theorem 2. The system (5-7) 1. has a unique steady state solution if R 0 > 1 and whenever cases 1, 5, and 7 are satisfied; 2. can have more than one steady state solution if R 0 > 1 and case 3 is satisfied; 3. can have two steady state solutions if R 0 < 1 and cases 2, 4, and 6 are satisfied.
The results of  Figure 1A shows the case of two endemic equilibria for R 0 < 1 and one endemic equilibrium for R 0 > 1. Figure 1B shows, on the other hand, the case of the existence for R 0 < 1 of only the disease-free equilibrium (no endemic solution) and one endemic equilibrium for R 0 > 1. The steady state multiplicity for R 0 < 1 (Figure 1A) indicates the possible existence of a backward bifurcation where a stable endemic equilibrium coexists with the disease-free solution.
Before proceeding to the stability analysis, another useful relation can be obtained for the case of R 0 < 1. The analysis of the coefficients a i (i = 0, 3) (Equation 12) shows that the following inequalities are always satisfied: Therefore, if the last term in Equation (14) is negative i.e., β 1 ≤ µ + γ 0 the other terms of Equation (14) are also negative and therefore a 3 and a 2 are always negative which means that no positive solution can exist under this condition.

STABILITY OF THE DISEASE-FREE EQUILIBRIUM
In this section we study the stability of the disease-free equilibrium point E 0 .
Proof: The Jacobian matrix J(E 0 ) of the system (Equations 5-6) at E 0 is given by The eigenvalues of the matrix Equation (15) Therefore, we can see that if R 0 < 1 then E 0 is an attracting node and is a saddle if R 0 > 1.
In the case of R 0 = 1 the first eigenvalue is zero (λ 1 = 0) and E 0 is a non-hyperbolic point. In order to examine the nature of E 0 we make use of the center manifold theory. We first transform the model (Equations 5-6) into the origin. Let u = s−1 and v = i then we have: Expanding the system (Equations 16-17) in a second order Taylor series around (u, v) = (0, 0) yields When R 0 = 1 then β 1 − γ 1 − µ = 0 and the system (Equations 18-19) becomes: dv dt The Jacobian matrix of the system is with eigenvalues We can therefore rewrite the system (Equations 20-21) as follows: with µm(γ 1 −γ 0 ) then it can be seen from (Equations 24-25) that E 0 is a saddle-node.
Therefore, E 0 is stable if a 11 = 0 and a 22 < 0. Substituting the condition a 11 = 0 into a 22 yields The condition a 22 < 0 requires that otherwise point E 0 is unstable.

Backward Bifurcation
In this section we set up the conditions under which the model exhibits a backward bifurcation. Proof: We examine the existence of a backward bifurcation using the methodology described in [26]. First, we make the transformation of variables as follows: x 1 = s, x 2 = i. The model (Equations 5-6) is rewritten as follows: We assume the bifurcation parameter to be β 1 . The condition R 0 = 1 corresponds to β 1 = γ 1 + µ: = β * 1 . For β 1 = β * 1 , the Jacobian matrix of the system (Equations 28-29) at the disease-free equilibrium E 0 = (1, 0) is given by The right eigenvector of this matrix at the zero eigenvalue is . Similarly, the left eigenvector is given by V = (v 1 , v 2 ) = (0, 1). The stability parameters (a) and (b) that determine the conditions for the existence of a backward bifurcation can be found by the methodology explained in [26]. The coefficient (a) is given by: which is reduced (since v 1 = 0 and v 2 = 1) to The different derivatives evaluated at E 0 are given by Substituting yields On the other hand, the bifurcation parameter b is given by and can be shown to be equal to 1. We conclude therefore that since b = 1 is always positive then the system exhibits a backward bifurcation when a>0. The condition a>0 is equivalent to It can be seen that in the absence of media effects i.e., β 2 = 0 the condition Equation (36) can be satisfied for some values of hospital-bed to population ratio (b). However, in the presence of media effects β 2 = 0 and if the right hand side of the inequality is negative i.e., then the inequality (36) is always satisfied for any values of b. This means that media awareness also plays a critical role in the occurrence/absence of a backward bifurcation. If, on other hand, the inequality (37) is not satisfied then a backward bifurcation occurs when

Hopf Bifurcation
In this section we study the Hopf bifurcation of Equations (5 − −6). The Jacobian matrix is where We use the fact that along the steady state solution we have Thus, we have We can write this equation as follows: with β 2 ) , The coefficients of Equation (44) are always positive when the following condition holds Simple analysis shows that the maximum value of γ 1 when the right side of the condition (45) is positive satisfies The other condition of Hopf point Det(J) > 0 can also be readily derived. However, both Hopf conditions are not easily amenable to analytical manipulations and therefore numerical simulations are used to illustrate the Hopf-induced dynamics of the model.

NUMERICAL SIMULATIONS
In order to carry out numerical simulations we needed first to choose carefully the values of model parameters that reflect the disease parameters in Saudi Arabia, a country that has recorded more than 362,000 cases as on December 25, 2020 [27]. Only two model parameters were assumed constant: the natural death rate µ which is 4.21 × 10 −5 [28] and the hospital-bed to population ratio which is 22 per 10,000 [29]. The rest of model parameters (m, β 1 , β 2 , γ 0 , γ 1 ) were fitted to COVID-19 data in the country.
The fitting was carried out with the help of the optimization function (fmincon) of MATLAB. The COVID-19 data used for fitting consist of infected cases and removed (recovered and deceased) cases. The data was collected for a period of 90 days starting from March 2 (the day of the detection of the first COVID-19 case in the country) until June 2. The data is publicly accessible via the COVID-19 dashboard of the Saudi Ministry of Health [27]. The reason we restricted the data to this period is to be able to compare the obtained reproduction number with similar results available in the literature [30,31]. Figure 2 shows the fitting results. Both the infected cases (Figure 2A) and removed cases ( Figure 2B) show a reasonable quality of fit. Since the optimization task led to local optima, we were guided by the results reported in [22,30,31] in selecting initial conditions and bounds on the fitted model parameters. The optimum values of model parameters are shown in Table 2.
For bifurcation studies the basic reproduction number R 0 was taken to be the main bifurcation parameter so the transmission rate β 1 can be calculated by β 1 = R 0 (γ 1 + µ) (Equation 13). The bifurcation diagrams were obtained using the continuation software AUTO [32]. Using the obtained values of model parameters, the condition (Equation 38) for a backward bifurcation to occur can be shown to be b < 2.01 × 10 −4 . Figure 3 shows an example of bifurcation diagrams in the (R 0 − i) plane for b = 6.322 × 10 −7 . This is the actual value of the dimensionless hospital-bed to population ratio, and can be seen to be much smaller than the aforementioned critical value. Saddle-node bifurcation can be seen to occur in this case. A static limit point LP appears at R 0 = 0.803. The diagram is also characterized by the presence of a Hopf point occurring at R 0 = 0.896. An unstable limit cycle can be seen to emerge from the Hopf point and disappears from homoclinic bifurcation. Here we have a situation where the disease-free equilibrium coexists with the endemic equilibrium for the range of basic reproduction number extending from HB (R 0 = 0.896) to (R 0 = 1). Figures 4A,B show numerical simulations for the Hopf point and the limit cycle. Figure 4A shows the birth of a limit cycle close to the Hopf point while Figure 4B shows an example of time trace of an unstable limit cycle at (R 0 = 0.899). The oscillations grow in magnitude until they terminate as they collide with the static branch (S = 1).
From an epidemiological perspective, the existence of a Hopf point in the model and the nature of the Hopf bifurcation (subcritical) has two effects on the dynamics of the disease transmission: (1) The endemic equilibrium coexists now with the disease-free solution for values of R 0 ranging from the Hopf point to R 0 = 1. This means that even if the basic reproduction number is reduced below unity, it is still possible for the pandemic to persist unless the reproduction number is reduced below the point associated with the Hopf point (HB, R 0 = 0.896). (2) The existence of a subcritical Hopf bifurcation (as opposed to supercritical when a stable limit  cycle surrounds an unstable equilibrium point) means that the pandemic cannot have a stable periodic behavior and this is consistent with the current situation in the country which showed a steady decline in the disease active cases. Even if initial conditions (S, I) would fall close to the limit cycle, the resulting unstable oscillations would increase in magnitude but they will ultimately collapse and lead to the disease suppression ( Figure 4B). Figure 5 shows, on the other hand, the bifurcation diagram for b = 2.5 × 10 −4 , a value just larger than the aforementioned critical point. The backward bifurcation disappears and we have the occurrence of a forward bifurcation where the endemic equilibrium is the only stable point. A forward bifurcation is of course the outcome to aim for in a pandemic since if the basic reproduction number is reduced to a value below unity, then the disease outbreak is suppressed. Figures 3, 5 therefore showed the important role the number of hospital beds plays in the control of the pandemic. However from a practical point of view, the current hospital-bed to population ratio (22 per 10,000 people) is much smaller than the critical value. This means that a forward bifurcation is unlikely.
Next, we examine the effects of variations in model parameters on the occurrence of the backward bifurcation. This can be suitably done by showing the location of the Hopf point when these parameters are varied, as shown in Figures 6A-E obtained also using AUTO. These figures are quite useful since they illustrate the range of the backward bifurcation, that is the range in term of R 0 of the stable coexistence of the endemic equilibrium with the disease-free equilibrium. This sensitivity analysis is also necessary given the unavoidable degree of uncertainty in the fitted values of model parameters.
We start with the effects of media shown in Figures 6A,B. The media effects are represented by two parameters: the awareness rate β 2 and the half-saturation constant m. It can be seen from Figure 6A that the range of backward bifurcation decreased in term of R 0 since the Hopf point moves to large values. The effect of m ( Figure 6B) shows that as m increases the Hopf point moves to larger values of R 0 which means again a smaller range of coexistence with the disease-free equilibrium. We can therefore conclude that the media also have an important role in attenuating the backward bifurcation. Figure 6C shows, on the other hand, the effects of the minimum removal rate γ 0 . It can be seen that larger values of γ 0 decrease the range of the coexistence since the Hopf point moves to larger values and closer to 1. The effect of the maximum removal rate γ 1 is shown in Figure 6D. Unlike γ 0 , it can be seen that as γ 1 decreases, the Hopf point moves to larger values of R 0 . While it is true that the values of γ 0 and γ 1 depend on the inherent properties of the virus causing the COVID-19 disease, they also depend on health care and basic clinical resources [14]. Their role is therefore important in attenuating the backward bifurcation. Finally, Figure 6E shows that as the hospital-bed to population ratio b increases, the location of the Hopf point moves to larger values of R 0 . This confirms the already obtained results that the backward bifurcation is reduced with the increase in available hospital beds.

CONCLUSIONS
The paper examined the dynamics of a SIR model applied to COVID-19 disease with two modifications: (1) the removal rate was assumed to be nonlinear and function of both the number of infectious individuals and hospital-bed to population ratio; (2) the model included public awareness due to media. The model was validated using COVID-19 data pertinent to Saudi Arabia. The model was shown to exhibit backward and Hopf bifurcations. The Hopf point is subcritical and unstable limit cycles terminate homoclinically. Therefore no stable periodic behavior occurred in the model for the chosen set of parameters. The main behavior observed in the model is a static coexistence of the disease-free equilibrium with the endemic point for a range of the basic reproduction number extending from the Hopf point to R 0 = 1. The analysis also showed the important role played by health care resources in term of hospital-bed to populations ratio as well as the minimum and maximum removal rates. In fact with the country's current health resources, a backward bifurcation is the most likely outcome and can be attenuated by the increase in the media awareness.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.