# A Mathematical Model to Investigate the Transmission of COVID-19 in the Kingdom of Saudi Arabia

Authors: Fehaid Salem Alshammari
Journal: Computational and Mathematical Methods in Medicine
 (2020)

---

## Abstract

Since the first confirmed case of SARS-CoV-2 coronavirus (COVID-19) on March 02, 2020, Saudi Arabia has not reported quite a rapid COVD-19 spread as seen in America and many European countries. Possible causes include the spread of asymptomatic COVID-19 cases. To characterize the transmission of COVID-19 in Saudi Arabia, a susceptible, exposed, symptomatic, asymptomatic, hospitalized, and recovered dynamical model was formulated, and a basic analysis of the model is presented including model positivity, boundedness, and stability around the disease-free equilibrium. It is found that the model is locally and globally stable around the disease-free equilibrium whenR0<1. The model parameterized from COVID-19 confirmed cases reported by the Ministry of Health in Saudi Arabia (MOH) from March 02 till April 14, while some parameters are estimated from the literature. The numerical simulation showed that the model predicted infected curve is in good agreement with the real data of COVID-19-infected cases. An analytical expression of the basic reproduction number R0 is obtained, and the numerical value is estimated as R0≈2.7.


---

## Body

## 1. Introduction

As of April 22, 2020, more than 12772 cases and 114 deaths of coronavirus disease 2019 (COVID-19) caused by the SARS-CoV-2 virus had been confirmed in the Kingdom of Saudi Arabia. Since March 04 [1], control measures have been implemented within Saudi Arabia to control the spread of the disease. Isolation of confirmed cases and contact tracing are crucial parts of these measures, which are common interventions for controlling infectious disease outbreaks [2–4]. For example, the severe acute respiratory syndrome (SARS) outbreak and the Middle East respiratory syndrome (MERS) were controlled through tracing suspected cases and isolating confirmed cases because the majority of transmission occurred concurrent or after symptom onset [3–5].However, it is unknown if transmission of COVID-19 can occur before symptom onset, which could decrease the effectiveness of isolation and contact tracing [2, 3, 5].In this paper, the impact of asymptomatic COVID-19 cases on the spread of the disease is considered using a modified version of the susceptible-exposed-infected-recovered (SEIR) dynamical model, along with the data of COVID-19 cases reported daily by the Ministry of Health in Saudi Arabia (MOH). Other main objectives of this paper include obtaining an analytical expression and numerical estimation of the basic reproduction numberR0 of COVID-19 in Saudi Arabia and estimating the maximum required number of hospital beds and intensive care units (ICU). The rest of this study is organized as follows: The model establishment is presented in Section 2. Basic analysis of the model, including local and global stability results of the disease-free equilibrium, is explored in Section 3. The model fitting to daily reported cases and estimation of parameters is given in Section 4. Numerical results and discussion are presented in Section 5. Finally, brief concluding remarks are given in Section 6.


## 2. Model Establishment

The Saudi populationQ will be divided into six categories: susceptible S, exposed E, symptomatic Y, asymptomatic N, hospitalized H, and recovered R individuals (SEYNHR). Individuals move from the susceptible compartment S to the exposed compartment E after interacting with infected individuals with transmission rates β1, β2, and β3 as shown in Figure 1. COVID-19 is known to have an incubation period, from 2 to 14 days, between exposure and development of symptoms [6, 7]. After this period, exposed individual transits from the compartment E to either compartment Y at a rate α or compartment N at a rate α1−γ. An individual could move from compartment N to Y at a rate K if they show symptoms. Once an individual becomes infected with the coronavirus that causes COVID-19, that individual develops immunity against the virus with a rate ΦY or the individual will be hospitalized with a rate of ε or dies because of the disease with a rate of μ1. When an individual becomes hospitalized, that individual receives treatment and develops immunity against the virus with a rate r or dies because of the disease with a rate μ2.Figure 1
Schematic diagram of theSEYNHR compartment model. The arrows, except the black ones, represent a progression from one compartment to the next.As shown in Figure1, the SEYNHR model has six compartments; therefore, a discrete dynamical system consisting of six nonlinear differential equations is formed as follows:
(1)dSdt=A−β1Y+β2N+β3HS−μS,(2)dEdt=β1Y+β2N+β3HS−α+μE,(3)dYdt=α1−γE−ΦY+ε+μ1+μY+KN,(4)dNdt=αγE−ΦN+K+μN,(5)dHdt=εY−r+μ2+μH,(6)dRdt=ΦYY+ΦNN+rH−μR,where Qt=St+Et+Yt+Nt+Ht+Rt and Qt is the Saudi population at time t. The next-generation matrix is used in the next section to derive an analytical expression of the basic reproduction number R0, for the compartmental model above. Calculating R0 is a useful metric for assessing the transmission potential of an emerging COVID-19 in Saudi Arabia.


## 3. Basic Analytical Results

In this section, the positivity and boundedness of the proposedSEYNHR model solution (6) and stability of the model around the disease-free equilibrium is investigated. Furthermore, an analytical expression of the basic reproduction number R0 is established.


### 3.1. Positivity and Boundedness

To show that theSEYNHR model (6) is biologically well-behaved or epidemiologically meaningful, we must show that all its state variables are nonnegative and bounded for t>0.
(7)dSdt=A−β1Yt+β2Nt+β3HtSt−μSt.Lemma 1.
Fort>0 and initial values P0≥0, where Pt=St,Et,Yt,Nt,Ht,Rt, the solution of the SEYNHR model (6) is nonnegative if they exist.Proof.
Assumet1=t>0:Pt>0∈0,t, then it follows from the first equation of (6) asIt can be written as(8)dSdtStexpμt+∫0tβ1Yζ+β2Nζ+β3HζSζdζ−Aexpμt+∫0tβ1Yζ+β2Nζ+β3HζSζdζ=0.Hence,(9)St1expμt1+∫0t1β1Yζ+β2Nζ+β3HζSζdζ−S0=∫0t1Aexpμp+∫0pβ1Yζ+β2Nζ+β3HζSζdζdp,which can be found as below
(10)St1=S0exp−μt1+∫0t1β1Yζ+β2Nζ+β3HζSζdζ+exp−μt1+∫0t1β1Yζ+β2Nζ+β3HζSζdζ×∫0t1Aexpμp+∫0pβ1Yζ+β2Nζ+β3HζSζdζdp>0.Hence,St is nonnegative for t>0. In a similar way, it can be shown that E>0, Y>0, N>0, H>0, and R>0 for t>0.
(11)dQdt=A−μQt−μ1Y−μ2H≤A−μQt.Lemma 2.
Ω=St,Et,Yt,Nt,Ht,Rt∈R+6∪0:0<St+Et+Yt+Nt+Ht+Rt≤A/μ is the positively invariant region of the model (6) with nonnegative initial conditions in R+6.Proof.
SupposeQt=St+Et+Yt+Nt+Ht+Rt holds for t≥0, then we getIt is obvious that(12)Qt≤Q0e−μt+Aμ1−e−μt.Thus,limt⟶∞supQt≤A/μ. Furthermore, dQt/dt<0 if Qt>A/μ. This shows that solutions of the model (6) point towards Ω. Hence, Ω is positively invariant and solutions of (6) are bounded.


### 3.2. Basic Reproduction NumberR0

It can be easily verified that system (6) always has a disease-free equilibrium (DFE), which we will denote by E0, and it is given by
(13)E0=Aμ,0,0,0,0,0.Next, we investigate an important concept in epidemiology which is the basic reproduction number, defined as “the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual” [8]. The next-generation method is used to calculate R0; more details about the method are in [8]. The system (6) can be rewritten as follows:
(14)dwdt=Φw−Ψw,where F≔F1,F2,F3,F4,F5,F6T and V≔V1,V2,V3,V4,V5,V6T, or more explicitly
(15)ĖẎṄḢṘṠ=β1Y+β2N+β3HS00000−α+μEY+ε+μ1+μY−α1−γ∗E−KNΦN+K+μN−αγEr+μ2+μH−εYP+μR−ΦY−ΦN−rHβ1Y+β2N+β3HS−A−PR.The Jacobian matrices ofF and V evaluated at the solution E0∗=0,0,0,0,0,A/μT are given by
(16)F=∂Φi∂wjE0∗1≤i,j≤4=0β1Aμβ2Aμβ3Aμ000000000000,V=∂Ψi∂wjE0∗1≤i,j≤4=α+μ000−α1−γΦY+ε+μ1+μ−K0−αγ0ΦN+K+μ00−ε0r+μ2+μ.Direct calculations show that(17)V−1=1α+μ000m1m4m60m20m70m3m5m81r+μ2+μ,where
(18)m1=α1−γK+α1−γμ+α1−γΦN+KαγΦY+ε+μ1+μΦN+K+μα+μ,m2=αγα+μΦN+K+μ,m3=εα1−γK+α1−γμ+α1−γΦN+KαγΦY+ε+μ1+μΦN+K+μα+μr+μ2+μ,m4=1ΦY+ε+μ1+μ,m5=εΦY+ε+μ1+μr+μ2+μ,m6=KΦN+K+μΦY+ε+μ1+μ,m7=1ΦN+K+μ,m8=εKΦY+ε+μ1+μΦN+K+μr+μ2+μ.Denoting the4×4 identity matrix by I, the characteristic polynomial Γλ of the matrix FV−1 is given by
(19)Γλ=detFV−1−λI,=−λ3ΦN+K+μβ1μ+r+μ2β1+β3εAα1−γ+B1+B2Aαγ+λ4μΦY+ε+μ1+μΦN+K+μα+μr+μ2+μ,where
(20)B1=β2μ2+μ1+μ2+ΦY+ε+rβ2+β1Kμ+Kr+μ2β1,B2=r+μ2ΦY+ε+μ1β2+β3εK.The solutionsλ1,2,3,4 are given by
(21)0,0,0,ΦN+K+μβ1μ+r+μ2β1+β3εAα1−γ+B1+B2AαγμΦY+ε+μ1+μΦN+K+μα+μr+μ2+μ.Therefore, the reproduction numberR0 for the SEYNHR model (6) is given by
(22)R0=maxλ1,λ2,λ3,λ4,=ΦN+K+μβ1μ+r+μ2β1+β3εAα1−γ+B1+B2AαγμΦY+ε+μ1+μΦN+K+μα+μr+μ2+μ.For simplicity, let us denote(23)ρ1=α+μ,ρ2=ΦY+ε+μ1+μ,ρ3=ΦN+K+μ,ρ4=r+μ2+μ.In order to interpret the basic reproduction numberR0, we split the above expression of R0 as follows:
(24)R0=RYH+REYNH1+REYNH2,where
(25)RYH=α1−γ×Aμ×1ρ1×1ρ1ρ4×β1ρ4+β3ε,REYNH1=αγ×Aμ×1ρ1ρ2ρ3ρ4×β3μ21ρ1ρ2ρ3ρ4×β1Kr+μ2+1ρ1ρ2ρ3ρ4×μ1+μ2+ΦY+ε+rβ2+β1Kμ,REYNH2=αγ×Aμ1ρ1ρ2ρ3ρ4×β2r+μ2ρ2−μ+1ρ1ρ2ρ3ρ4×β1εK.The first and second terms inRYH show the probability of the total population becoming symptomatic upon infection multiplied by mean symptomatic, exposed, and hospitalized infectious periods multiplied by symptomatic and hospitalized contact rates. The first and second terms in REYNH1 show the probability of the total population becoming asymptomatic multiplied by mean exposed, symptomatic, asymptomatic, and hospitalized infectious periods multiplied by symptomatic, asymptomatic, and hospitalized contact rates. Furthermore, the first and second terms in REYNH2 show the probability of the total population becoming asymptomatic multiplied by mean exposed, symptomatic, asymptomatic, and hospitalized infectious periods multiplied by asymptomatic and hospitalized contact rates. The obtained analytical expression of R0 shows that asymptomatic as well as symptomatic cases of COVID-19 could drive the growth of the pandemic in Saudi Arabia.


### 3.3. Stability of Disease-Free Equilibrium

In this section, the local and global stability of the model (6) around the disease-free equilibrium is explored. The epidemiological implication of the local stability result of DFE is that a small influx of COVID-19 infection cases will not generate a COVID-19 outbreak if R0<1 while the global stability result of DFE indicates that any influx of COVID-19 infection cases will not generate a COVID-19 outbreak if R0<1.
(26)JE0=−μ0−β1Aμ−β2Aμ−β3Aμ00−α+μβ1Aμβ2Aμβ3Aμ00α1−γ−ΦY+ε+μ1+μK000αγ0−ΦN+K+μ0000ε0−r+μ2+μ000ΦYΦNr−μ.Theorem 3.
The DFEE0 of the system (6) is locally asymptotically stable if R0<1 and unstable otherwise.Proof.
The Jacobian matrixJE0 obtained at the DFE E0 is as follows:Clearly, the eigenvalues−μ and −μ are negative real numbers. The remaining eigenvalues can be obtained through the following equation:
(27)λ4+C1λ3+C2λ2+C3λ+C4=0,where
(28)C1=ΦY+ΦN+K+4μ+α+ε+μ1+μ2+r,C2=−Aαβ11−γμ+β1ρ1ρ2ρ32ρ4ψ,C3=−2Aβ1α+Aα1−γμ+ρ1ρ2ρ32ρ4ψ+2β1ρ1ρ2ρ32ρ4ψ,C4=−Aβ1α+Aα1−γμ+ρ1ρ2ρ32ρ4ψ+β1ρ1ρ2ρ32ρ4ψ+positiveterms,where ψ=β1μ+r+μ2β1+β3ε and C1 is obviously positive. We can rewrite R0 as
(29)R0=Aα1−γ+positivetermsμρ1ρ2ρ32ρ4ψ.Clearly ifR0<1, then C2 is positive. Similarly, we can show that C3 and C4 are all positive if R0<1. Further, it is easy to show the remaining Routh-Hurtwiz condition for the fourth-order polynomial (27). Thus, the DFE is locally asymptotically stable if R0<1.The global stability of DFEE0 of the COVID-19 transmission model is studied in the following result.
(30)Πt=W1QE+W2Y+W3N+W4H,where Wi for i=1,2,3,4, are positive constants. Differentiating the function Πt with respect to t and using the solutions of system (6), we obtain
(31)Πtdt=W1β1Y+β2N+β3HSQ−ρ1EQ+W2α1−γE−ρ2Y+KY+W3αγ−ρ3N+W4εY−ρ4H≤W1β1Y+β2N+β3H−ρ1E+W2α1−γE−ρ2Y+KY+W3αγ−ρ3N+W4εY−ρ4H,asS+EQ<1.=W2α1−γ−ρ1+W3αγE+W1β1−W2ρ2+W4εY+W1β2−W3ρ3N+W1β3−W4ρ4H=ρ1W2α1−γ+W3αγρ1−1E+W1β1−W2ρ2+W4εYρ1W2α1−γ+W3αγρ1−1E+W1β1−W2ρ2+W4εY+W1β2−W3ρ3N+W1β3−W4ρ4HTheorem 4.
The DFEE0 of the system (6) is globally asymptotically stable if R0<1 and unstable otherwise.Proof.
Let us consider the following Lyapunov function, described in [9]:Now choosing(32)W1=1,W2=Aβ1μ+r+μ2β1+β3εμρ1ρ2ρ4,W3=AB1+B2μρ1ρ2ρ3ρ4,W4=1.We obtain(33)Πtdt≤ρ1R0−1E.Hence, ifR0<1, then Πt/dt<0. Therefore, the largest compact invariant set in Ω is the singleton set E0, and using LaSalle’s invariant principle [10], E0 is globally asymptotically stable in Ω.


## 3.1. Positivity and Boundedness

To show that theSEYNHR model (6) is biologically well-behaved or epidemiologically meaningful, we must show that all its state variables are nonnegative and bounded for t>0.
(7)dSdt=A−β1Yt+β2Nt+β3HtSt−μSt.Lemma 1.
Fort>0 and initial values P0≥0, where Pt=St,Et,Yt,Nt,Ht,Rt, the solution of the SEYNHR model (6) is nonnegative if they exist.Proof.
Assumet1=t>0:Pt>0∈0,t, then it follows from the first equation of (6) asIt can be written as(8)dSdtStexpμt+∫0tβ1Yζ+β2Nζ+β3HζSζdζ−Aexpμt+∫0tβ1Yζ+β2Nζ+β3HζSζdζ=0.Hence,(9)St1expμt1+∫0t1β1Yζ+β2Nζ+β3HζSζdζ−S0=∫0t1Aexpμp+∫0pβ1Yζ+β2Nζ+β3HζSζdζdp,which can be found as below
(10)St1=S0exp−μt1+∫0t1β1Yζ+β2Nζ+β3HζSζdζ+exp−μt1+∫0t1β1Yζ+β2Nζ+β3HζSζdζ×∫0t1Aexpμp+∫0pβ1Yζ+β2Nζ+β3HζSζdζdp>0.Hence,St is nonnegative for t>0. In a similar way, it can be shown that E>0, Y>0, N>0, H>0, and R>0 for t>0.
(11)dQdt=A−μQt−μ1Y−μ2H≤A−μQt.Lemma 2.
Ω=St,Et,Yt,Nt,Ht,Rt∈R+6∪0:0<St+Et+Yt+Nt+Ht+Rt≤A/μ is the positively invariant region of the model (6) with nonnegative initial conditions in R+6.Proof.
SupposeQt=St+Et+Yt+Nt+Ht+Rt holds for t≥0, then we getIt is obvious that(12)Qt≤Q0e−μt+Aμ1−e−μt.Thus,limt⟶∞supQt≤A/μ. Furthermore, dQt/dt<0 if Qt>A/μ. This shows that solutions of the model (6) point towards Ω. Hence, Ω is positively invariant and solutions of (6) are bounded.


## 3.2. Basic Reproduction NumberR0

It can be easily verified that system (6) always has a disease-free equilibrium (DFE), which we will denote by E0, and it is given by
(13)E0=Aμ,0,0,0,0,0.Next, we investigate an important concept in epidemiology which is the basic reproduction number, defined as “the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual” [8]. The next-generation method is used to calculate R0; more details about the method are in [8]. The system (6) can be rewritten as follows:
(14)dwdt=Φw−Ψw,where F≔F1,F2,F3,F4,F5,F6T and V≔V1,V2,V3,V4,V5,V6T, or more explicitly
(15)ĖẎṄḢṘṠ=β1Y+β2N+β3HS00000−α+μEY+ε+μ1+μY−α1−γ∗E−KNΦN+K+μN−αγEr+μ2+μH−εYP+μR−ΦY−ΦN−rHβ1Y+β2N+β3HS−A−PR.The Jacobian matrices ofF and V evaluated at the solution E0∗=0,0,0,0,0,A/μT are given by
(16)F=∂Φi∂wjE0∗1≤i,j≤4=0β1Aμβ2Aμβ3Aμ000000000000,V=∂Ψi∂wjE0∗1≤i,j≤4=α+μ000−α1−γΦY+ε+μ1+μ−K0−αγ0ΦN+K+μ00−ε0r+μ2+μ.Direct calculations show that(17)V−1=1α+μ000m1m4m60m20m70m3m5m81r+μ2+μ,where
(18)m1=α1−γK+α1−γμ+α1−γΦN+KαγΦY+ε+μ1+μΦN+K+μα+μ,m2=αγα+μΦN+K+μ,m3=εα1−γK+α1−γμ+α1−γΦN+KαγΦY+ε+μ1+μΦN+K+μα+μr+μ2+μ,m4=1ΦY+ε+μ1+μ,m5=εΦY+ε+μ1+μr+μ2+μ,m6=KΦN+K+μΦY+ε+μ1+μ,m7=1ΦN+K+μ,m8=εKΦY+ε+μ1+μΦN+K+μr+μ2+μ.Denoting the4×4 identity matrix by I, the characteristic polynomial Γλ of the matrix FV−1 is given by
(19)Γλ=detFV−1−λI,=−λ3ΦN+K+μβ1μ+r+μ2β1+β3εAα1−γ+B1+B2Aαγ+λ4μΦY+ε+μ1+μΦN+K+μα+μr+μ2+μ,where
(20)B1=β2μ2+μ1+μ2+ΦY+ε+rβ2+β1Kμ+Kr+μ2β1,B2=r+μ2ΦY+ε+μ1β2+β3εK.The solutionsλ1,2,3,4 are given by
(21)0,0,0,ΦN+K+μβ1μ+r+μ2β1+β3εAα1−γ+B1+B2AαγμΦY+ε+μ1+μΦN+K+μα+μr+μ2+μ.Therefore, the reproduction numberR0 for the SEYNHR model (6) is given by
(22)R0=maxλ1,λ2,λ3,λ4,=ΦN+K+μβ1μ+r+μ2β1+β3εAα1−γ+B1+B2AαγμΦY+ε+μ1+μΦN+K+μα+μr+μ2+μ.For simplicity, let us denote(23)ρ1=α+μ,ρ2=ΦY+ε+μ1+μ,ρ3=ΦN+K+μ,ρ4=r+μ2+μ.In order to interpret the basic reproduction numberR0, we split the above expression of R0 as follows:
(24)R0=RYH+REYNH1+REYNH2,where
(25)RYH=α1−γ×Aμ×1ρ1×1ρ1ρ4×β1ρ4+β3ε,REYNH1=αγ×Aμ×1ρ1ρ2ρ3ρ4×β3μ21ρ1ρ2ρ3ρ4×β1Kr+μ2+1ρ1ρ2ρ3ρ4×μ1+μ2+ΦY+ε+rβ2+β1Kμ,REYNH2=αγ×Aμ1ρ1ρ2ρ3ρ4×β2r+μ2ρ2−μ+1ρ1ρ2ρ3ρ4×β1εK.The first and second terms inRYH show the probability of the total population becoming symptomatic upon infection multiplied by mean symptomatic, exposed, and hospitalized infectious periods multiplied by symptomatic and hospitalized contact rates. The first and second terms in REYNH1 show the probability of the total population becoming asymptomatic multiplied by mean exposed, symptomatic, asymptomatic, and hospitalized infectious periods multiplied by symptomatic, asymptomatic, and hospitalized contact rates. Furthermore, the first and second terms in REYNH2 show the probability of the total population becoming asymptomatic multiplied by mean exposed, symptomatic, asymptomatic, and hospitalized infectious periods multiplied by asymptomatic and hospitalized contact rates. The obtained analytical expression of R0 shows that asymptomatic as well as symptomatic cases of COVID-19 could drive the growth of the pandemic in Saudi Arabia.


## 3.3. Stability of Disease-Free Equilibrium

In this section, the local and global stability of the model (6) around the disease-free equilibrium is explored. The epidemiological implication of the local stability result of DFE is that a small influx of COVID-19 infection cases will not generate a COVID-19 outbreak if R0<1 while the global stability result of DFE indicates that any influx of COVID-19 infection cases will not generate a COVID-19 outbreak if R0<1.
(26)JE0=−μ0−β1Aμ−β2Aμ−β3Aμ00−α+μβ1Aμβ2Aμβ3Aμ00α1−γ−ΦY+ε+μ1+μK000αγ0−ΦN+K+μ0000ε0−r+μ2+μ000ΦYΦNr−μ.Theorem 3.
The DFEE0 of the system (6) is locally asymptotically stable if R0<1 and unstable otherwise.Proof.
The Jacobian matrixJE0 obtained at the DFE E0 is as follows:Clearly, the eigenvalues−μ and −μ are negative real numbers. The remaining eigenvalues can be obtained through the following equation:
(27)λ4+C1λ3+C2λ2+C3λ+C4=0,where
(28)C1=ΦY+ΦN+K+4μ+α+ε+μ1+μ2+r,C2=−Aαβ11−γμ+β1ρ1ρ2ρ32ρ4ψ,C3=−2Aβ1α+Aα1−γμ+ρ1ρ2ρ32ρ4ψ+2β1ρ1ρ2ρ32ρ4ψ,C4=−Aβ1α+Aα1−γμ+ρ1ρ2ρ32ρ4ψ+β1ρ1ρ2ρ32ρ4ψ+positiveterms,where ψ=β1μ+r+μ2β1+β3ε and C1 is obviously positive. We can rewrite R0 as
(29)R0=Aα1−γ+positivetermsμρ1ρ2ρ32ρ4ψ.Clearly ifR0<1, then C2 is positive. Similarly, we can show that C3 and C4 are all positive if R0<1. Further, it is easy to show the remaining Routh-Hurtwiz condition for the fourth-order polynomial (27). Thus, the DFE is locally asymptotically stable if R0<1.The global stability of DFEE0 of the COVID-19 transmission model is studied in the following result.
(30)Πt=W1QE+W2Y+W3N+W4H,where Wi for i=1,2,3,4, are positive constants. Differentiating the function Πt with respect to t and using the solutions of system (6), we obtain
(31)Πtdt=W1β1Y+β2N+β3HSQ−ρ1EQ+W2α1−γE−ρ2Y+KY+W3αγ−ρ3N+W4εY−ρ4H≤W1β1Y+β2N+β3H−ρ1E+W2α1−γE−ρ2Y+KY+W3αγ−ρ3N+W4εY−ρ4H,asS+EQ<1.=W2α1−γ−ρ1+W3αγE+W1β1−W2ρ2+W4εY+W1β2−W3ρ3N+W1β3−W4ρ4H=ρ1W2α1−γ+W3αγρ1−1E+W1β1−W2ρ2+W4εYρ1W2α1−γ+W3αγρ1−1E+W1β1−W2ρ2+W4εY+W1β2−W3ρ3N+W1β3−W4ρ4HTheorem 4.
The DFEE0 of the system (6) is globally asymptotically stable if R0<1 and unstable otherwise.Proof.
Let us consider the following Lyapunov function, described in [9]:Now choosing(32)W1=1,W2=Aβ1μ+r+μ2β1+β3εμρ1ρ2ρ4,W3=AB1+B2μρ1ρ2ρ3ρ4,W4=1.We obtain(33)Πtdt≤ρ1R0−1E.Hence, ifR0<1, then Πt/dt<0. Therefore, the largest compact invariant set in Ω is the singleton set E0, and using LaSalle’s invariant principle [10], E0 is globally asymptotically stable in Ω.


## 4. Model Fitting and Parameter Estimation

For the parameterizations of the model (6), we consider some of the parameter values from the literature and the rest are fitted to the data of daily total COVID-19 cases in Saudi Arabia from March 02, 2020, till April 14, 2020, using the nonlinear least square method implemented in MATLAB.Considering the time unit to be days, we can estimate the following parameters:(i)
Natural death rateμ: the average life span of Saudi people is 74.87 years; therefore, we have μ=1/74.87×365=3.6593×10−5 per day [11](ii)
The birth rateA is obtained from A/μ=Q0, where the total Saudi population Q0 is 34218169; therefore, A=1252 per day [11]Mean asymptomatic infectious periodΦN is assumed to be the same as the mean asymptomatic infectious period ΦY because there is no estimation available in the literature [5, 12]. The remaining system parameters are either estimated from literature or obtained using a nonlinear least-square curve fitting in MATLAB and given in Table 1. The basic steps in the parameter estimation procedure are described in [13, 14] and summarized as follows: system (6) can be expressed as
(34)dudt=Zt,u,θ,ut0=u0.Table 1
Parameters used in the simulations.
ParameterDescriptionDimensionValueA 95% confidence intervalSourceNPopulation of Saudi Arabia in 2019.Individual34218169—[11]ABirth rate in Saudi Arabia in 2019.Individual/day1252—[11]μDeath rate in Saudi Arabia in 2019.Day−13.6593x10−5—[11]β1Effective contact rate from symptomatic to susceptible.Day−10.35[0.267, 0.5]Fittedβ2Effective contact rate from asymptomatic to susceptible.Day−10.7[0.53, 1]Fittedβ3Effective contact rate from hospitalized to susceptible.Day−10.18[0.13, 0.26]Fittedα−1Mean latent period.Days5.1—[23]γProbability of becoming asymptomatic upon infection.n/a0.86—[24]ΦY−1Mean symptomatic infectious period.Days8—[25]εRate of symptomatic becomes hospitalized.Day−10.125—[1]μ1Death rate of symptomatic patients.Day−10.392—[26]KRate of asymptomatic becomes symptomatic.Day−10.12[0.1, 0.125]FittedΦN−1Mean asymptomatic infectious period.Days8—AssumedrRate of recovered hospitalized patients.Day−10.1—[6]μ2Death rate of hospitalized patients.Day−10.392—[1]The functionZ is dependent on t, the vectors of the dependent variables u, and unknown parameters θ to be estimated. The least-square technique estimates the best values of the model parameters by minimizing the error between the reported data points y^ti and the solution of the model yti associated with the model parameters θ. The objective function used in the minimization procedure is given as n where it denotes the available real data points. To obtain the model parameters, we minimize the following objective function:
(35)θ^=∑i=1nyti−y^ti2.(36)minθ^subjecttoequation4.The model (6) is solved using MATLAB solver ode45 which uses the Runge-Kutta methods to solve the initial value problem. Then, the lsqcurvefit package was implemented to fit model (6) to COVID-19 data of confirmed cases from March 02 till April 14, to estimate the parameters using the described approach above. In Figure 2 and Figure 3(a), the numerical simulation showed that the model predicted infected curve is in good agreement to the real data of infected cases.Figure 2
Data fitting to the reported cases in Saudi Arabia from March 02, 2020, till April 14, 2020, using model (6) and the nonlinear least-square method.Figure 3
Numerical results of theSEYNHR model with the parameter listed in Table 1. (a) Shows the estimated number of symptomatic COVID-19 cases, with the published data by the Ministry of Health in Saudi Arabia of confirmed COVID-19 cases (blue circles). (b) Shows the estimated susceptible, exposed, symptomatic, asymptomatic, hospitalized, and recovered subpopulations. (c) Shows estimations of the hospitalized cases and the required ICU beds (black dashed line). The red line represents the number of hospital beds available, while the red dashed line represents the number of ICU beds available in Saudi Arabia.
(a)(b)(c)


## 5. Numerical Results and Discussion

To illustrate the numerical results, we take the initial population and subpopulations asQ0=34218169, Y0=2, S0=Q0−Y0=34218167, and E0=N0=H0=R0. The estimated values of parameters fitted from COVID-19 cases from March 02 till April 14 published by the Ministry of Health in Saudi Arabia (MOH) are presented in Table 1. Figure 3 shows the simulation of the model (6) with parameter values in Table 1. Figure 3(b) shows that about 18% of the entire Saudi population will be asymptomatic in the last week of May 2020, and about 17% will be exposed in the third week of May. The percentage of the entire population being symptomatic at any time will not exceed 1%, which is estimated to occur in the third week of May.Moreover, as shown in Figure3(c), about 60000 hospital beds and 18000 ICU beds are required (assuming 30% of the hospitalized cases need ICU [6]) immediately after the second week of May. As of April 2020, the Ministry of Health in Saudi Arabia designated 25 hospitals for COVID-19-infected patients with up to 80000 beds and 8000 intensive care unit (ICU) beds [15] and therefore extra 10000 ICU beds could be required.The parametersβ1, β2, and β3 represent the impact of effective contact on the disease transmission as shown in Figures 4, 5, and 6. The effective contact rate from asymptomatic to susceptible β2 has the greatest impact on the number of new COVID-19 symptomatic cases Yt as shown in Figure 5. This is because of the large number of asymptomatic cases in comparison to symptomatic cases in the same timescale as shown in Figure 3(b). Thus, social distancing (reduction of β2) is an effective strategy to control the spread of the disease. The effective contact rate from hospitalized to susceptible β3 has no impact on the number of new COVID-19 symptomatic cases Yt as shown in Figure 6 because of the small number of hospitalized cases in comparison to asymptomatic and symptomatic cases. Figure 7 shows that increasing the rate of asymptomatic becomes symptomatic K decreases the number of new COVID-19 symptomatic cases Yt because decreasing K implies increasing the contact time between individuals with COVID-19 (with no treatment) and susceptible individuals.Figure 4
Effective contact rate from symptomatic to susceptibleβ1 (a measure of social-distancing effectiveness) on the number of new COVID-19-infected cases Yt. The baseline value of β1 is as given in Table 1.Figure 5
Effective contact rate from asymptomatic to susceptibleβ2 (a measure of social-distancing effectiveness) on the number of new COVID-19-infected cases Yt. The baseline value of β2 is as given in Table 1.Figure 6
Effective contact rate from hospitalized to susceptibleβ3 on the number of new COVID-19-infected cases Yt. The baseline value of β3 is as given in Table 1.Figure 7
Impact of the rate of asymptomatic becomes symptomaticK on the number of new COVID-19-infected cases Yt.Based on estimated and measured parameter values, the basic reproduction number is estimated asR0≈2.7. The variation of the basic reproduction number R0 for different values of β1 and β2 and β1 and K are shown in the heat maps in Figure 8. The upper heat map shows that decreasing the effective contact rate from symptomatic to susceptible β1 and the effective contact rate from asymptomatic to susceptible β1 will decrease the spread of the disease as R0 decreases and vice versa. The lower heat map of Figure 8 shows that increasing K increases R0 and vice versa. In reality, R0 is not a biological constant; it could fluctuate daily depending on environmental and social factors such as the percentage of the entire susceptible population wearing a suitable medical mask and practicing physical distancing. In the literature, estimates of R0 vary greatly: from 1 to 6 [13, 16–22] up to 26.5 [12]. This variation is because of the different assumptions and factors they had considered in their calculations.Figure 8
Heat maps showing the variation ofR0 for different parameter values: the upper heat map shows the variation of R0 for different values for β1 and β2, and the lower heat map shows the variation of R0 for different values for β1 and K.Parameters with the highest degree of uncertainty are the effective contact rates from symptomatic to susceptibleβ1, from asymptomatic to susceptible β2, and from hospitalized to susceptible β3 (expected to be a fraction of β1 because of the protective measures in hospitals); the rate of asymptomatic becomes symptomatic K as well as the mean asymptomatic infectious period ΦN−1.The maximum estimated value ofβ1 is 0.5 per day which is one half of the value reported by Li et al. [16]. This could be a reasonable estimation as we have not seen a similar scenario in Saudi Arabia after 5 weeks since reaching 100 confirmed cases on the 14th of March (week 7 since the first case) as we have seen in many other countries like China, America, and different European countries in the same timescale. This could be a result of the precautionary measures taken by the Saudi authorities, including the closure of schools and universities that started as early as March 08 (six days after the first confirmed COVID-19 case in Saudi Arabia).


## 6. Conclusion

The COVID-19 pandemic has rapidly spread out to most countries around the world. As of April 22, 2020, more than 12772 cases had been confirmed in Saudi Arabia alone. In this study, a mathematical model is formulated in order to study the transmission of COVID-19 in Saudi Arabia. Some basic properties of the model are presented in Section3 including model positivity, boundedness, and local and global stability around the disease-free equilibrium. It is proven that the disease-free equilibrium is locally and globally stable when R0<1. The model parameterized from COVID-19 confirmed cases reported by the Ministry of Health in Saudi Arabia (MOH) from March 02 till April 14, while some parameters are estimated from the literature. The numerical simulation showed that the model predicted infected curve is in good agreement to the real data of infected cases. An analytical expression of the basic reproduction number R0 is obtained in Section 3, and the numerical value is estimated as R0≈2.7.The contribution of undocumented COVID-19 infections (asymptomatic cases) on the transmission of the disease deserves further studies and investigations. The obtained analytical expression ofR0, along with the numerical result in Figure 5, shows that asymptomatic transmission of COVID-19 could drive the growth of the pandemic in Saudi Arabia. Therefore, more testing and social distancing are needed to identify COVID-19 patients and to contain the spread of the disease.



---

*Source: 9136157-2020-10-12.xml*