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


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][3][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][4][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 R 0 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.

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 Φ 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 .
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 R 0 , for the compartmental model above. Calculating R 0 is a useful metric for assessing the transmission potential of an emerging COVID-19 in Saudi Arabia.

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 R 0 is established.

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.  (6) is nonnegative if they exist.
Proof. Assume t 1 = ft > 0 : PðtÞ > 0 ∈ ½0, tg, 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.
Proof. 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, lim t⟶∞ sup QðtÞ ≤ A/μ. Furthermore, ðdQðtÞÞ/d t < 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. 2 Computational and Mathematical Methods in Medicine 3.2. Basic Reproduction Number R 0 . It can be easily verified that system (6) always has a disease-free equilibrium (DFE), which we will denote by E 0 , 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 R 0 ; more details about the method are in [8]. The system (6) can be rewritten as follows: where F ≔ ðF1, F2, F3, F4, F5, F6Þ T and V ≔ ðV1, V2, V3, V4, V5, V6Þ T , or more explicitly The Jacobian matrices of F and V evaluated at the solution E * 0 = ð0, 0, 0, 0, 0, A/μÞ T are given by Direct calculations show that where Denoting the 4 × 4 identity matrix by I , 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 R 0 for the SEYNHR model (6) is given by

Computational and Mathematical Methods in Medicine
For simplicity, let us denote In order to interpret the basic reproduction number R 0 , we split the above expression of R 0 as follows: where The first and second terms in R YH 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 EYNH1 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 EYNH2 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 R 0 shows that asymptomatic as well as symptomatic cases of COVID-19 could drive the growth of the pandemic in Saudi Arabia.

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 R 0 < 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 R 0 < 1.
Theorem 3. The DFE E 0 of the system (6) is locally asymptotically stable if R 0 < 1 and unstable otherwise.
Proof. The Jacobian matrix J E 0 obtained at the DFE E 0 is as follows: Clearly, the eigenvalues −μ and −μ are negative real numbers. The remaining eigenvalues can be obtained through the following equation: where ψ = ðβ 1 μ + ðr + μ 2 Þβ 1 + β 3 εÞ and C 1 is obviously positive. We can rewrite R 0 as Clearly if R 0 < 1, then C 2 is positive. Similarly, we can show that C 3 and C 4 are all positive if R 0 < 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 R 0 < 1.
The global stability of DFE E 0 of the COVID-19 transmission model is studied in the following result.
where W i 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 4. The DFE E 0 of the system (6) is globally asymptotically stable if R 0 < 1 and unstable otherwise.
Proof. Let us consider the following Lyapunov function, described in [9]: Now choosing We obtain Hence, if R 0 < 1, then ΠðtÞ/dt < 0. Therefore, the largest compact invariant set in Ω is the singleton set E 0 , and using LaSalle's invariant principle [10], E 0 is globally asymptotically stable in Ω.

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 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 Φ 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 leastsquare 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 The function Z is dependent on t, the vectors of the dependent variables u, and unknown parameters θ to be The model (6) is solved using MATLAB solver ode45 which uses the Runge-Kutta methods to solve the initial value problem. Then, the lsqcurvef it 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.
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-19infected 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-  Table 1.  Table 1.

Computational and Mathematical Methods in Medicine
Number of symptomatic infected individuals Y (t) 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. 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. Based on estimated and measured parameter values, the basic reproduction number is estimated as R 0 ≈ 2:7. The variation of the basic reproduction number R 0 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 R 0 decreases and vice versa. The lower heat map of Figure 8 shows that increasing K increases R 0 and vice versa. In reality, R 0 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 R 0 vary greatly: from 1 to 6 [13,[16][17][18][19][20][21][22] up to 26.5 [12]. This variation is because of the different assumptions and factors they had considered in their calculations.
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 N . 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).

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 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 to the real data of infected cases. An analytical expression of the basic reproduction number R 0 is obtained in Section 3, and the numerical value is estimated as R 0 ≈ 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 R 0 , 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.

Data Availability
The data used in my paper is the total number of confirmed cases of COVID-19 reported by the Ministry of Health in Saudi Arabia from the 2nd of March until the 21st of April 2020. Data supporting this paper are open and can be accessed through (in Arabic) Ministry of Health COVID 19 Dashboard. Saudi Arabia; 2020 (https://covid19.moh.gov .sa). Alternatively, data supporting this paper are open and can be accessed through (in English) Our World in Data; Oxford; 2020 (https://covid.ourworldindata.org/data/owidcovid-data.csv).

Conflicts of Interest
The author declares that he has no competing interests.