An Extended SEIR Model with Vaccination for Forecasting the COVID-19 Pandemic in Saudi Arabia Using an Ensemble Kalman Filter

Abstract

In this paper, an extended SEIR model with a vaccination compartment is proposed to simulate the novel coronavirus disease (COVID-19) spread in Saudi Arabia. The model considers seven stages of infection: susceptible (S), exposed (E), infectious (I), quarantined (Q), recovered (R), deaths (D), and vaccinated (V). Initially, a mathematical analysis is carried out to illustrate the non-negativity, boundedness, epidemic equilibrium, existence, and uniqueness of the endemic equilibrium, and the basic reproduction number of the proposed model. Such numerical models can be, however, subject to various sources of uncertainties, due to an imperfect description of the biological processes governing the disease spread, which may strongly limit their forecasting skills. A data assimilation method, mainly, the ensemble Kalman filter (EnKF), is then used to constrain the model outputs and its parameters with available data. We conduct joint state-parameters estimation experiments assimilating daily data into the proposed model using the EnKF in order to enhance the model’s forecasting skills. Starting from the estimated set of model parameters, we then conduct short-term predictions in order to assess the predicability range of the model. We apply the proposed assimilation system on real data sets from Saudi Arabia. The numerical results demonstrate the capability of the proposed model in achieving accurate prediction of the epidemic development up to two-week time scales. Finally, we investigate the effect of vaccination on the spread of the pandemic.

1. Introduction

The novel coronavirus disease (COVID-19) first appeared in Wuhan city, and due to its high transmission rate, the virus has been rapidly spreading all over the world. On 8 March 2020, the WHO announced COVID-19 as a global pandemic [1]. As of 1 January 2021, there are over 100 million reported cases, and a death toll exceeding 2 million.

In the absence of a ready-to-use vaccine, and besides medical and biological research, mathematical models can play an important role in understanding and predicting disease transmission. Moreover, it helps to implement appropriate measures and efficient strategies to control the pandemic’s spread and mitigate its impact.

Mathematical epidemic models based on the SIR model [2] are widely used to study the spread of a disease in a population [3,4,5,6,7]. Recently, several mathematical models have been developed to study the transmission dynamics of COVID-19 [8,9,10,11]. However, these models suffer from various sources of uncertainties, due to the incomplete description of the biological processes governing the disease spread, and also due to some involved parameters being poorly known. One way to mitigate these uncertainties is to constrain epidemic forecasting models with available data. These data can be combined with the model output to improve its prediction and reduce uncertainties [12,13,14]. This process, known as data assimilation, sequentially adjusts the model state and parameters once data become available to keep the model as close as possible to the real track [15].

The ensemble Kalman filter (EnKF) is a popular data assimilation technique introduced by [16] to tackle large-scale, sequential, and non-linear problems. The EnKF is a Monte Carlo implementation of the traditional Kalman filter [17], operating in cycles of two steps: (i) A forecast step in which a set of state samples, or ensemble members, are integrated forward in time with the model, and (ii) an analysis step to update the forecasted members with available data based on the Bayes rule in a Gaussian framework [15,16]. Because of its affordable computational cost, robust performance, nonintrusive formulation, and ease of implementation, the EnKF have been successfully implemented in many geophysical applications, including weather forecasting [18,19,20], hydrology and reservoir applications [21,22,23,24,25,26], ocean modeling [27,28,29,30,31], and flood modeling [32,33,34].

In epidemiology, previous studies have used data assimilation in the context of both variational and Kalman filter methods [35,36,37,38]. The Reference [39] used the EnKF to estimate the parameters of a stochastic SEIR model. The Reference [40] applied a deterministic variant of the EnKF, the Ensemble Adjustment Kalman Filter (EAKF), in combination with a network of SEIR models, simulating different connected cities. The Reference [41] used an ensemble smoother with a multiple data assimilation (ESMDA) approach for parameter estimation in a SEIR model. The Reference [42] also implemented the ESMDA to effectively estimate the parameters of a more sophisticated SEIR model accounting for age-classes and compartments of the sick, hospitalized, and dead. The Reference [43] used an EnKF approach for parameter estimation and shot-term forecasts with a simple SIR model. The Reference [44] proposed an EnKF for estimating unmeasurable state variables and unknown parameters of a SIRV model, which takes into account the circulation of free coronavirus in the environment. In the context of variational data assimilation methods, the Reference [45] performed parameter estimation and predictions using a SIR model, while [46] proposed a modified SEIR model that distinguishes between symptomatic and asymptomatic, and conducted observation sensitivity experiments to identify suitable observing strategies. However, these studies did not consider the impact of vaccination on the pandemic spread.

This work therefore aims to study the impact of vaccination on the COVID-19 spread. An extended SEIR model comprising of seven compartments—susceptible, exposed, infectious, quarantined, recovered, deaths, and vaccinated—is then proposed. First, we conduct a mathematical analysis to illustrate the non-negativity, boundedness, epidemic equilibrium, existence, and uniqueness of the endemic equilibrium, and the basic reproduction number of the model. A Joint-EnKF is implemented for state-parameters estimation by assimilating COVID-19 data into the model. By considering uncertainties in all model states and poorly known parameters, our ensemble data assimilation approach computes updates of both state and parameter ensembles. The estimated set of model state and parameters is then used to assess the model prediction skill. The proposed approach is validated with the cases reported from the most-infected Gulf country, Saudi Arabia, showing a great ability to provide accurate and consistent forecasts up to two weeks. We finally investigate the effect of vaccination on the spread of the pandemic, showing a notable decrease in the number of cases.

The remainder of this paper is organized as follows. Section 2 presents the extended SEIR model and briefly describes the Joint-EnKF algorithm. Section 3 is devoted to the mathematical analysis of the model. The numerical results are presented and discussed in Section 4. Section 5 concludes the study, with a summary of the main results.

2. Model and Assimilation Scheme

2.1. Extended SEIR Model

We extend the SEIR model to seven compartments to simulate the epidemic of COVID-19. Seven state variables are considered within a population, that is, $S\left(t\right)$, $E\left(t\right)$, $I\left(t\right)$, $Q\left(t\right)$, $R\left(t\right)$, $D\left(t\right)$, and $V\left(t\right)$, denoting the number of susceptible, exposed (infected, but not yet infectious), infectious (not yet quarantined), quarantined (confirmed and infected), recovered, dead, and vaccinated cases, respectively. The disease transmission flow of the proposed model is sketched in Figure 1. The model is then governed by the following set of nonlinear ordinary differential equations:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS\left(t\right)}{dt}}& =\Lambda -\beta S\left(t\right)I\left(t\right)-\alpha S\left(t\right)-\mu S\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dE\left(t\right)}{dt}}& =\beta S\left(t\right)I\left(t\right)-\gamma E\left(t\right)+\sigma \beta V\left(t\right)I\left(t\right)-\mu E\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dI\left(t\right)}{dt}}& =\gamma E\left(t\right)-\delta I\left(t\right)-\mu I\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dQ\left(t\right)}{dt}}& =\delta I\left(t\right)-(1-\kappa )\lambda Q\left(t\right)-\kappa \rho Q\left(t\right)-\mu Q\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dR\left(t\right)}{dt}}& =(1-\kappa )\lambda Q\left(t\right)-\mu R\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dD\left(t\right)}{dt}}& =\kappa \rho Q\left(t\right),\hfill \\ \hfill {\displaystyle \frac{dV\left(t\right)}{dt}}& =\alpha S\left(t\right)-\sigma \beta V\left(t\right)I\left(t\right)-\mu V\left(t\right),\hfill \end{array}$$

(1)

with non-negative initial conditions $S\left(0\right)=S_0$, $E\left(0\right)=E_0$, $I\left(0\right)=I_0$, $Q\left(0\right)=Q_0$, $R\left(0\right)=R_0$, $D\left(0\right)=D_0$, and $V\left(0\right)=V_0$. The coefficients $\Lambda$, $\beta$, $\alpha$, $\mu$, $\gamma$, $\sigma$, $\delta$, $\kappa$, $\lambda$, and $\rho$ represent the new births and new residents per unit of time, transmission rate divided by N, vaccination rate (rate of people who are vaccinated), natural death rate, average latent time, average quarantine time, mortality rate, average days until recovery, and average days until death, respectively. $\sigma$ is the vaccine inefficacy $(0\le \sigma \le 1)$. So $1-\sigma$ represents the vaccine efficacy. If $\sigma =0$, the vaccine offers $100\%$ protection against the disease. $N=S+E+I+Q+R+D+V$ is the total population size. Finally, system (1) is solved using a fourth-order Runge–Kutta method.

2.2. Ensemble Kalman Filter (EnKF)

Consider the following nonlinear data assimilation problem:

$$\begin{array}{cc}\hfill {\mathbf{x}}_{k+1}& =\mathcal{M}\left({\mathbf{x}}_k,\theta \right)+{\eta }_{k+1},\hfill \\ \hfill {\mathbf{y}}_{k+1}& =\mathcal{H}\left({\mathbf{x}}_{k+1}\right)+{ϵ}_{k+1},\hfill \end{array}$$

(2)

where ${\mathbf{x}}_k$ is the state vector of the model at time $t_k$, $\mathcal{M}$ is the dynamical operator (system (1)) that integrates the model state from time $t_k$ to time $t_{k+1}$, $\theta$ is the parameters vector, ${\mathbf{y}}_{k+1}$ is the observation vector (active, recovered, and death cases) at time $t_{k+1}$, $\mathcal{H}$ is the observation operator that projects the state space onto the observation space, and ${\eta }_{k+1}$ and ${ϵ}_{k+1}$ are model and observation noises assumed to follow Gaussian distributions with zero mean and covariances ${\mathbf{Q}}_{k+1}$ and ${\mathbf{R}}_{k+1}$, respectively.

At every forecast step, we integrate the ensemble members with the dynamical model (1). The predicted state mean and the associated error covariance matrix are then estimated

$${\widehat{\mathbf{x}}}_{k+1}^f={\displaystyle \frac{1}{N_e}}\sum _{i=1}^{N_e}{\mathbf{x}}_{k+1}^{f,i},$$

(3)

$${\widehat{\mathbf{P}}}_{k+1}^f={\displaystyle \frac{1}{N_e-1}}\sum _{i=1}^{N_e}\left({\mathbf{x}}_{k+1}^{f,i}-{\widehat{\mathbf{x}}}_{k+1}^f\right){\left({\mathbf{x}}_{k+1}^{f,i}-{\widehat{\mathbf{x}}}_{k+1}^f\right)}^T={\mathbf{X}}_{k+1}^f{\left({\mathbf{X}}_{k+1}^f\right)}^T,$$

(4)

where $N_e$ is the number of ensemble members.

The Joint-EnKF consists of estimating the augmented state-parameters vector ${\mathbf{z}}_{k+1}=\left[{\mathbf{x}}_{k+1}^T,{\theta }^T\right]$ at each time $t_{k+1}$ given the observation ${\mathbf{y}}_{k+1}$. The prediction step of the Joint-EnKF becomes:

$${\mathbf{z}}_{k+1}^f=\tilde{\mathcal{M}}\left({\mathbf{z}}_k,\theta \right)+{\eta }_{k+1}=\left[\begin{array}{c}\mathcal{M}\left({\mathbf{x}}_k,\theta \right)+{\eta }_{k+1}\\ \theta \end{array}\right],$$

(5)

where $\tilde{\mathcal{M}}$ is the augmented nonlinear dynamic operator for the joint state-parameters. With similar computations as in (4), the joint state-parameters error covariance matrix for the augmented ensemble vectors ${\mathbf{z}}_{k+1}^{f,i}$ is given by

$${\widehat{\mathbf{P}}}_{k+1}^f=\left[\begin{array}{c}{\widehat{\mathbf{P}}}_{{\mathbf{x}}_{k+1}{\mathbf{x}}_{k+1}}^f{\widehat{\mathbf{P}}}_{{\mathbf{x}}_{k+1}{\theta }_{k+1}}^f\\ {\widehat{\mathbf{P}}}_{{\theta }_{k+1}{\mathbf{x}}_{k+1}}^f{\widehat{\mathbf{P}}}_{{\theta }_{k+1}{\theta }_{k+1}}^f\end{array}\right]=\left[\begin{array}{c}{\mathbf{X}}_{{\mathbf{x}}_{k+1}}^f{\left({\mathbf{X}}_{{\mathbf{x}}_{k+1}}^f\right)}^T{\mathbf{X}}_{{\mathbf{x}}_{k+1}}^f{\left({\mathbf{X}}_{{\theta }_{k+1}}^f\right)}^T\\ {\mathbf{X}}_{{\theta }_{k+1}}^f{\left({\mathbf{X}}_{{\mathbf{x}}_{k+1}}^f\right)}^T{\mathbf{X}}_{{\theta }_{k+1}}^f{\left({\mathbf{X}}_{{\theta }_{k+1}}^f\right)}^T\end{array}\right],$$

(6)

where ${\mathbf{X}}_{\theta }$ denotes the parameter ensemble perturbations. The updated joint state and parameters are then computed using the current observation ensemble ${\mathbf{y}}_{k+1}^i$:

$${\mathbf{z}}_{k+1}^{a,i}={\mathbf{z}}_{k+1}^{f,i}+{\widehat{\mathbf{P}}}_{k+1}^f{\tilde{\mathcal{H}}}^T{\left(\tilde{\mathcal{H}}{\widehat{\mathbf{P}}}_{k+1}^f{\tilde{\mathcal{H}}}^T+{\mathbf{R}}_{k+1}\right)}^{-1}\left[{\mathbf{y}}_{k+1}^i-\tilde{\mathcal{H}}{\mathbf{z}}_{k+1}^{f,i}\right],$$

(7)

where $\tilde{\mathcal{H}}=\left[\mathcal{H}\mathbf{0}\right]$ is the augmented observational operator with the zero matrix.

3. Mathematical Analysis of the COVID-19 Model

In this section, we discuss the boundedness, non-negativity, epidemic equilibrium, endemic equilibrium, and basic reproduction number of the COVID-19 model (1).

3.1. Non-Negativity of the Model

Theorem 1. 

If $S_0\ge 0$, $E_0\ge 0$, $I_0\ge 0$, $Q_0\ge 0$, $R_0\ge 0$, $D_0\ge 0$, and $V_0\ge 0$, then the solutions of system (1) remain non-negative for all $t>0$.

Proof of Theorem 1.

From the first equation of system (1), we have:

$${\displaystyle \frac{dS}{dt}}\ge -\mu S.$$

(8)

Integrating Equation (8), we get

$$S\left(t\right)\ge S_0{e}^{-\mu t}\ge 0.$$

(9)

Hence, $S\left(t\right)$ remains non-negative for all $t>0$. Similarly, using the remaining equations of system (1), we obtain $E\left(t\right)\ge 0$, $I\left(t\right)\ge 0$, $Q\left(t\right)\ge 0$, $R\left(t\right)\ge 0$, $D\left(t\right)\ge 0$, and $V\left(t\right)\ge 0$. Therefore, Theorem 1 is proved. □

3.2. Boundedness of the Model

Theorem 2. 

All solutions of the proposed model with non-negative initial conditions are bounded and $N\left(t\right)\le {\displaystyle \frac{\Lambda }{\mu }}$ for all $t>0$.

Proof of Theorem 2.

Adding all the equations in the system (1), the growth rate of the total population can be written as

$${\displaystyle \frac{dN}{dt}}={\displaystyle \frac{dS}{dt}}+{\displaystyle \frac{dE}{dt}}+{\displaystyle \frac{dI}{dt}}+{\displaystyle \frac{dQ}{dt}}+{\displaystyle \frac{dR}{dt}}+{\displaystyle \frac{dD}{dt}}+{\displaystyle \frac{dV}{dt}}.$$

(10)

From Equation (10), we have the following equation:

$${\displaystyle \frac{dN}{dt}}=\Lambda -\mu N.$$

(11)

Integrating both sides of Equation (11) gives the solution for $N\left(t\right)$ as follows

$$N\left(t\right)={\displaystyle \frac{\Lambda }{\mu }}+\left(N_0-{\displaystyle \frac{\Lambda }{\mu }}\right)e^{-\mu t}.$$

(12)

Thus, when $t\to \infty$, we find that

$$N\left(t\right)\le {\displaystyle \frac{\Lambda }{\mu }}.$$

(13)

Finally, Theorem 1 and Equation (13) give

$$0\le N\left(t\right)\le {\displaystyle \frac{\Lambda }{\mu }}.$$

(14)

Thus, $S\left(t\right)$, $E\left(t\right)$, $I\left(t\right)$, $Q\left(t\right)$, $R\left(t\right)$, $D\left(t\right)$, and $V\left(t\right)$ are bounded and Theorem 2 is proved. □

3.3. Epidemic Equilibrium and Basic Reproduction Number of the Model

The epidemic equilibrium $X^0$ of system (1) is obtained by setting all the derivatives to zero with $I=0$, that yields to:

$$X^0=\left({\displaystyle \frac{\Lambda }{\mu +\alpha }},0,0,0,0,0,{\displaystyle \frac{\alpha \Lambda }{\mu (\mu +\alpha )}}\right).$$

(15)

$R_0$ is computed using the next-generation matrix concept [47]. Let $X={(E,I)}^T$. System (1) can be written as

$$X^{\prime }=\mathbb{F}\left(X\right)-\mathbb{W}\left(X\right),$$

(16)

where $\mathbb{F}\left(X\right)={\left(\beta SI+\sigma \beta IV,0\right)}^T$ and $\mathbb{W}\left(X\right)={\left((\mu +\delta )E,-\gamma E+(\mu +\delta )\right)}^T$. Their corresponding Jacobian matrices at the disease-free equilibrium are

$$F=\left(\begin{array}{cc}0& \beta S^0+\sigma \beta V^0\\ 0& 0\end{array}\right),W=\left(\begin{array}{cc}\mu +\gamma & 0\\ -\gamma & \mu +\delta \end{array}\right).$$

(17)

According to [47], the basic reproduction number $R_0$ is the spectral radius of the next-generation matrix $F{W}^{-1}$. That gives

$$R_0={\displaystyle \frac{\beta \gamma }{\left(\mu +\gamma \right)\left(\mu +\delta \right)}}\left(S^0+\sigma V^0\right).$$

(18)

Substituting $S^0$ and $V^0$ from Equation (15) into Equation (20) yields the following formula for the basic reproduction number of system (1):

$$R_0={\displaystyle \frac{\beta \gamma \Lambda (\mu +\alpha \sigma )}{\mu \left(\mu +\gamma \right)\left(\mu +\delta \right)(\mu +\alpha )}},$$

(19)

used to measure the transmission potential of a disease. It is the number of secondary cases or new infections produced by a single infectious individual in a completely susceptible population [47,48]. An epidemic is expected to increase exponentially if $R_0>1$, and to end if $R_0<1$ [47].

Theorem 3. 

The disease-free equilibrium $X^0$ is locally asymptotically stable if $R_0<1$, and unstable if $R_0>1$.

Proof of Theorem 3.

System (1) can be written in a simplified form by eliminating the sixth equation for the death compartment (D). This is due to the fact that the state D only appears in the corresponding differential equation and has no significance on the overall system. The number of deaths can be determined from $D=N-S-E-I-Q-R-V$. The Jacobian matrix of the reduced model (1) at the disease-free equilibrium is given by

$$J\left(X^0\right)=\left(\begin{array}{cccccc}-{ϵ}_1& 0& -\beta S^0& 0& 0& 0\\ 0& -{ϵ}_2& \beta \left(S^0+\sigma V^0\right)& 0& 0& 0\\ 0& \gamma & -{ϵ}_3& 0& 0& 0\\ 0& 0& \delta & -{ϵ}_4& 0& 0\\ 0& 0& 0& \lambda \left(1-\kappa \right)& -\mu & 0\\ \alpha & 0& -\sigma \beta V^0& 0& 0& -\mu \end{array}\right),$$

(20)

where ${ϵ}_1=-(\alpha +\mu )$, ${ϵ}_2=-(\gamma +\mu )$, ${ϵ}_3=-(\delta +\mu )$, and ${ϵ}_4=-(\lambda (1-\kappa )+\kappa \rho +\mu )$. $J\left(X^0\right)$ has ${\lambda }_1=-{ϵ}_1<0$, ${\lambda }_2=-{ϵ}_4<0$, ${\lambda }_3={\lambda }_4=-\mu <0$. The remaining two eigenvalues are

$$\begin{array}{cc}\hfill {\lambda }_5& =-{\displaystyle \frac{1}{2}}\left({ϵ}_2+{ϵ}_3+\sqrt{{\left({ϵ}_2-{ϵ}_3\right)}^2+4{ϵ}_2{ϵ}_3{R}_0}\right),\hfill \\ \hfill {\lambda }_6& =-{\displaystyle \frac{1}{2}}\left({ϵ}_2+{ϵ}_3-\sqrt{{\left({ϵ}_2-{ϵ}_3\right)}^2+4{ϵ}_2{ϵ}_3{R}_0}\right).\hfill \end{array}$$

(21)

when $R_0<1$, we have ${\lambda }_5<0$ and ${\lambda }_6<0$. Since all the eigenvalues are negative, the disease-free equilibrium $X^0$ is locally asymptotically stable. If $R_0>1$, we have ${\lambda }_5<0$ and ${\lambda }_6>0$, and the disease-free equilibrium $X^0$ is locally asymptotically unstable [49]. □

3.4. Existence and Uniqueness of the Endemic Equilibrium

In this section, we investigate the existence and uniqueness of the endemic equilibrium when $R_0>1$. Let $X^{*}=(S^{*},E^{*},I^{*},Q^{*},D^{*},V^{*})$ be the endemic equilibrium of model (1). The endemic equilibrium is obtained by setting all the derivatives in model (1) to zero, giving

$$\begin{array}{c}\hfill \Lambda -\beta S^{*}{I}^{*}-\alpha S^{*}-\mu S^{*}=0,\\ \hfill \beta S^{*}{I}^{*}-\gamma E^{*}+\sigma \beta V^{*}{I}^{*}-\mu E^{*}=0,\\ \hfill \gamma E^{*}-\delta I^{*}-\mu I^{*}=0,\\ \hfill \delta I^{*}-(1-\kappa )\lambda Q^{*}-\kappa \rho Q^{*}-\mu Q^{*}=0,\\ \hfill (1-\kappa )\lambda Q^{*}-\mu R^{*}=0,\\ \hfill \alpha S^{*}-\sigma \beta V^{*}{I}^{*}-\mu V^{*}=0.\end{array}$$

(22)

Solving the first, third, and fifth equations of (22) in terms of $I^{*}$ gives

$$E^{*}={\displaystyle \frac{{ϵ}_3}{\gamma }}{I}^{*},Q^{*}={\displaystyle \frac{\delta }{{ϵ}_4}}{I}^{*},S^{*}={\displaystyle \frac{\Lambda }{\beta I^{*}+{ϵ}_1}},$$

(23)

where ${ϵ}_1=\mu +\alpha$, ${ϵ}_3=\mu +\delta$, and ${ϵ}_4=\mu +\lambda (1-\kappa )+\kappa \rho$. Adding the second and sixth equation in (22) gives

$$V^{*}={\displaystyle \frac{\Lambda \beta \gamma I+\Lambda \alpha \gamma -{ϵ}_2{ϵ}_3\left(\beta I+{ϵ}_1\right)I}{\mu \gamma \left(\beta I+{ϵ}_1\right)}}.$$

(24)

Now substituting the expressions of $S^{*}$, $E^{*}$, and $V^{*}$ in the second equation of (22), we get the following quadratic equation for $I^{*}$:

$$a_2{I^{*}}^2+a_1{I}^{*}+a_0=0,,$$

(25)

where

$$\begin{array}{c}\hfill a_0=\mu {ϵ}_1{ϵ}_2{ϵ}_3\left(1-R_0\right),\\ \hfill a_1={ϵ}_1{ϵ}_2{ϵ}_3+\mu \beta {ϵ}_2{ϵ}_3-\Lambda {\beta }^2\gamma \sigma ,\\ \hfill a_2=\beta {ϵ}_2{ϵ}_3.\end{array}$$

(26)

It is clear that $a_2$ is always positive and $a_0$ is negative when $R_0>1$. Therefore, there is a unique positive value for $I^{*}$, and consequently, a unique endemic equilibrium $X^{*}$ when $R_0>1$.

4. Results and Discussion

The proposed modeling system was applied to forecast the COVID-19 pandemic in Saudi Arabia. Data are available from the Saudi Center for Diseases Prevention and Control [50]. These include the number of deaths, as well as recovered and confirmed cases. The total number of deaths, active cases, and recovered cases were assimilated daily, and model parameters were also updated using the Joint-EnKF.

4.1. Assimilation Settings

In this section, we cover the period before the vaccine becomes available. Therefore, state and parameters related to vaccination are excluded from the model at this stage. The state vector x to be estimated is then composed of six elements

$$x=\left[SEIQRD\right],$$

(27)

and the parameters to be estimated become

$$\theta =\left[\Lambda \beta \mu \gamma \delta \kappa \lambda \rho \right].$$

(28)

The augmented state vector to be estimated is then:

$$z=\left[x^{T}log\left({\theta }^T\right)\right],$$

(29)

where the parameters are in a logarithmic sense in order to avoid negative values in the model.

We set the initial numbers of infected and exposed people to $I\left(0\right)=E\left(0\right)=3$, three times the number of accumulated positive tests $\left(Q_0=1\right)$ provided by March 2nd [50]. The initial numbers of recovered, and death cases are $R\left(0\right)=D\left(0\right)=0$. The total population and the total number of susceptible cases are $N\left(0\right)=$ 34,218,000 and $S\left(0\right)=N\left(0\right)-E\left(0\right)-I\left(0\right)-Q\left(0\right)-R\left(0\right)-D\left(0\right)$, respectively. We consider the following two periods: the first period concerns the period of the beginning of the pandemic when the intervention had not been applied using a basic reproduction number greater than one, $R_0=2.5$, which gives ${\beta }_1=8.58\times {10}^{-9}$, while the second period corresponds to the period when the intervention was applied by the government until 30 June. During this period, the numbers of new cases are relatively constant. Therefore, we use a basic reproduction number $R_0=1$, which gives ${\beta }_2=3.43\times {10}^{-9}$. To reflect the uncertainty in the transmission rate, we use a standard deviation of $25\%$. The initial values of the remaining parameters are listed in

Table 1. Initial model parameters and biological interpretations.

Parameter	Initial Value	Description	Reference
$\Lambda$	2300 persons/day	New births and new residents	[51]
${\beta }_1$	$8.58\times {10}^{-9}$ day${}^{-1}$	Transmission rate before intervention	Assumed
${\beta }_2$	$3.43\times {10}^{-9}$ day${}^{-1}$	Transmission rate during and after intervention	Assumed
$\alpha$	$3.5\times {10}^{-4}$ day${}^{-1}$	Vaccination rate	[52]
$\mu$	$3\times {10}^{-5}$ persons/day	Natural death rate	[51]
${\gamma }^{-1}$	5.5 days	Incubation period	[42]
$\sigma$	0.05	Vaccine inefficacy	[53]
${\delta }^{-1}$	3.8 days	Infection time	[42]
$\kappa$	0.014	Case fatality rate	[50]
${\lambda }^{-1}$	10 days	Recovery time	[42]
${\rho }^{-1}$	15 days	Time until death	[42]

, with a standard deviation of $10\%$. In addition, we assume the observation-error standard deviations to be $5\%$ of the data value.

In order to initialize the assimilation system, we perturb the initial state using a Gaussian noise as follows:

$$x_0^i=\max \left[x_0\left(1+\eta ·\sigma \right),0\right]i=1,2,\cdots ,N_e,$$

(30)

where $\sigma \sim \mathcal{N}(0,1)$ and $\eta$ is set to $10\%$. Similarly, the parameters’ ensemble is generated by perturbing the initial values, given in Table 1, using a log-normal distribution.

To evaluate the accuracy of the assimilation system solution, we compute the Relative Mean Absolute Error (RMAE), defined as follows:

$$\mathrm{RMAE}={\displaystyle \frac{1}{n}}\sum _{j=1}^n{\displaystyle \frac{\left|y\right(j)-x(j\left)\right|}{\left|y\right(j\left)\right|}},$$

(31)

where $y\left(j\right)$ represents the observed value, $x\left(j\right)$ is the model simulation, and n is the sample size of the observed data.

4.2. Sensitivity to Ensemble Size

We first study the sensitivity of the filter, to the ensemble size, $N_e$. Computing accurate state and parameter estimates with small ensembles is desirable to avoid an important increase in computational time. We run our assimilation system using seven ensemble sizes: $N_e=10$, 20, 50, 100, 200, 500, and 1000. Figure 2 presents the performance of the EnKF in terms of RMAE for the different ensemble sizes. As expected, increasing the ensemble size reduces the discrepancy between the estimates and the data. A clear improvement over the case with 20 ensemble members is observed. However, increasing the ensemble size beyond 200 members does not notably enhance the performance of the filter, although it dramatically increases the computational cost.

Figure 3 summarizes the assimilation results for the active, recovered, and confirmed cases obtained by the EnKF using 200 ensemble numbers. All results show a good fit between the predicted and observed values, with generally minor differences. The values of the model parameters estimated by the EnKF during the assimilation process are sketched in Figure 4. The parameters do not appear to be dramatically different from the initial values, except for the reproduction number and the recovery time. We observe a sharp increase in the reproduction number between March 1st and 15th. This number remains relatively high, almost $R_0=4$, during the pre-intervention period. During and after the intervention, the value of $R_0$ remains close to 1. That suggests that the intervention was very fruitful in controlling the pandemic. The recovery time increases significantly during the first two months, peaking at $16.1$, and remains constant after that. We also see a slight decrease in the time until death to $13.75$. The other parameters show slight deviations from the proposed initial values.

4.3. Predictability Experiment

To assess the system relevance for making short-term predictions, we daily performed two-week forecasts, starting from 1st July until 17th December (the day where the vaccination campaign started in Saudi Arabia). The estimated state and parameters during the assimilation process were used to conduct these forecasts. The RMAE of the deaths, and recovered and confirmed cases were then computed and plotted in Figure 5. The model performed very well. The values of the RMAE for the forecasted recovered and confirmed cases are always less than $5\%$, within the margin of the observation errors. Slightly higher RMAE are observed during the month of July for the death cases, reaching $12\%$. However, the model achieves a lower RMAE in the next few months, with an average of $3\%$.

Figure 6 shows the final predictions performed on 3rd December for the first two weeks of December. The predicted numbers of deaths, and recovered and confirmed cases are all in very good agreement with the data. While they perfectly match the data during the first 10 days, they tend to slightly underestimate or overestimate the observed values at the end of the second week. Overall, our numerical results demonstrate the capability of the proposed system in achieving an excellent two-week prediction skill.

4.4. Impact of Vaccination

On 17 December, Saudi Arabia started a three-phase COVID-19 vaccination program, thus becoming the first Arab country to roll out the Pfizer-BioNTech vaccine. The first phase targets people aged 65 and above, as well as people most exposed to the disease. The second phase will include people aged over 50, as well as those with specific chronic diseases, such as asthma and diabetes. The rest of the population will get the vaccine at the third stage. However, the Ministry has identified those who should not get the vaccine. They are women who are pregnant or breastfeeding, women who plan on getting pregnant in the next two months, people who suffer from severe allergic reactions, and those who have been infected by the virus in the past 90 days.

As of 17 January 2021, nearly 295,000 people have been vaccinated against COVID-19, according to the Saudi Health Council [52]. That gives $\alpha =3\times {10}^{-4}$. The Pfizer-BioNtech vaccine offers up to $95\%$ protection against Covid-19 after a second dose [53], which gives $\sigma =0.05$. We start our numerical simulation on 17 December for a period of six months with and without vaccination. Numerical results are presented in Figure 7, where $\alpha =0$ corresponds to the model output without vaccination. We notice that the current pace of vaccination does not show a significant decrease in the number of cases. We observe only a reduction of $9.02\%$ in the total confirmed cases, and $7.47\%$ in the total deaths. However, the vaccination rate is expected to increase in the upcoming weeks. Therefore, we consider three more different values of $\alpha$$\left(6\times {10}^{-4},1.2\times {10}^{-3},\mathrm{and}2.4\times {10}^{-3}\right)$. As we can see in Figure 7, intensifying the vaccination campaign starts to significantly decrease the number of cases. For $\alpha =0.0006$, we observe a reduction up to $16.08\%$ and $13.47\%$ for the total confirmed cases and deaths, respectively. Finally, a notable reduction is observed for $\alpha =0.0024$. Results show a decline by $38.67\%$ and $33.65\%$, respectively.

5. Conclusions

This paper presented an extended SEIR model with seven stages of infection, including vaccination, to simulate and predict the development of the novel COVID-19 outbreak. Initially, a mathematical analysis was carried out to illustrate the non-negativity, boundedness, epidemic equilibrium, existence, and uniqueness of the endemic equilibrium, and basic reproduction number of the proposed model. To mitigate the imperfection of SEIR models and the poorly known parameters, we first applied a data assimilation framework based on a Joint-EnKF for state-parameters estimation in order to enhance the model forecast skill. We assimilated real data of the number of deaths, the number of recovered cases, and the number of active cases. In the second step, the updated state and parameters were used to perform predictions two weeks ahead in order to assess the relevance of the system of making short-term predictions of the COVID-19 pandemic. Finally, we investigated the impact of vaccination on the spread of the disease.

The proposed methodology was validated with cases reported from Saudi Arabia. In the assimilation process, the estimated numbers of active, recovered, deaths, and confirmed cases were accurately predicted by the model, with generally minor differences. In the predictability experiments, the model provided encouraging results for short-term predictions, with a relatively small ensemble size. Moreover, model predictions can be updated weekly or every two weeks to accommodate short- and medium-term predictions. Results also show that intensifying the vaccination campaign significantly decreases the number of confirmed cases and deaths. The proposed model can serve as a tool for health authorities to plan, prepare, and take appropriate measures and decisions to control the pandemic.