Study of Dynamics of a COVID-19 Model for Saudi Arabia with Vaccination Rate, Saturated Treatment Function and Saturated Incidence Rate

: This paper proposes, validates and analyzes the dynamics of the susceptible exposed infectious recovered (SEIR) model for the propagation of COVID-19 in Saudi Arabia, which recorded the largest number of cases in the Arab world. The model incorporates a saturated incidence rate, a constant vaccination rate and a nonlinear treatment function. The rate of treatment is assumed to be proportional to the number of infected persons when this number is low and reaches a ﬁxed value for large number of infected individuals. The expression of the basic reproduction number is derived, and the model basic stability properties are studied. We show that when the basic reproduction number is less than one the model can predict both a Hopf and backward bifurcations. Simulations are also provided to ﬁt the model to COVID-19 data in Saudi Arabia and to study the effects of the parameters of the treatment function and vaccination rate on disease control.


Introduction
Multiple waves of COVID-19 are still attacking many countries around the world. The severity of the pandemic is accentuated by the emergence of new and more contagious strains of the virus. Many countries are still struggling to vaccinate their population in order to curve the waves of infections and reduce the hospitalizations that have affected their health systems. The use of mathematical models to predict disease transmission is useful for analysis and prediction particularly when the model is validated against the available disease data.
A mathematical model can predict future situations of the disease spread and evaluate the best strategy to be invoked depending on the epidemiological situation. Since the outbreak of the disease, there have been a large number of studies [1][2][3][4][5][6][7][8][9] that used compartment models based on SIR (susceptible (S), infectious (I) and recovered (R)) and SEIR (susceptible (S), exposed (E), infectious (I) and recovered (R)) epidemic models to predict the spread of the disease.
These models consist generally of deterministic ordinary differential equations or difference equations and are amenable sometimes to analytical manipulations that can extract some important parameters, such as the basic reproduction number that can quantify the spread of the disease. These models are also flexible. They can accommodate a variety of expressions for incidence rate [1], can include the effect of a variable recovery rate [2] and can simulate various scenarios of pandemic management, such as lockdown [3,4], closure of sectors of the economy, compliance of people with protection measures [4], social distancing [5], travel restrictions [6], media coverage [7], quarantine [8], super-spreaders individuals [9] and vaccination [10][11][12]. Once validated, these models can be used as a public health decision support tool [4].
In this paper, we develop, validate and examine the dynamics of a SEIR model for COVID-19. The proposed model includes a number of characteristics. In many epidemic models, the bilinear incidence rate βSI is used. In this paper, we adopt a nonlinear expression of the incidence rate βSI/(1 + αI) that includes an inhibition 1/(1 + αI) effect reflecting the behavior of susceptible persons as the infection spreads. This form of incidence rate was shown to be more realistic than the bilinear incidence expression [1]. The proposed model also includes a treatment function. A number of epidemic models considered the treatment rate to be proportional to the number of infected persons. However, the recovery rate is known to be strongly dependent on health resources and the treatment efficiency. Since, generally, the treatment capacity of each country is limited, it seems logical to consider a more realistic treatment function. Various expressions for treatment functions have been proposed in the literature [13][14][15][16][17]. Here, we adopt the saturated treatment cI/(1 + bI) function suggested by Zhang and Liu [17]. This form has the advantage of producing a constant value when the number of infected persons (I) is very low and reaches an asymptotic constant value if the number of infected individuals is large. Furthermore, this form of treatment function has a finite and continuous value. The model also includes a vaccination rate. Vaccine administration is a highly effective method of preventing and reducing viral infections [18]. Vaccination and optimal control are the key points to controlling an epidemic situation as discussed in [11,[19][20][21], and the majority of countries have embarked on aggressive vaccination campaigns against COVID-19.
The novelty in the proposed SEIR model is that it includes a saturated incidence rate and a saturated treatment function in addition to a vaccination rate. Some models in the literature either included some of these elements or included all of these elements but not in the current expressions adopted in this paper. Zhang et al. [16], for instance, included, in their SEIR model, all the elements except the vaccination rate. Annas [10], in their study of the spread of the disease in Indonesia included, in their SEIR model, all the aforementioned elements except the treatment function. Deng et al. [22], proposed a non-smooth Filippov epidemic system to examine the effects of some control strategies, such as media coverage, treatment and vaccination. However, the authors considered only a SIR and not a SEIR model.
In the first step of our analysis, the model is validated using real COVID-19 data in Saudi Arabia. This allows the extraction of a number of model parameters. In the next step, simulations are carried out to study the bifurcation behavior of the model using the basic reproduction number as the main bifurcation parameter. A numerical sensitivity analysis is also performed for the effect of model parameters associated with the treatment function and the vaccination rate.
The organization of the rest of the paper is as follows. In Section 2, we present and explain the model. Then, in Section 3, we prove the model positivity and boundedness of its solutions and derive the expression of the basic reproduction number. In Section 4, the backward bifurcation is examined. Section 5 of the paper includes numerical simulations to support the theoretical proofs.

The Model
The proposed SEIR model for COVID-19 transmission is composed of four compartmentsnamely, susceptible (S), exposed (E), infected (I) and recovered individuals (R). The changes that occur in each compartment can be interpreted by Figure 1. The model consists of the following ordinary differential equations, A is the population recruitment rate, µ is the population natural death rate per time, and α is the transmission rate. The term 1 1+α 1 I reflects the inhibitory effect from the behavioral change of the susceptible when the number of incidence increases. β is the rate of transformation to infective persons, δ is the natural recovery rate of the infected persons, ε is the mortality due to the disease, and v is the rate of vaccination of susceptible population. The term T(I) = cI (1+bI) represents the treatment function where b represents the saturation factor measuring the effect of the infected individual when delayed for treatment, and c is the maximal supplied medical resources. T(I) can be seen to approach cI when I is small; On the other hand, T(I) approaches c/b when I is large. It should be noted that, compared to other forms proposed in the literature [13][14][15], this form of treatment function has also the advantage of being continuous and differential.

Positivity of the Model Solutions
Proof. The proof follows the general approach used by a number of authors, including [23]. We have from (Equation (1)) that, where Integrating Equation (5) yields Thus, S(t) > 0 for all t.
For the proofs that E(t) and I(t) are positive, we divide this issue into four cases: (2) and (3), we can see is continuous at t = 0 and since dI dt (t = 0) = βE 0 > 0, we can conclude that E(t) > 0 and I(t) > 0 for all t ≥ 0. If this is not true, then we can choose which contradicts the assumption that E(t 1 ) = 0. If, on the other hand, I(t 1 ) = 0, then there exists θ such that t 1 > θ and 0 This gives which also contradicts the assumption that I(t 1 ) = 0. The same analysis can be carried out for the other cases: (E 0 = 0, I 0 > 0) and (E 0 > 0, I 0 > 0). We conclude, therefore, that E(t) and I(t) are positive for all t ≥ 0. Next, since we have shown that S(t) > 0 and I(t) > 0, we can deduce from Equation (4) Applying Grönwall's lemma [24] to Equation (9) with initial condition R 0 yields, This proves that R(t) > 0.
Proof. Denote by N(t) = S(t) + E(t) + I(t) + R(t), and then we have Therefore, where N(0) is the initial condition of N(t). Therefore, 0 < N(t) < A µ as t goes to +∞ and This shows that Ψ is positively invariant.

Existence of Equilibria, Stability and Bifurcation Analysis
In the rest of analysis, we consider only Equations (1) where I is the positive root of the quadratic equation, where a 0 , a 1 and a 2 are defined by In these equations, R 0 is the basic reproduction number, obtained with the techniques in [25]: Let ∆ be the discriminant of the quadratic Equation (16). Solving for ∆ = 0 in terms of R 0 , we find R 0 = R c 0 , where The existence of the equilibrium can be analyzed by the following results: 1. Case b = 0: Equation (16) has a unique solution I = −a 0 /a 1 . In this case, the Equations (1)-(3) has a unique steady state solution whenever R 0 > 1 and no steady state solution if R 0 ≤ 1.

Local Stability Analysis of the Disease-Free Solution
In this part of the paper, the local stability of the disease-free solution is studied using the eigenvalues of the Jacobian matrices of (1)-(3) at the disease-free equilibrium E 0 . Thus, we have The eigenvalues of Jacobian matrices J(E 0 ) are: with By simple Algebraic calculation, we have λ 2 = 0 when R 0 = 1, λ 2 < 0 for R 0 < 1 and λ 2 > 0 when R 0 > 1. Therefore, the following result is obtained: Lemma 1. The disease-free equilibrium E 0 is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1. Proof. We consider the methodology detailed in [26] to prove the existence of a backward bifurcation. We define new variables: x 1 = S, x 2 = E, x 3 = I. We can rewrite Equations (1)-(3) as,

Backward Bifurcation
Let us assume α to be the bifurcation parameter. The condition R 0 = 1 leads to The Jacobian matrix of (21)-(23) at E 0 when α = α * is The right eigenvector of the Jacobian matrix for the zero eigenvalue is The left eigenvector is (v 1 , v 2 , v 3 ) = (0, β µ+β , 1). The two stability parameters a and b [26] that define the conditions for the appearance of a backward bifurcation are: which yields The other parameter b is: which is reduced to Since b is always positive, a backward bifurcation exists when

Numerical Simulations
The first step in numerical simulations is to extract fitted values of model parameters using COVID-19 data for Saudi Arabia, a country with a population of 34,800,000 that reported 543,000 cases by the first of September 2021 [27]. The following model parameters were assumed to have constant and known values. These are the recruitment rate A = 1252 and the natural death rate µ = 4.21 × 10 −5 [28]. The death rate due to the disease is taken to be ε = 0.012 [27]. The other model parameters (α 1 , β, δ, b, c) were fitted to the COVID-19 data, which are freely available on the Saudi Ministry of Health website [27].
The fitting task was achieved using (FMINCON), a MATLAB optimization software. FMINCON is used to find minimum of a constrained nonlinear multivariable function using a sequential quadratic programming (SQP) method. The fitting was done for a six month period starting from 24 March 2020. During this period, the vaccination campaign had not started yet, and thus the vaccination rate was taken to be zero. Figure 2 shows the results of fitting both infectious cases (I) and recovered cases (R). The fitting seems reasonable especially after an initial period of time. Table 1 summarizes the obtained values of the fitted parameters. Using the fitted values (Table 1) Next, numerical results are shown for the bifurcation behavior of the model using R 0 as the main bifurcation parameter. The numerical analysis of the model is carried out using the software AUTO [29]. AUTO can be used to solve continuation and bifurcation problems in ordinary differential equations. These include continuation of equilibria and periodic orbits. First, it can be noted that using the values of Table 1 yields the following condition b > b cr = 0.0186 (Equation (30)) for the emergence of a backward bifurcation. Figure 3a illustrates the bifurcation behavior for b = 0.0376, which is the actual fitted value and which is larger than the critical value. On the diagram, the occurrence of a static limit point is at R 0 = 0.62721, and a Hopf point is at R 0 = 0.62722. The limit point and Hopf point are very close to each other. Normally (not shown in the figure), an unstable periodic branch would bifurcate from the subcritical Hopf point and collide homoclinically with the static branch. However, the details are too small to have a practical impact given the closeness of the Hopf and limit points. It can be concluded from the diagram that the disease-free solution coexists with the stable endemic solution for R 0 extending from HB (R 0 = 0.62722) to (R 0 = 1). This diagram shows how important the role the parameter b of the treatment function is in controlling the disease (b accounts for the infected individuals when their treatment is delayed). For values of b larger than the critical value, the eradication of the disease does not happen when the reproduction number is reduced below one, rather it is necessary to reduce it to values below the R 0 value corresponding to the Hopf point. Figure 3b shows the bifurcation behavior when b is chosen smaller than the critical value, i.e., b = 0.018. The backward bifurcation no longer exists and, instead, there is the appearance of a forward bifurcation, and the disease-free equilibrium is the only stable solution for R 0 < 1. The disease can be suppressed if the basic reproduction number falls below the value of one.  The vaccination rate does not greatly affect the backward bifurcation, but it does affect the time evolution of the disease. Figure 5 shows the effect of the vaccination rate v taking the values of 0, 0.0003, 0.0006, 0.0012 to 0.0024. It can be observed that, in the window of one year, the effect of doubling the vaccination rate each year has a highly nonlinear effect on the number of susceptible and recovered individuals. This shows that the giving of vaccines to the susceptible group is influential in increasing the total population of recovered COVID-19 patients. These results are pertinent to the vaccination rate and in line with previous work carried out in Indonesia using a simpler model that did not consider a treatment function [10].

Conclusions
This paper studied the dynamical behavior of a SEIR model for the propagation of COVID-19. The model included a constant vaccination rate and a nonlinear treatment function in addition to a saturated incidence rate. The dynamics of the model were analyzed, and the basic reproduction number was derived. The model was fitted to the COVID-19 data of Saudi Arabia.
The main behavior found in the model was Hopf and backward bifurcations. The analysis also showed that the capacity of the treatment and the saturation factor measuring the effect of the infected individual, when delayed for treatment, played important roles in the existence of backward bifurcation. The remaining model parameters, including the vaccination rate, had a minimal effect on the occurrence of backward bifurcation. The vaccination rate, on the other hand, had an important effect on the time evolution of the disease.