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

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 when R 0 < 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 R 0 is obtained, and the numerical value is estimated as R 0 ≈ 2.7.

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 number R0 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 population Q 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 Φ 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 the SEYNHR compartment model. The arrows, except the black ones, represent a progression from one compartment to the next.

As shown in Figure 1, the SEYNHR model has six compartments; therefore, a discrete dynamical system consisting of six nonlinear differential equations is formed as follows: where Q(t) = S(t) + E(t) + Y(t) + N(t) + H(t) + R(t) and Q(t) 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 proposed SEYNHR 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 the SEYNHR model (6) is biologically well-behaved or epidemiologically meaningful, we must show that all its state variables are nonnegative and bounded for t > 0.

Lemma 1 .

For t > 0 and initial values P(0) ≥ 0, where P(t) = (S(t), E(t), Y(t), N(t), H(t), R(t)), the solution of the SEYNHR model (6) is nonnegative if they exist.

Proof

Assume t1 = {t > 0 : P(t) > 0 ∈ [0, t]}, then it follows from the first equation of (6) as

It can be written as

Hence, which can be found as below

Hence, S(t) 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.

Lemma 2 .

Ω = {(S(t), E(t), Y(t), N(t), H(t), R(t)) ∈ R+6 ∪ {0}: 0 < S(t) + E(t) + Y(t) + N(t) + H(t) + R(t) ≤ A/μ} is the positively invariant region of the model (6) with nonnegative initial conditions in R+6.

Suppose Q(t) = S(t) + E(t) + Y(t) + N(t) + H(t) + R(t) holds for t ≥ 0, then we get

It is obvious that

Thus, limsup Q(t) ≤ A/μ. Furthermore, (dQ(t))/dt < 0 if Q(t) > 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 Number R0

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

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: where F≔(F1, F2, F3, F4, F5, F6) and V≔(V1, V2, V3, V4, V5, V6), or more explicitly

The Jacobian matrices of F and V evaluated at the solution E0∗ = (0, 0, 0, 0, 0, A/μ) are given by

Direct calculations show that where

Denoting the 4 × 4 identity matrix by ℐ, the characteristic polynomial Γ(λ) of the matrix FV−1 is given by where

The solutions λ1,2,3,4 are given by

Therefore, the reproduction number R0 for the SEYNHR model (6) is given by

For simplicity, let us denote

In order to interpret the basic reproduction number R0, we split the above expression of R0 as follows: where

The first and second terms in R 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 R 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 R 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.

Theorem 1 .

The DFE E0 of the system (6) is locally asymptotically stable if R0 < 1 and unstable otherwise.

The Jacobian matrix J 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: where where ψ = (β1μ + (r + μ2)β1 + β3ε) and C1 is obviously positive. We can rewrite R0 as

Clearly if R0 < 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 DFE E0 of the COVID-19 transmission model is studied in the following result. where W 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

Theorem 2 .

The DFE E0 of the system (6) is globally asymptotically stable if R0 < 1 and unstable otherwise.

Let us consider the following Lyapunov function, described in [9]:

Now choosing

We obtain

Hence, if R0 < 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:

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]

The birth rate A is obtained from A/μ = Q(0), where the total Saudi population Q(0) is 34218169; therefore, A = 1252 per day [11]

Mean asymptomatic infectious period Φ is assumed to be the same as the mean asymptomatic infectious period Φ 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

Table 1

Parameters used in the simulations.

Parameter	Description	Dimension	Value	A 95% confidence interval	Source
N	Population of Saudi Arabia in 2019.	Individual	34218169	—	[11]
A	Birth rate in Saudi Arabia in 2019.	Individual/day	1252	—	[11]
μ	Death rate in Saudi Arabia in 2019.	Day −1	3.6593x10−5	—	[11]
β 1	Effective contact rate from symptomatic to susceptible.	Day −1	0.35	[0.267, 0.5]	Fitted
β 2	Effective contact rate from asymptomatic to susceptible.	Day −1	0.7	[0.53, 1]	Fitted
β 3	Effective contact rate from hospitalized to susceptible.	Day −1	0.18	[0.13, 0.26]	Fitted
α −1	Mean latent period.	Days	5.1	—	[23]
γ	Probability of becoming asymptomatic upon infection.	n/a	0.86	—	[24]
Φ Y −1	Mean symptomatic infectious period.	Days	8	—	[25]
ε	Rate of symptomatic becomes hospitalized.	Day −1	0.125	—	[1]
μ 1	Death rate of symptomatic patients.	Day −1	0.392	—	[26]
K	Rate of asymptomatic becomes symptomatic.	Day −1	0.12	[0.1, 0.125]	Fitted
Φ N −1	Mean asymptomatic infectious period.	Days	8	—	Assumed
r	Rate of recovered hospitalized patients.	Day −1	0.1	—	[6]
μ 2	Death rate of hospitalized patients.	Day −1	0.392	—	[1]

The function Z 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 and the solution of the model y 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:

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 the SEYNHR 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.

5. Numerical Results and Discussion

To illustrate the numerical results, we take the initial population and subpopulations as Q(0) = 34218169, Y(0) = 2, S(0) = Q(0) − Y(0) = 34218167, and E(0) = N(0) = H(0) = R(0). 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 Figure 3(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 Y(t) 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 Y(t) 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 Y(t) 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 Y(t). 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 Y(t). 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 Y(t). The baseline value of β3 is as given in Table 1.

Figure 7

Impact of the rate of asymptomatic becomes symptomatic K on the number of new COVID-19-infected cases Y(t).

Based on estimated and measured parameter values, the basic reproduction number is estimated as R0 ≈ 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 of R0 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 Φ−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 Section 3 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 of R0, 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.