A mathematical model for the spread of COVID-19 and control mechanisms in Saudi Arabia

In this work, we develop and analyze a nonautonomous mathematical model for the spread of the new corona-virus disease (COVID-19) in Saudi Arabia. The model includes eight time-dependent compartments: the dynamics of low-risk \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$S_{L}$\end{document}SL and high-risk \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$S_{M}$\end{document}SM susceptible individuals; the compartment of exposed individuals E; the compartment of infected individuals (divided into two compartments, namely those of infected undiagnosed individuals \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$I_{U}$\end{document}IU and the one consisting of infected diagnosed individuals \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$I_{D}$\end{document}ID); the compartment of recovered undiagnosed individuals \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{U}$\end{document}RU, that of recovered diagnosed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$R_{D}$\end{document}RD individuals, and the compartment of extinct Ex individuals. We investigate the persistence and the local stability including the reproduction number of the model, taking into account the control measures imposed by the authorities. We perform a parameter estimation over a short period of the total duration of the pandemic based on the COVID-19 epidemiological data, including the number of infected, recovered, and extinct individuals, in different time episodes of the COVID-19 spread.


Introduction
In the last century several compartmental models in epidemiology are derivatives of the type Susceptible, Infected, and Recovered individuals (SIR), see for example [1,2] for further details and recent modeling contributions with more references therein [3][4][5][6][7]. The most widely used models, and among the simplest ones, proved to be those proposed by Kermack and McKendrick [8][9][10]. One of the shortcomings of this work is the oversimplification of two different effects: (i) susceptible; and (ii) infectious, which are both lumped into a single effect. This kind of simplification does not take into account an individual's degree of exposure to the infection. Hospital staff members or individuals with weakened immune systems or with chronic diseases, for example, should be placed in a category different from that of the individuals who, due to their occupation, age, or physical condition are less susceptible to contracting the infection. In particular, we quantitatively and qualitatively distinguish these two groups by considering the category of susceptible in-dividuals as the union of two different compartments: that of higher-risk individuals and the compartment of lesser-risk individuals.
Extensions of SIR-type models have been proposed by several authors, see for example [11][12][13][14] with more reference therein, trying to capture the changes in the dynamics of the model interaction, particularly with respect to the reproductive number R 0 . An overview of some of the mathematical models in epidemiology which appeared in the literature is given in [15], enhanced with a recent publication on COVID-19 modeling a series of the parameter values (rates) given in [3,[16][17][18]. Such parameters will provide good guidelines for our parameter estimation procedure.
Given the complexity of the biological systems and the interactions between given subsystems, where we have data limitations, mathematical modeling and analysis are needed to quantify such iterations and understand the behavior of such subsystems. The theoretical result such as the existence of positive solutions and their boundedness are important by nature, such as the number of individuals or the biomass dependent on the given model, where they should be positive and bounded. In addition to studying the stability of such a system of equations, the models must be as simple as possible and complex enough as necessary to reflect the reality, see for example [3,19]. Consequently, real life data collection is very important for the parameter estimations in order to understand the interaction between the given subsystems. A broad overview of model validation is presented in [19][20][21], where more details about parameter estimations are given.
Recently, Lin et al. in [22] proposed a compartmental model with time-dependent transmission rate, which incorporates the impact of governmental action that could be modeled by step functions. On the other hand, Lin's model should be in correspondence with the understanding of the quantity and quality of the available data. In addition to dividing the susceptible individual into higher-and lesser-risk compartments, we suggest splitting the compartment of infected individuals into the subcompartments consisting of diagnosedand undiagnosed-infected individuals. This method will be helpful in controlling the outbreak of COVID-19 (see the model proposed by Chowell et al. [23], which concerns the SARS outbreak in Ontario, Hon Kong, and Singapore, where they show that diagnosing and isolation help control the spread of the virus). Due to the large scale of the model, we conduct a careful sensitivity analysis to determine which model parameters have a strong influence on the dynamics of the system and therefore can be determined accurately in the parameter estimations, as was demonstrated in [19,20].
In Sect. 2, we develop a nonautonomous mathematical model for the spread of the new coronavirus disease  in Saudi Arabia, taking into account the control measures imposed by the authorities, which are time dependent. The model includes eight time-dependent compartments: the dynamics of low-risk S L and high-risk S M susceptible individuals; the compartment of exposed individuals E; the compartment of infected individuals (divided into two compartments, namely that of infected undiagnosed individuals I U and the one consisting of infected diagnosed individuals I D ); the compartment of recovered undiagnosed individuals R U , that of recovered diagnosed R D individuals, and the compartment of extinct Ex individuals. In Sect. 3, we investigate the qualitative aspect of the model including the model persistence and local stability, and Sect. 4 is devoted to the quantitative dynamics of the spread and the parameter estimation over a short period of the total duration of the pandemic based on the COVID-19 data, including the number of infected, recovered, and extinct individuals in different episodes of the COVID-19 spread. Finally, concluding remarks and future directions of our research are offered in Sect. 5.

Model development
As mentioned above, the compartment of susceptible individuals is divided into two further subcompartments: that consisting of the most susceptible individuals, denoted by S M , with higher risk of contracting the COVID-19 disease and the group of individuals with lower likelihood of getting infected. The latter will be called the compartment of less susceptible individuals and will be denoted by S L . These groups of individuals will be considered according to how closely they follow the sanitary recommendations of the authorities in their daily life, namely, whether they often wash their hands, regularly disinfect their surroundings, and wear face masks. In general, susceptible individual are not infected by the disease-causing pathogen. These individuals are socially active and within certain limits follow the sanitary protection measures suggested by the US Center for Disease Control and Prevention (CDC), such as social distancing of 1.5-2 m, regular hand washing, and wearing of face masks in their daily activities, in order to limit the airborne transmission of COVID-19, see [24]. Further measures might be needed for reducing airborne transmission of COVID-19, such as room ventilation and regular disinfection, which lead to effectively limiting the concentration of SARS-CoV-2 in aerosols, see [25]. In practice, however, most of the time social distancing is not respected. Yet, the sanitary protection provided by the universal wearing of masks represents a promising practice for reducing the transmission of the COVID-19 infection, as was demonstrated in [26]. The epidemiological data provided at the beginning of the pandemic by a number countries, such as Taiwan, Japan, Hong Kong, Singapore, and South Korea, show the efficacy of universal masking as a control measure, even in the absence of a severe lockdown during the pandemic. In [26] the necessity of masking and testing in the fight against asymptomatic spread in aerosols and droplets is shown, and evidence is provided of the fact that no masking maximizes exposure, whereas the least level of exposure is achieved through universal masking. The universal use of face masks, it is shown in the work in point, helps reduce the size of the spread particles from 100 μm to 1-0.1 μm, especially in asymptomatic people and those with mild symptoms, see [27] and Fig. 1. Prather et al. [26] demonstrate that the aerosol filtering efficiency of different materials, thickness, and layers used in properly fitted homemade masks was similar to that of medical masks, see [28]. Tellier et al. [29] show that in still air, a 100 μm droplet will drop to the ground from 2.4384 m in 4.6 s, whereas a 1 μm aerosol particle will take 12.4 hours. The (CDC) recommends 1.8288 m for social distancing, but this recommendation might not be sufficient for indoor conditions, where aerosols can remain airborne and accumulate for hours, see [30,31]. The authors in [31] show that breezes and winds often transport infected droplets and aerosols long distances. Further experimental research is needed in order to more precisely characterize the transport and to achieve a better understanding of the relevance of airborne transmission of the COVID-19 infection [26]. For the simulation of COVID-19, the epidemiological EIISSRREx-model, for all t ≥ t 0 , would be formulated from Fig. 2 as follows: Figure 1 Graphic of V. Altounian [26], showing that masks reduce airborne transmission and those who will recover R U ; R D represents those diagnosed, hospitalized individuals that will recover; E represents the group of those diagnosed hospitalized individuals who die Here, E represents the number of individuals in the incubation period, that is, the number of persons exposed to the virus with no visible symptoms. I U represents the number of infected, undiagnosed individuals, asymptomatic, and symptomatic individuals with mild symptoms that have not been identified by the authorities. I D represents the number of infected, diagnosed individuals, that is, the number of persons that have been officially identified as infected and are either hospitalized or in quarantine at home. Those individuals are registered in an official applications provided by the local Saudi Authority, called "Tawakkalna ( )", see "https://ta.sdaia.gov.sa/en/index". Such application will help with contact tracing. Therefore as we increase the number of testing, we will be able to localize clustering of infected individuals which will help to reduce their numbers. S M represents the number of most susceptible individuals, with higher risk of contracting the COVID-19 infection. S L represents the number of individuals that are less susceptible to the COVID-19 infection. R U represents the number of recovered, undiagnosed individuals, who are not being officially identified as such. R D represents the number of infected, recovered individuals officially identified as such and finally Ex represents the number of individuals who are deceased from COVID-19 (called extinct individuals). The parameter μ L represents the reduction risk factor of the COVID-19 infection for the less susceptible individuals S L . C M (t) represents the efficiency of the control measures imposed by the authorities, such as the lockdowns. This quantity can be modeled as a step function in the following fashion: where  Table 2 α I U Contact disease rate of a person in compartment I U at time t Estimated day -1 Table 2 α I D Contact disease rate of a person in compartment I Diag at time t Estimated day -1 Table 2 β E Transition rate of a person in compartment E at time t Estimated day -1 Table 2 β I U Transition rate of a person in compartment I U at time t Estimated day -1 Table 2 β R U Rate at which an undiagnosed infected person recovers at time t Estimated day -1 Table 2 β R D Rate at which a diagnosed infected person recovers at time t Estimated day -1 Table 2 β Ex Rate at which a diagnosed infected person dies at time t Estimated day -1 Table 2 μ L Reduction risk factor of infection in compartment S L at time t -T a b l e 2 ρ Proportion of the population size N that is initially at higher risk of contracting the infection 0.4 [23] N Total size of the population 30 · 10 6 - The time episode of the lockdown imposed by the authorities is denoted by an increasing positive sequence (t i ) 0≤i≤n+1 , and the authorities-controlled parameters are given by κ i ∈ (0, 1), i = 0, . . . , n. The positive initial conditions are given by where t 0 ≥ 0 represents any initial time and ρ is the proportion of the population size N that is initially at higher risk of contracting the infection, see [23]. The disease constant rates α E , α I , α I D , α I U ∈ (day -1 ) of individuals in the corresponding compartment are supposed to be constant, without taking into account the control measures. Without loosing any generality, we can suppose that all parameters are piecewise constants, as in (9), which means represents the time-dependent model parameters, and are the model parameters constant over each time episode of the lockdown at the interval Table 1 displays a detailed description of the parameters (including their units).
The positive initial conditions E(t 0 ) = E 0 , I D (t 0 ) = I 0 D , and I U (t 0 ) = I 0 U are given in [0, ∞). The transition rate β E from the exposed compartment E to the infected compartment I is supposed to be constant at the beginning of the pandemic according to [32].

Model analysis
The model provides a description of the mass transfer property between the compartments. It follows that each solution of the equations involved in it (1)-(8) should be positive, bounded, and persistent. The sum of the states (total population) is constant and equal to N , i.e., In the sequel, we set C(R + , R 8 ) to denote the space of all continuous functions from R + into R 8 . It is a well-established mathematical fact [33] that the autonomous system of differential equations (4)- (8), which depends on the initial conditions  8) is continuously differentiable almost everywhere in C(R + , R 8 ). We refer the reader to [33] for the mathematical result that guarantees that under such circumstances the system of equations (1)-(8) with the given initial conditions is uniquely solvable for all t ≥ 0. Denote the solution of (4) by S M (t). If for any which is a contradiction. Therefore, S M (t) ≥ 0 for all t ≥ t 0 . In a similar way, we can prove that S L (t) is positive for all t ≥ 0. Denote the solution of equation (1) by E(t). Assuming the existence of τ 1 such that such that E(τ 1 ) = 0, we would have dE(t) dt |τ 1 < 0. From equation (1), we have and it follows that Using equations (2) and (3), we obtain d dt Then we have Then, for any positive initial condition, I U (τ 1 ) and I D (τ 1 ) are positive. Again, equation (10) yields This is a contradiction and thus it must hold that E(t) ≥ 0 for all t ≥ t 0 . It can be shown analogously that I U (t), I D (t), R U (t), R D (t), and Ex(t) are positive for all t ≥ 0, hence the positiveness of the solution of system (1)- (8).
Next, we show the boundedness of S M (t), S L (t), E(t), I U (t), I D (t), R U (t), R D (t), and Ex(t) for all t > 0. We tackle first the function S M (t). Equation (4) yields, for t ≥ t 0 , and from equation (5) S L (t) = e Thus, S M (t) and S L (t) are bounded from above. Now, we show that any given solution of (1)-(8) is bounded. It is clear from the model equations (4)-(8) that we have and by using the positivity of each state equation, we can conclude that, for all t ≥ 0, all solutions are bounded from above.
Throughout this paper, we use the following notations for the lim sup and the lim inf of output of the model in each compartment:

Theorem 3.2 Model (1)-(8) is persistent. This is to say that the components of each solution are eventually uniformly bounded from above and away from zero.
Proof Since S 0 M and S 0 L are a finite, using the positivity and boundedness of the operator C M (t)α E E(t) + C M (t)α I U I U (t) + C M (t)α I D I D (t) in integrals (11) and (12), we conclude that 0 < S M < S M < ∞ and 0 < S L < S L < ∞.

Again, since E, I U , and I D are positive and bounded, and from equations (2)-(3) and (6)-(8), we have
Then Now, we will show that E = lim inf t→∞ E(t) > 0. Suppose that E < E. Then, by using the fluctuation lemma, see for example Hirsch [34], there exists a sequence {t k } ∞ k=1 such that, for all k ≥ 1, It is clear from the model equation (1) that we have where m = min t∈[t 0 ,∞) (t) > 0, := C M , S M , S L , I U , I D , μ L , and β E = max t∈[t 0 ,∞) β E (t). Therefore, as k → ∞, we obtain If we have E = 0, then we get a contradiction as we have If we have lim t→∞ E = 0, the model equation (1) implies that dE(t) dt > γ > 0 as t → ∞, and we have lim t→∞ E = ∞, which will be a contradiction. Finally, if we suppose E = E, then lim t→∞ E exists and E is eventually uniformly bounded from above and away from zero, which is the desired conclusion.
It is clear that the last three equations (6)-(8) can be solved explicitly and that the solutions will depend on the infected compartments (6)- (7). We denote the susceptible compartment by x S = (S M , S L ) T and the infected compartment by x I = (E, I U , I D ) T . Then the model equations (1)-(5), can be written as a feedback control structure and the remaining recovered and extinct compartments (8)-(6) can be denoted as x REx = (R U , R D , Ex) T , which depend on the infected compartment x I as follows: where To find the basic reproduction number R 0 , the average number of individuals that one single infected individual is capable of infecting during each time episode of the lockdown in the interval [t i , t i+1 ], 0 ≤ i ≤ n, should be kept in mind. In other words, we are going to protect uninfected individuals as we have no control measure. In this case, the number of uninfected individuals will depend on the size of the population and not on the time period of the spread. Therefore, to study the stability, and for the sake of simplicity, we will suppose that all parameters are constant. The possible equilibrium of the system of equations (1)-(8) is given by Therefore only most susceptible, low susceptible, received, and deceased individuals are present, which means that the epidemic is over. It is clear that the rate of appearance F = (F i ) 1≤i≤8 of new infections in each compartment i and the rate of transfer V + = (V + i ) 1≤i≤8 (resp. V + = (Vi ) 1≤i≤8 of individuals into (resp. out of ) in each compartment i of the model equations can be written as follows: The model equations given in (1)-(8) can be rewritten as the following system of equations: To discuss the local asymptotic stability of the endemic steady state X * = (0, 0, 0, S M , S L , R U , R D , Ex) T , we linearize around the equilibrium X * . Using the notation for the corresponding linearized model around X = (E l , I l U , I l D , S l M , S l L , R l U , R l D , Ex l ) T , we have then To understand the qualitative behavior of the model equations (1) x I x S ), where x includes the infected subjects with susceptible individuals for the infection, it is clear that the new Jacobian matrix J has three negative eigenvalues given by -β E , -(β I U + β R U ), -(β R D + β Ex ) and that the characteristic polynomial can be written as In the proof of the local asymptotic stability, we borrow the powerful tools from the theory of M-matrices [35][36][37][38][39]. For the reader's convenience, we state the characterization of Mmatrices with regard to the spectral radius (SR) and to the spectral abscissa (SA).
Recall that given a matrix A, its spectral radius, denoted by ρ(A), is given by (SR) [38,39]; If A ∈ R n×n can be written in the form [38,39]; A matrix A is said to be an M-matrix if A ∈ R n×n has the Z n×n pattern, that is, if A belongs to the class Clearly, an M-matrix A is nonsingular if and only if the following condition holds: s = ρ(B) for some B ≥ 0 in (SR), so that (λ i (A)) > 0 in (SA).
The following lemma will be useful throughout the characterization of the reproduction number R 0 . We refer the reader to [39, p. 159] and [40, p. 127] for further details.  [41] are satisfied for the model equations (13)- (16). Therefore the rate of appearance of new infections in each compartment i, F = ( F i ) 1≤i≤5 , and the rate of transfer of individuals into (resp. out of ) each compartment i, V + = ( V + i ) 1≤i≤5 (resp. V + = ( Vi ) 1≤i≤5 , of the model equations (13)- (16) satisfy Following the same arguments as those of Lemma 1 in van Driessche and Watmough [41], we have and It is clear that all diagonal elements of V are positive, whence V is a nonsingular M-matrix, see [35][36][37][38][39] and F is a nonnegative matrix. As in [2,42], we refer to F V -1 as the next generation matrix of the model equations (13)- (16) and we define the basic reproduction number R 0 as follows: where ρ(A) denotes the spectral radius of the matrix A given in (17). The free infected compartments including S M , S L are monotonically decreasing, whereas the compartments of recovered and defunct individuals R U , R D , Ex are monotonically increasing and converge to their asymptotic equilibrium if and only if the compartments of infected individuals (E, I U , I D ) T converge to the zero equilibrium, taking into account the time varying of the compartments of susceptible individuals, and converge to the steady state value. Thus we have the following theorem about the asymptotic stability of (13)- (16). Now we are ready to state the main theorem concerning the local asymptotic stability of the endemic equilibrium.  (19). In other words, if and only if Proof To study the local asymptotic stability of (1)-(8) around the equilibrium X * in each time episode of the lockdown interval [t i , t i+1 ], we will focus first on the Jacobian subsystems matrices F and V including the compartments E, I U , I D , S M , S L around the equilibrium (0, 0, 0, S M , S L ) T of the model equations (13)- (14). By using the inverse-positivity characterization of nonsingular M-matrices given in [39] and taking into account that V , given in (18) is nonnegative, it is readily concluded that F V -1 is nonnegative. Let us consider the following matrix difference: It is clear that -J 1 and -J 1 V -1 = I -F V -1 have the Z sign pattern. Then, by using Lemma 3.3, we conclude that -J 1 is a nonsingular M-matrix if and only if I -F V -1 is a nonsingular M-matrix. The spectral radius (SR) characterization of nonsingular M-matrices guarantees that I -F V -1 is a nonsingular M-matrix if and only if ρ( F V -1 ) < 1. Again, using (SA) and the basic reproduction number R 0 given in (19), it follows that (λ i (J 1 )) < 0 if and only if R 0 = ρ( F V -1 ) < 1.
Now, we can define the effective reproduction number R e as the number of cases that one single infected individual is capable of infecting during different time episodes of the lockdown interval [t i , t i+1 ], 0 ≤ i ≤ n. In other words, we are going to protect uninfected individuals taking into account the efficiency of the control measures C M imposed by the authorities proposed in (9). Furthermore, where R 0 is the basic reproduction number given in (19).

Data collection and parameter estimations
Looking at the response to the COVID-19 spread in Saudi Arabia, we are going to present different stages of the actions taken by the authorities, taking into account the intensity of the restrictive measures in different time episodes. In the interest of simplicity and considering various control measures within the country, we will divide our time period in four different episodes, namely, from day 0 to day 30, from day 30 to day 90, from day 90 to day 120, and from day 120 to day 260. The data were collected from the Saudi Arabian Ministry of Health, see "https://covid19.moh.gov.sa/", and present the number of diagnosed COVID-19 infections, recovered and extinct individuals every day over a 9-month period. These data will be used to estimate the parameters over a given period of the duration of the restriction. It will be assumed that all parameters of our model are given in Table 1. We conduct a sensitivity analysis to determine which model parameters have a strong influence on the dynamics of the model equations (1)- (8), see [19,20] for more details, on the parameter estimates methodology using differential equations. We are going to estimate the control parameters and the parameters of recovered and the dead compartments over different time episodes by using multi-objective optimization. We will perform several fits for each time episode using the Nelder-Mead search algorithm to find a local minimum. We set  where d j represents the data measurement available at the time t j , N stands for the number of measurements, and x(t j , θ ) represents the (COVID-19) model equations (1)- (8). The best fit relative estimation of the control parameters α E , α I U , α I D , β I U , β E , and the estimated parameters β R U , β R D , β Ex of recovered, and the dead compartments are given in Table 2, with different time episodes. The model estimations compared to the given data of the number of diagnosed COVID-19 infections, recovered and extinct individuals every day over a 9-month period are shown in Fig. 3.

Concluding remarks and suggestions
It is necessary to achieve a better understanding of the COVID-19 dynamics, taking into account different control measures in order to reduce the number of infections and to decrease the mortality rate. The model proposed in the present manuscript has the advantage of describing the best way of controlling the pandemic. As was demonstrated in the parameter estimation, the ability to reduce the parameters α I U , α I D , and β I U in this case will result in the detection of the number of infected individuals and in the reduction of the number of undiagnosed individuals by enhancing the testing. This will facilitate the control of the virus. Again, considering the effective reproduction number R e (t) given in (20), it is clear that if we reduce the number of the more susceptible individuals and the parameter μ L by respecting the social distancing and wearing masks, it will help with the control of the pandemic. We aim at improving the model by including the number of individuals in intensive care and by adding a modification that allows for the effect of a vaccine, as given in [43] for measles epidemics. Another important effect that should be included will be the incubation time delays of receiving COVID-19 positive test from the authority.
The study of persistence and stability in conjunction with the reproduction number is of paramount importance for the parameter estimation. We believe that this model provides a good basis for the study of more general cases, in which these parameters vary over each time interval. It will be of interest to develop a parallel algorithm, see [44,45], in order to reduce the time needed to find the best fit of the multi-objective cost functional.