Mathematical COVID-19 model with vaccination: a case study in Saudi Arabia

Introduction

The outbreak of several pandemics such as COVID-19 requires the development of mathematical models in order to exhibit key epidemiological features, investigate transmission dynamics, and develop adequate control policies. Mathematical modelling when dealing with infectious diseases allows revealing inherent patterns and underlying structures that govern outbreaks. Simple models that contain the essential components and interactions are powerful tools to test different hypotheses and understand disease control for both short and long time. The stability analysis near the free disease equilibrium will show if the apparition of new infection cases will yield to disease outbreak. Some countries such Tunisia and Jordan registered zero cases for days in Summer 2020 but the introduction of new cases resulted in critical endemic situation by Autumn.

The complex spreading patterns of COVID-19 and the various spread speed of its variants make its containing and mitigating real challenges. The existing models vary in form and complexity, but the common objective is to provide important information for global health decision makers about the disease dynamics. The first control measure was lockdown and then health authorities imposed mask wearing and social distancing. Khajanchi et al. (2021) showed, by extending the classical SEIR model, that social distance is an important factor to reduce the reproduction number and this to reduce the virus spread. Despite the herd immunity acquired via vaccine or infection, the social distancing is still recommended as a public health measure. Driven by the observed characteristics of COVID-19, we propose a mathematical model with two infectious states. It was reported by World Health Organization that one in three people who get COVID-19 do not show any symptoms. This is a challenging problem for health authorities as the asymptotic individuals carry the virus and may infect other people without knowing it. Moreover, consequent efforts were made worldwide since the authorisation of new vaccines by the end of 2020. By the end of November 2021, more than 50% of the world population has first dose administered and only 40% has second dose administered. In order to study the efficacy of vaccination to contain the virus spread and its negative consequences, our model include vaccinated state. The objective is to provide efficient public health policies in determining optimal vaccination strategies. Some questions have raised since the beginning of vaccination campaigns: how many individuals should be vaccinated? Is the vaccine a solution to get rid of the disease permanently? These questions are related to financial and moral costs associated with the chosen governmental policy. This paper gives theoretical and numerical analysis associated with COVID-19 epidemic dynamics in order to answer these critical questions. Although we focus mainly on the Saudi case, the model structure is general and numerically adapted to any specific context without loss of validity of the qualitative results here shown. The main contributions of this research are given as follows:

• Developing a novel mathematical model to predict the spread of COVID-19, with the presence of vaccinated and asymptotic compartments

• Analyzing the existence of endemic equilibrium point and the stability of disease-free equilibrium

• Investigating a real case study in Saudi Arabia, discussing the impact of vaccination on disease dynamics

In ‘Related Works’, we present the related works regarding epidemic modeling with a focus on COVID-19 control strategies and particularly population vaccination. ‘Proposed Model and Effective Reproduction Number’ and ‘Model Analysis’ include model description and analysis, respectively. The numerical results are given in ‘Numerical Simulations’ and ‘Conclusion’ concludes this paper.

Related Works

The mathematical modeling in epidemiology started in England, in the 18th century, when Bernoulli analyzed the mortality of smallpox. Since then, a large variety of epidemiological models have been developed (Gumel et al., 2021; Mandal et al., 2020; Batistela et al., 2021; Samui, Mondal & Khajanchi, 2020). In this section, we present recent works proposed in this century, impacted by several outbreaks such as Ebola, Zika, and the swine flu (Bekiros & Kouloumpou, 2020). Alexander & Moghadas (2005) developed a Susceptible-Infected-Recovered-Susceptible (SIRS) epidemic model. The authors considered that the immunity acquired by the population after infection decreases over time. The dynamical behavior of the model is investigated using different types of bifurcation, including saddle–node, Hopf, and Bogdanov-Takens. The stability analysis based on the basic reproductive number and the rate of loss of natural immunity demonstrated the coexistence of two concentric limit cycles. These theoretical results have epidemiological implications such the determination of epidemic outbreak and the control the disease spread.

Yi, Zhao & Zhang (2009) investigated the Susceptible, Exposed, Infectious, Quarantine, Susceptible (SEIQS) epidemic model, with a nonlinear incidence rate. This model takes into consideration the communal sanitation measure of quarantine, aiming at avoiding broad infection. The authors provided a stability analysis using codimension-1 (transcritical, saddle–node, and Hopf) and codimension-2 bifurcations (Bogdanov-Takens).

Recently, Lu et al. (2019) studied the the SIRS epidemic model, the same considered in Alexander & Moghadas (2005) but with a generalized non-monotone incidence rate. The incidence rate is a function of the infection force of a disease and the number of susceptible individuals. The given formula for the incidence rate models the psychological pressure of some epidemic disease. The government is, in general, led to take some protective measures like lockdown when the infection number becomes very high. The authors showed that the model has both repelling and attracting Bogdanov-Takens bifurcations. Moreover, from the super-critical Hopf bifurcation, the authors concluded that a disease following this model presents periodic outbreak, which is very important to understand its dynamics, in the real world.

The impact of treatment function was investigated in Xiao et al. (2015) using the SIS model, where recovered individuals become again susceptible and the incidence rate is bi-linear. In the considered model, the treatment function is saturated, which results in the existence of backward bifurcation. Thus, the eradication of the disease is not only related to the reproduction number but also to other biological or epidemiological mechanisms, such as imperfect vaccine. The bifurcation analysis outlines the necessary conditions to eliminate the disease. Zhang, Ge & Lin (2020) discussed the impact of the number of hospital beds on SIS epidemic model, by considering a nonlinear recovery rate. The authors calculated the basic reproduction number corresponding to their model. This number determines the condition for the disease-free equilibrium to be globally asymptotically stable.

The limitations of medical resources, mainly the availability of vaccinations, is modeled using a piecewise-defined function for patient treatment in Wang, Xiao & Smith (2019). This function admits a backward bifurcation with limited available medical resources. The variation of vaccination threshold affects the existence of multiple steady states,crossing cycle, and generalized endemic equilibria. Similarly, Perez, Avila-Vales & Garcia-Almeida (2019) considered nonlinear incidence rate for a generalized SIR model. Besides, the authors assumed that the model has saturated Holling type II treatment rate and logistic growth. Non linear and saturated functions allows to represent more accurately the dynamics epidemic diseases. Similar to previous stated works, the authors revealed the importance of the basic reproduction number R₀, whose value determines the existence of endemic equilibrium and the stability of the disease-free equilibrium. Under some conditions related to the disease transmission rate and the treatment rate, the model may undertake a backward bifurcation and a Hopf bifurcation. The above-mentioned articles considered general disease models. In the literature, we can also find specific models targeting particular disease such as avian influenza (Jiang et al., 2020) and bacterial meningitis (Asamoah et al., 2020). Since the declaration of World Health Organization (WHO) of the Severe Acute Respiratory Syndrome Coronavirus (SARS-CoV-2) as a pandemic on March 2020, the scientific community has been trying to understand the dynamics of this virus.

One of the measures to control the virus spread in to declare total or partial lockdown, forcing social distancing. The scientific community believes that the main cause of infection is the inhalation of virus droplets (Jayaweera et al., 2020). Gevertz et al. (2021) modeled social distancing as a flow rate between susceptible and asymptomatic individuals. The model reveals the existence of of a critical implementation delay, when implementing social distancing mandates. A delay of two weeks is the critical threshold between infection containment and infection expansion.

Nadim & Chattopadhyay (2020) investigated the effect of imperfect lockdown. In the adopted model, when the basic reproduction number, R₀ is less than unity, the stable disease free equilibrium coexists with a stable endemic equilibrium. This means that COVID-19 undergoes backward bifurcation. This phenomenon was observed in the Kingdom of Saudi Arabia where the new cases were decreasing to reach 97 in 06, January 2021. Unfortunately, this rate reached 386 new cases, after one month, which obliged the Ministry of Health to declare partial lockdown for 10 days. The infection force is so high that the disease cannot be totally eradicated. The authors showed that under perfect lockdown, this backward bifurcation does not exist, but such condition is not possible in the real world. In Di Giamberardino & Iacoviello (2021), the authors included in their mathematical model, based on the classical SEIR, several prevention actions such as test campaign on the population and quarantining infected persons. The model took in consideration infection treatment efforts, such as vaccination and the therapy of induced cardio-respiratory complications. Besides the usual classes of the population, the authors considered two new classes, driven by specific characteristics of the virus: infected but asymptomatic patients and suspected infected individuals. The theoretical results, tuned using the Chinese case, were compared to United Kingdom case and the Italian case, showing the similarity between the model dynamics and the real epidemic behaviour. Lü et al. (2021) proposed an epidemic model that distinguishes between the first and the second waves od COVID-19 in China. The two-stage model includes a Contacts compartment, besides the usual Susceptible, Infectious and Recovered compartments. The authors of Khajanchi, Sarkar & Banerjee (2022) raised the issue of undetected symptomatic and asymptomatic individuals. They also investigated the effect of two control strategies that include the improvement of treatment for infected and isolated individuals. Another scientific aspect of COVID-19 is the possible transmission of the virus through contaminated surfaces. It is believed that the virus can survive several days on the surfaces depending on the material (wood, glass or plastic). Another issue faced by the governments is the awareness level of the population. Some individuals, deliberately, decide not to apply the precautionary measures, mainly wearing mask and respecting social distancing (Kassa, Njagarah & Terefe, 2020).

The issue of the efficiency of social distancing and rapid testing strategies against the pandemic was examined in Aldila et al. (2020), where the authors extended the standard SEIR model. The authors considered also the problem of undetected asymptomatic individuals, who have no symptoms but participate actively to virus spread. Furthermore, the limitation of medical resources was incorporated to the model. The theoretical findings emphasized the role of the basic reproduction number R₀ in the existence of stable COVID-19 free and COVID-endemic equilibrium points. This conclusion is contested by Mohd & Sulayman (2020), who studied the SIRS model with limited medical resources and false detection issues. The authors showed that the condition of reducing the basic reproduction number under the unity value is necessary to eliminate the disease but not sufficient. Since the authorization COVID-19 vaccines, several research works focused on giving insight to mathematical characteristic of virus spread after population vaccination. Algehyne & Ibrahim (2021) used nonlinear functional analysis and fractal fractional derivative to model the evolution over time of four compartments: susceptible, infected, infected positive tested, and recovered. The Spanish case was investigated in Kumar, Erturk & Murillo-Arcila (2021) using also fractional derivatives. It is important to highlight that these works do not consider vaccinated state as a separate compartment. They rather consider that the vaccinated individuals are moved from susceptible to recovered compartment. The vaccinated individuals are considered to move also from exposed state in Wintachai & Prathom (2021).

Different mathematical tools are used by Rajaei et al. (2021) to compare the effect of vaccination with social distancing and hospitalization. The extended Kalman filter (EKF) is used for state estimation under uncertainty.

Most of the existing research works developing a relationship between infectious and asymptotic individuals focus on estimating the model parameters using actual data (Asamoah et al., 2021; Gevertz et al., 2021; Sasmita et al., 2020; Giordano et al., 2020). To the best of our knowledge, our work is the first to provide to study mathematical stability of endemic and disease free equilibrium. Most relevant works in COVID-19 mathematical modeling are compared in Table 1, in terms of disease control strategy, considered model and country of the case study.

Proposed model and effective reproduction number

Our objective is to derive the mathematical equations that better present the dynamics of COVID-19 virus. The population is divided into five compartments: susceptible, vaccinated, infectious, asymptotic, and recovered; the numbers in these states are denoted by S(t), V(t), I(t), A(t), and R(t), respectively. Figure 1 depicted the flow diagram of the disease spread.

Table 2 and Table 3 summarize the different model parameters and variables.

All newborns are assumed to be susceptible. The natural recruitment and the natural death are denoted by Λ and µ, respectively. The disease-induced death rate is ignored. Susceptible individuals are vaccinated at rate constant ψ. The parameters α and β are the infecting rates of asymptotic and infectious individuals, respectively. γ₁ and γ₂ are the rates that the infectious and asymptotic individuals become recovered and acquire temporary immunity, respectively. The vaccinated individuals need a period of time to develop their immunity against the virus, represented by $\frac{1}{\eta }$.

The virus may infect vaccinated individuals but at a lower rate than susceptible individuals who are not unvaccinated. Thus in this case, the transmission rates β and α are multiplied by a scaling factor ɛ₁ and ɛ₂ (0 ≤ɛ₁, ɛ₂ ≤1).

Based on the above assumptions and Fig. 1, we formulate the following model of differential equations.

(1a)${\displaystyle \frac{dS}{dt}=\Lambda -\mu S-\psi S-\beta IS-\alpha AS}$ (1b)${\displaystyle \frac{dV}{dt}=-\mu V+\psi S-\eta V-{\varepsilon }_1\beta IV-{\varepsilon }_2\alpha AV}$ (1c)${\displaystyle \frac{dI}{dt}=-\mu I+\beta IS+{\varepsilon }_1\beta IV+\sigma A-{\gamma }_{1}I}$ (1d)${\displaystyle \frac{dA}{dt}=-\mu A+\alpha AS-\sigma A+{\varepsilon }_2\alpha AV-{\gamma }_{2}A}$ (1e)${\displaystyle \frac{dR}{dt}=-\mu R+\eta V+{\gamma }_{1}I+{\gamma }_{2}A}$

The basic reproduction number is defined as the number of secondary infections produced by a single infectious individual during his or her entire infectious period. Since we introduce a vaccination program in our model, it is called the effective reproduction number. The system (1) has always a disease-free equilibrium, which is obtained by setting all the derivatives to zero with I = A = 0, that yields to: $P_0=\left(S_0,I_0,A_0,R_0,V_0\right)=\left(\frac{\Lambda }{\mu +\psi },0,0,\frac{\eta \psi \Lambda }{\mu \left(\mu +\psi \right)\left(\mu +\eta \right)},\frac{\psi \Lambda }{\left(\mu +\psi \right)\left(\mu +\eta \right)}\right)$

Let x = (I, A, V, R, S)^T. System (1) can be rewritten as $x^{\prime }=\mathcal{F}\left(x\right)-\mathcal{N}\left(x\right),$, where $\mathcal{F}$ be the rate of appearance of new infections in each compartment. The progression from A to I is not considered to be new infection, but rather the progression of an infected individual through various infectious compartments.

$\begin{array}{c}\mathcal{F}\left(x\right)=\left(\begin{array}{c}\hfill \beta IS+{\varepsilon }_1\beta IV\hfill \\ \hfill \alpha AS+{\varepsilon }_2\alpha AV\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \\ \hfill 0\hfill \end{array}\right),\mathcal{N}\left(x\right)=\left(\begin{array}{c}\hfill \left(\mu +{\gamma }_1\right)I-\sigma A\hfill \\ \hfill \left(\mu +\sigma +{\gamma }_2\right)A\hfill \\ \hfill \left(\mu +\eta \right)V-\psi S+{\varepsilon }_1\beta IV+{\varepsilon }_2\alpha AV\hfill \\ \hfill \mu R-\eta V-{\gamma }_{1}I-{\gamma }_{2}A\hfill \\ \hfill -\Lambda +\left(\mu +\psi \right)S+\beta IS+\alpha AS\hfill \end{array}\right).\hfill \end{array}$

The infected compartments are A and I, giving m = 2. With A =I =0, the Jacobian matrices of $\mathcal{F}\left(x\right)$ and $\mathcal{N}\left(x\right)$ at the disease-free equilibrium P₀ are, respectively,

$\begin{array}{c}D\mathcal{F}\left(P_0\right)=\left(\begin{array}{cccc}\hfill F\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right),D\mathcal{N}\left(P_0\right)=\left(\begin{array}{cccc}\hfill N\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\varepsilon }_1\beta V_0{\varepsilon }_2\alpha V_0\hfill & \hfill \mu +\eta \hfill & \hfill 0\hfill & \hfill -\psi \hfill \\ \hfill -{\gamma }_1-{\gamma }_2\hfill & \hfill -\eta \hfill & \hfill \mu \hfill & \hfill 0\hfill \\ \hfill \beta S_0\alpha S_0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \mu +\psi \hfill \end{array}\right),\hfill \end{array}$

$\begin{array}{c}F=\left(\begin{array}{cc}\hfill \beta S_0+{\varepsilon }_1\beta V_0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \alpha S_0+{\varepsilon }_2\alpha V_0\hfill \end{array}\right),N=\left(\begin{array}{cc}\hfill \mu +{\gamma }_1\hfill & \hfill -\sigma \hfill \\ \hfill 0\hfill & \hfill \mu +\sigma +{\gamma }_2\hfill \end{array}\right).\hfill \end{array}$

Our developed model is similar to the two-strain model in Van den Driessche & Watmough (2002) with two infectious compartments. FN⁻¹, the next generation matrix of system (1) has the two eigenvalues.

$R_1=\frac{\beta \left(S_0+{\varepsilon }_1{V}_0\right)}{\mu +{\gamma }_1}=\frac{\beta \Lambda \left(\frac{1}{\mu +\psi }+{\varepsilon }_1\frac{\psi }{\left(\mu +\psi \right)\left(\mu +\eta \right)}\right)}{\mu +{\gamma }_1}=\frac{\beta \Lambda \left(1+{\varepsilon }_1\frac{\psi }{\mu +\eta }\right)}{\left(\mu +{\gamma }_1\right)\left(\mu +\psi \right)}$

$R_2=\frac{\alpha \left(S_0+{\varepsilon }_2{V}_0\right)}{\mu +{\gamma }_2+\sigma }=\frac{\alpha \Lambda \left(1+{\varepsilon }_2\frac{\psi }{\mu +\eta }\right)}{\left(\mu +{\gamma }_2+\sigma \right)\left(\mu +\psi \right)}=\frac{\alpha \Lambda \left(\mu +\eta +{\varepsilon }_2\psi \right)}{\left(\mu +{\gamma }_2+\sigma \right)\left(\mu +\psi \right)\left(\mu +\eta \right)}$

The effective reproduction number for the system is the maximum of the two.

Model analysis

Numerical simulations

In this paper, we focus on vaccination analysis in Saudi Arabia. The presented numerical simulations provide also general results that can be applied to any region. The data set is provided by King Abdullah Petroleum Studies and Research Center (KAPSARC). It includes five classes: Tested, Cases, Recoveries, Critical, Mortalities and Active and it spans the period from 04/03/2020 to 08/11/2021. It includes also important events and measures such as international flights suspension and lockdown. We used the Simulink Tool in order to simulate different scenarios.

The death and birth rate for Saudi Arabia are estimated to be equal to 3.39 and 14.56 for 1,000 per year, respectively. The vaccination campaign started on 18/12/2020 with a vaccination coverage of the total population of 0.02% to reach about 65% of the adult population fully vaccinated in November 2021. The vaccination rate is considered a the percentage of the total population that get vaccinated per day. With approximately 45,000 administrated doses per day and a total population of 35339000 in 2021, this rate is about 0.00 127. Is this rate enough to eradicate the disease? This what we are trying to answer is this work.

The estimation of model parameters is very important to have accurate numerical results. One approach consists in calibrating the model by fitting it with reported data using, for example, least square method (Rai et al., 2022). In our paper, we used scientific reported facts about COVID-19 transmission mechanisms. The research report (Nogrady, 2020) provides information about asymptotic individuals for COVID-19. Most people, with no symptoms at the beginning, develop symptoms in 7–13 days, which corresponds to the σ⁻¹. Recall that γ₁ is the recovery rate of infectious individuals. Interpreted as the expected value of a Poisson process, γ₁ can be related to the needed time from the beginning of the infection till recovery (Gevertz et al., 2021). With average recovery duration equal to 10 days (Chae et al., 2020), the recovery rate of infectious individuals is γ₁ = 0.1.

Let ω denote the fraction of asymptotic individuals among positive cases. According to Nogrady (2020), and based on 13 studies involving 21,708 people in 2020, ω = 0.17. Using the same methodology as in Gevertz et al. (2021) ${\gamma }_2=\frac{\omega }{1-\omega }\sigma \approx 0.2\sigma$. The asymptotic people are estimated to be 42% less contagious than symptomatic individuals (Nogrady, 2020). Thus, α = 0.42β. Table 4 summarizes selected values for the model parameters.

By 18/12/2020, considered as time 0 in this model, the number of recoveries is equal to 351722, the number of active cases is equal to 3,014. Assuming that ω = 0.17, the number of initial infectious with and without symptoms is equal to 2,501 and 513, respectively.

First, we investigate the effect of the vaccination rate on the effective reproduction number defined as the maximum of the two entities R₁ and R₂. According to Chae et al. (2020), COVID-19 transmission rate β ranges between 0.233 and 0.462. Figures 2A and 2B show the evolution of both R₁ and R₂ as a function of the vaccination rate ψ with virus transmission rate equal to 0.233 and 0.462, respectively. In both cases, R₁ corresponding to the strain of infectious individuals with symptoms is greater than R₂, corresponding to the strain of asymptotic individuals. Thus the number of individuals infected by one person currying the virus is mainly affected by individuals showing usual symptoms. Mathematical theoretical result confirms that the vaccine reduces the spread of the virus among the population. We would like to highlight the fact in our model that a vaccinated individual may be infectious with or without symptoms. This result is very important as, till the end of 2021, an important portion of worldwide population is still opposed to vaccine. In the case of high transmission rate and low vaccination rate, R₁ is higher than 1. The disease free equilibrium is consequently unstable according to theorem 1. For the Saudi case, when beta is equal to 0.233, R₁ and R₂ are equal to 0.0797 and 0.0195, respectively. When beta is equal to 0.462, R₁ and R₂ are equal to 0.1580 and 0.0387, respectively. For Saudi Arabia, the effective reproduction number is less than 1, even for high transmission rate. This result is explained by the high vaccination rate.

Figure 3 shows the weekly number of new active cases and recovered after infection in KSA, starting from 18/12/2020, the date when the vaccination starts. We can see that after 12 weeks, the number of cases raises. This behaviour was surprising for a population waiting to see the effect of vaccination. It’s only after 31 weeks that the number of new cases start to decrease. The same phenomena was observed in both Figs. 4C and 5C. The theatrical results is conform to actual statistics. The effect of vaccination is not immediate; it needs several weeks to observe a decrease in the number of new infectious cases.

We compare the evolution in time of the five compartments (S, V, I, A, and R) presented in our model, for two different transmission rate and with the Saudi vaccination rate. With different dynamics at the beginning, both scenarios show a convergence to a stable state. We observe almost similar patterns for S, V. I and R. The number of susceptible individuals S decreases slowly at the beginning and then, we observe a drastic decline. Obviously, the number of recovered follows the same slow and then fast pace but in decrease. The number of vaccinated individuals V increases gradually at the beginning and then it begins to fall down. The number of infectious individuals I remain stable for a short period to witness an expansion followed by a decline. The number of asymptotic individuals show different evolution patterns for two considered scenarios. When we set a low value for the virus transmission rate, this number immediately shrinks. However, when we set a high value for the virus transmission rate, this number increases before shrinking.

The effect of the transmission rate can also be observed in the amplitude of each category of individuals. When the models converge, the number of infectious and asymptotic individuals are zero. We emphasize here our theoretical result, mentioned in Theorem 1, that states that if both R₁ and R₂ are less than one, the disease free equilibrium is stable. This result is consistent with the simulation results. The difference between the two considered scenarios lies in the percentage of susceptible and vaccinated individuals in the equilibrium. This percentage is very low when the transmission rate is high. Although the percentages of vaccinated individuals are close, we observe a remarkable difference in the number of infectious individuals. When the transmission rate is high, almost 40% of the population is infected, which rises public health issues.

Conclusion

In this article, we present a mathematical model for COVID-19, based on the virus behaviour. Our main target is to evaluate the effect of vaccination on the population. The presence of individuals presenting no symptoms and the immunity loss are the main characteristics that make this virus different from other known and already modeled diseases. We provide analytical expression of the effective reproduction number with is a key factor to determine necessary conditions for endemic and disease free equilibrium. We supported our theoretical findings with the numerical analysis applied to the Saudi case. The main findings of this work are as follows.

• People can observe that the vaccine helps to reduce the severity of symptoms. We gave a mathematical proof that the vaccine reduces the transmission rate.

• The vaccination campaign in Saudi Arabia was immediately followed by the rise in the number of infections. We have showed that this observation is mathematically justified and this rise is a necessary transition before the increase of new infections.

• By adjusting the model parameters based on collected data, we provide the decision-makers with the vaccination rate necessary for virus spread control.

Recently, the scientific community is observing the new variants that show each time different patterns. We aim, in the future, to develop a new model with the new observed characteristics of variants such as beta and omicron. As discussed in the previous section, several countries including Saudi Arabia did not face medical resource shortages. However, as highlighted by Lotfi et al., (2021), the management of medical waste is critical during COVID-19 pandemic. Moreover, the number of asymptomatic and untested symptomatic infections is uncertain. As future work, we propose to capture disease dynamic uncertainty and incorporate risk assessment, to alleviate the impacts of pandemic peak. Hybrid fuzzy, data-driven, robust optimization, and stochastics (Lotfi et al., 2022a; Lotfi et al., 2022b) are examples of uncertainty analysis methods. Conditional value at risk (CVaR), which is the average shortfall beyond the VaR point, is a consistent and coherent risk assessment measure.

Supplemental Information