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 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 those 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. 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-7]. The most widely used models, and among the simplest ones, proved to be those proposed by Kermack and McKendrick [8-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 individuals 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-14] with more reference therein, trying to capture the changes in the dynamics of the model interaction, particularly with respect to the reproductive number . 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–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-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 diagnosed- and 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 (COVID-19) 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 and high-risk susceptible individuals; the compartment of exposed individuals E; the compartment of infected individuals (divided into two compartments, namely that of infected undiagnosed individuals and the one consisting of infected diagnosed individuals ); the compartment of recovered undiagnosed individuals , that of recovered diagnosed 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 , 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 . 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 , would be formulated from Fig. 2 as follows: 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. represents the number of infected, undiagnosed individuals, asymptomatic, and symptomatic individuals with mild symptoms that have not been identified by the authorities. 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. represents the number of most susceptible individuals, with higher risk of contracting the COVID-19 infection. represents the number of individuals that are less susceptible to the COVID-19 infection. represents the number of recovered, undiagnosed individuals, who are not being officially identified as such. 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 represents the reduction risk factor of the COVID-19 infection for the less susceptible individuals . 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 The time episode of the lockdown imposed by the authorities is denoted by an increasing positive sequence , and the authorities-controlled parameters are given by , . The positive initial conditions are given by and , where 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 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 where 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).

Figure 1

Graphic of V. Altounian [26], showing that masks reduce airborne transmission

Figure 2

Diagram of the schematic representation of the flow of individuals among different stages of the infection. The mathematical model is called EIISSRREx-model and considers two distinct susceptible compartments: the most susceptible individuals and the less susceptible individuals. Also, there are uninfected I, infected individuals, asymptomatic, and symptomatic individuals. stands for the infected individuals and there are asymptomatic/undiagnosed (undetected) and those who will recover ; represents those diagnosed, hospitalized individuals that will recover; E represents the group of those diagnosed hospitalized individuals who die

Table 1

Parameters used in model (1–8)

Parameters	Meaning	Values	Ref.
α_E	Contact disease rate of a person in compartment E at time t	Estimated day^−1	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	—	Table 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⋅106	—

The positive initial conditions , , and are given in . The transition rate 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 to denote the space of all continuous functions from into . It is a well-established mathematical fact [33] that the autonomous system of differential equations (4)–(8), which depends on the initial conditions , , , , , , , , is uniquely solvable.

Theorem 3.1

There exists a unique solution to the system of equations (1)–(8) in . Furthermore, the solution is positive and bounded from above.

Proof

The right-hand side of each of equations (1)–(8) is continuously differentiable almost everywhere in . 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 .

Denote the solution of (4) by . If for any then we would have . From equation (4), we get which is a contradiction. Therefore, for all . In a similar way, we can prove that is positive for all .

Denote the solution of equation (1) by . Assuming the existence of such that such that , we would have . From equation (1), we have and it follows that Using equations (2) and (3), we obtain Then we have where and Then, for any positive initial condition, and are positive. Again, equation (10) yields This is a contradiction and thus it must hold that for all .

It can be shown analogously that , , , , and are positive for all , hence the positiveness of the solution of system (1)–(8).

Next, we show the boundedness of , , , , , , , and for all . We tackle first the function . Equation (4) yields, for , and from equation (5)

Thus, and 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 , 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.

Since and are a finite, using the positivity and boundedness of the operator in integrals (11) and (12), we conclude that Again, since E, , and are positive and bounded, and from equations (2)–(3) and (6)–(8), we have Then Now, we will show that . Suppose that . Then, by using the fluctuation lemma, see for example Hirsch [34], there exists a sequence such that, for all , It is clear from the model equation (1) that we have then where , , and . Therefore, as , we obtain If we have , then we get a contradiction as we have If we have , the model equation (1) implies that as , and we have , which will be a contradiction. Finally, if we suppose , then 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 and the infected compartment by . 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 , which depend on the infected compartment as follows: where To find the basic reproduction number , the average number of individuals that one single infected individual is capable of infecting during each time episode of the lockdown in the interval , , 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 with 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 of new infections in each compartment i and the rate of transfer (resp. of individuals into (resp. out of) in each compartment i of the model equations can be written as follows: and where . The model equations given in (1)–(8) can be rewritten as the following system of equations: where . To discuss the local asymptotic stability of the endemic steady state , we linearize around the equilibrium . Using the notation for the corresponding linearized model around , we have then where To understand the qualitative behavior of the model equations (1)–(8), we observe that the infected compartments E, , in the equilibrium are zero and that the remaining compartments , , , , Ex are at the given equilibrium , , , , .

Now, if we set for all , where includes the infected subjects with susceptible individuals for the infection, it is clear that the new Jacobian matrix has three negative eigenvalues given by , , and that the characteristic polynomial can be written as where In the proof of the local asymptotic stability, we borrow the powerful tools from the theory of M-matrices [35-39]. For the reader’s convenience, we state the characterization of M-matrices with regard to the spectral radius (SR) and to the spectral abscissa (SA).

Recall that given a matrix A, its spectral radius, denoted by , is given by Clearly, an M-matrix A is nonsingular if and only if the following condition holds: for some in (SR), so that in (SA).

If can be written in the form , where and is nonnegative (, if , , ) and , then A is called an M-matrix.

A matrix A is said to be an M-matrix if has the pattern, that is, if A belongs to the class and if each eigenvalue of A satisfies

The following lemma will be useful throughout the characterization of the reproduction number . We refer the reader to [39, p. 159] and [40, p. 127] for further details.

Lemma 3.3

([39, 40])

Let V be a nonsingular M-matrix such that both J and have the sign pattern. Then J is a nonsingular M-matrix if and only if is a nonsingular M-matrix.

Now, conditions – in [41] are satisfied for the model equations (13)–(16). Therefore the rate of appearance of new infections in each compartment i, , and the rate of transfer of individuals into (resp. out of) each compartment i, (resp. , of the model equations (13)–(16) satisfy where and Following the same arguments as those of Lemma 1 in van Driessche and Watmough [41], we have where and where and It is clear that all diagonal elements of V̂ are positive, whence V̂ is a nonsingular M-matrix, see [35-39] and F̂ is a nonnegative matrix. As in [2, 42], we refer to as the next generation matrix of the model equations (13)–(16) and we define the basic reproduction number as follows: where denotes the spectral radius of the matrix A given in (17).

The free infected compartments including , are monotonically decreasing, whereas the compartments of recovered and defunct individuals , , Ex are monotonically increasing and converge to their asymptotic equilibrium if and only if the compartments of infected individuals 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.

Theorem 3.4

The endemic equilibrium of the model equations (1)–(8) exists and is locally asymptotically stable during each time episode of the lockdown interval , where is an increasing positive sequence if and only if , where is the basic reproduction number given in (19). In other words, if and only if where

To study the local asymptotic stability of (1)–(8) around the equilibrium in each time episode of the lockdown interval , we will focus first on the Jacobian subsystems matrices F̂ and V̂ including the compartments E, , , , around the equilibrium 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 a nonsingular M-matrix and is nonnegative, it is readily concluded that is nonnegative. Let us consider the following matrix difference: It is clear that and have the sign pattern. Then, by using Lemma 3.3, we conclude that is a nonsingular M-matrix if and only if is a nonsingular M-matrix. The spectral radius (SR) characterization of nonsingular M-matrices guarantees that is a nonsingular M-matrix if and only if . Again, using (SA) and the basic reproduction number given in (19), it follows that if and only if . □

Now, we can define the effective reproduction number as the number of cases that one single infected individual is capable of infecting during different time episodes of the lockdown interval , . In other words, we are going to protect uninfected individuals taking into account the efficiency of the control measures imposed by the authorities proposed in (9). Furthermore, where 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 represents the data measurement available at the time , N stands for the number of measurements, and represents the (COVID-19) model equations (1)–(8). The best fit relative estimation of the control parameters , , , , , and the estimated parameters , , 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.

Table 2

Estimated parameters of model (1)–(8), taking into account different control stage measures

Parameters	0 to 30 days	30 to 90 days	90 to 120 days	120 to 260 days.
α_E	0.970103	0.745319	25.285168	60.624593.
α_I_U	0.244387	0.000002	0.000688	0.000007.
α_I_D	0.014381	0.000061	0.000091	0.000131.
β_E	0.1818	0.1818	0.190241	0.076462.
β_I_U	0.003330	0.000645	0.000464	0.000958.
β_R_U	0.0752	0.0752	0.0352	0.0352.
β_R_D	0.18554	0.779554	0.726325	1.1825.
β_Ex	0.010549	0.006583539	0.009839	0.02299.
μ_L	0.05	0.05	0.01	0.00333333.

Figure 3

Model prediction taking into account different stages of the control measures imposed by the authorities and the corresponding epidemiological data, including the number of individuals diagnosed with (COVID-19), recovered and dead, collected from the Saudi Arabian Ministry of Health, see “https://covid19.moh.gov.sa/”

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 , , and 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 given in (20), it is clear that if we reduce the number of the more susceptible individuals and the parameter 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.