key: cord-0289433-57vaiail
authors: Mahrouf, Marouane; Boukhouima, Adnane; Zine, Houssine; Lotfi, El Mehdi; Torres, Delfim F. M.; Yousfi, Noura
title: Modeling and Forecasting of COVID-19 Spreading by Delayed Stochastic Differential Equations
date: 2021-02-04
journal: nan
DOI: 10.3390/axioms10010018
sha: 0a7f44538afde1826b3ba362149f9e6388a7fe44
doc_id: 289433
cord_uid: 57vaiail

The novel coronavirus disease (COVID-19) pneumonia has posed a great threat to the world recent months by causing many deaths and enormous economic damage worldwide. The first case of COVID-19 in Morocco was reported on 2 March 2020, and the number of reported cases has increased day by day. In this work, we extend the well-known SIR compartmental model to deterministic and stochastic time-delayed models in order to predict the epidemiological trend of COVID-19 in Morocco and to assess the potential role of multiple preventive measures and strategies imposed by Moroccan authorities. The main features of the work include the well-posedness of the models and conditions under which the COVID-19 may become extinct or persist in the population. Parameter values have been estimated from real data and numerical simulations are presented for forecasting the COVID-19 spreading as well as verification of theoretical results.

Coronaviruses are a large family of viruses that cause illnesses, ranging from the common cold to more serious illnesses such as Middle Eastern Respiratory Syndrome (MERS-CoV) and Severe Acute Respiratory Syndrome (SARS-CoV). The new coronavirus COVID-19 corresponds to a new strain that has not previously been identified in humans. On 11 March 2020, COVID-19 was reclassified as a pandemic by the World Health Organization (WHO). The disease has spread rapidly from country to country, causing enormous economic damage and many deaths around the world, prompting governments to issue a dramatic decree, ordering the lockdown of entire countries.

Since the confirmation of the first case of COVID-19 in Morocco on 2 March 2020 in the city of Casablanca, numerous preventive measures and strategies to control the spread of diseases have been imposed by the Moroccan authorities. In addition, Morocco declared a health emergency during the period from 20 March to 20 April 2020 and gradually extended it until 10 June 2020 in order to control the spread of the disease. In this paper, we report the assessment of the evolution of COVID-19 outbreak in Morocco. Besides shedding light on the dynamics of the pandemic, the practical intent of our analysis is to provide officials with the tendency of COVID-19 spreading, as well as gauge the effects of preventives measures using mathematical tools. Several other papers developed mathematical models for COVID-19 for particular regions in the globe and particular intervals of time, e.g., in [1] a Susceptible-Infectious-Quarantined-Recovered (SIQR) model to the analysis of data from the Brazilian Department of Health, obtained from 26 February 2020 to 25 March 2020 is proposed to better understand the early evolution of COVID-19 in Brazil; in [2] , a new COVID-19 epidemic model with media coverage and quarantine is constructed on the basis of the total confirmed new cases in the UK from 1 February 2020 to 23 March 2020; arXiv:2102.04260v1 [q-bio.PE] 4 Feb 2021

Based on the epidemiological feature of COVID-19 and the several strategies imposed by the government, with different degrees, to fight against this pandemic, we extend the classical SIR model to describe the transmission of COVID-19 in the Kingdom of Morocco. In particular, we divide the population into eight classes, denoted by S, I s , I a , F b , F g , F c , R and M, where S represents the susceptible individuals; I s the symptomatic infected individuals, which have not yet been treated; I a the asymptomatic infected individuals who are infected but do not transmit the disease; F b , F g and F c denote the patients diagnosed, supported by the Moroccan health system and under quarantine, and subdivided into three categories: benign, severe and critical forms, respectively. Finally, R and M are the recovered and fatality classes. This model satisfies the following assumptions:

(1) all coefficients involved in the model are positive constants;

(2) natural birth and death rate are not factors; (3) true asymptomatic patients will stay asymptomatic until recovery and do not spread the virus; (4) patients who are temporarily asymptomatic are included on symptomatic ones; (5) the second infection is not considered in the model; (6) the Moroccan health system is not overwhelmed.

According to the above assumptions and the actual strategies imposed by the Moroccan authorities, the spread of COVID-19 in the population is modeled by the following system of delayed differential equations (DDEs):

where t ∈ R + , N represents the total population size and u ∈ [0, 1] denotes the level of the preventive strategies on the susceptible population. The parameter β indicates the transmission rate and ∈ [0, 1] is the proportion for the symptomatic individuals. The parameter α denotes the proportion of the diagnosed symptomatic infected population that moves to the three forms: F b , F g and F c , by the rates γ b , γ g and γ c , respectively. The mean recovery period of these forms are denoted by 1/r b , 1/r g and 1/r c , respectively. The latter forms die also with the rates µ b , µ g and µ c , respectively. Asymptomatic infected population, which are not diagnosed, recover with rate η a and the symptomatic infected ones recover or die with rates η s and µ s , respectively. The time delays τ 1 , τ 2 , τ 3 and τ 4 denote the incubation period, the period of time needed before the charge by the health system, the time required before the death of individuals coming from the compartments I s , F b , F g , and F c , respectively. At each instant of time,

population is here assumed to be constant, that is, N(t) ≡ N during the period under study. This assumption is also reinforced by the fact that the Moroccan authorities have closed geographic borders.

The initial conditions of system (1) are

where τ = max{τ 1 , τ 2 , τ 3 , τ 4 }. Let C = C([−τ, 0], R 8 ) be the Banach space of continuous functions from the interval [−τ, 0] into R 8 equipped with the uniform topology. It follows from the theory of functional differential equations [16] that system (1) with initial conditions (ϕ 1 , ϕ 2 , ϕ 3 , ϕ 4 , ϕ 5 , ϕ 6 , ϕ 7 , ϕ 8 ) ∈ C has a unique solution. On the other hand, due to continuous fluctuation in the environment, the parameters of the system are actually not absolute constants and always fluctuate randomly around some average value. Hence, using delayed stochastic differential equations (DSDEs) to model the epidemic provide some additional degree of realism compared to their deterministic counterparts. The parameters β and α play an important role in controlling and preventing COVID-19 spreading and they are not completely known, but subject to some random environmental effects. We introduce randomness into system (1) by applying the technique of parameter perturbation, which has been used by many researchers (see, e.g., [17] [18] [19] ). In agreement, we replace the parameters β and α by β → β + σ 1Ḃ1 (t) and α → α + σ 2Ḃ2 (t), where B 1 (t) and B 2 (t) are independent standard Brownian motions defined on a complete probability space (Ω, F , P) with a filtration {F t } t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F 0 contains all P-null sets) and σ i represents the intensity of B i for i = 1, 2. Therefore, we obtain the following model governed by delayed stochastic differential equations:

where the coefficients are locally Lipshitz with respect to all the variables, for all t ∈ R + . Let us denote R 8

We have the following result. (3), there is a unique solution

to the COVID-19 stochastic model (4) that remains in R 8 + with a probability of one.

Proof. Since the coefficients of the stochastic differential equations with several delays (4) are locally Lipschitz continuous, it follows from [20] that for any square integrable initial value x(0) ∈ R 8 + , which is independent of the considered standard Brownian motion B, there exists a unique local solution x(t) on t ∈ [0, τ e ), where τ e is the explosion time. For showing that this solution is global, knowing that the linear growth condition is not verified, we need to prove that

For each integer k ≥ k 0 , we define the stopping time

Define the twice differentiable function W on R * 3 + → R + as follows:

T with the superscript "T " representing transposition, and L is the differential operator of function W defined by

Thus,

We now apply the elementary inequality 2xy ≤ x 2 + y 2 , valid for any x, y ∈ R, by firstly taking x = β (1 − u) and y = S(t) + I s (t) + I a (t) and, secondly,

In this way, we easily increase the right-hand side of inequality (5) to obtain that

between t 0 and t ∧ τ k and acting the expectation, which eliminates the martingale part, we get

EW(x(s)))ds and Gronwall's inequality implies that

For ω ∈ {τ k ≤ T}, x i (τ k ) equals k or 1 k for some i = 1, 2, 3. Hence,

It follows that

Letting k → ∞, we get P(τ e ≤ T) = 0. Since T is arbitrary, we obtain P(τ e = ∞) = 1. By defining the stopping

, k for some i = 4, . . . , 8 , and considering the twice differentiable functionW on

we deduce, with the same technique, that all the variables of the system are positive on [0, ∞).

The basic reproduction number, as a measure for disease spread in a population, plays an important role in the course and control of an ongoing outbreak [21] . This number is defined as the expected number of secondary cases produced, in a completely susceptible population, by a typical infective individual. Note that the calculation of the basic reproduction number R 0 does not depend on the variables of the system but depends on its parameters. In addition, the R 0 of our model does not depend on the time delays. For this reason, we use the next-generation matrix approach outlined in [22] to compute R 0 . Precisely, the basic reproduction number R 0 of system (1) is given by

where ρ is the spectral radius of the next-generation matrix FV −1 with

Noting that the classes that are directly involved in the spread of disease are only I s , I a , F b , F g and F c , we can reduce the local stability of system (1) to the local stability of

The other classes are uncoupled to the equations of system (1) and the total population size N is constant. Then, we can easily obtain the following analytical results:

Let E = (I s , I a , F b , F g , F c ) be an arbitrary equilibrium, and consider into system (7), the following change of unknowns:

By substituting U i (t), i = 1, 2, . . . , 5, into system (7) and linearizing around the free equilibrium, we get a new system that is equivalent to

where X(t) = (U 1 (t), U 2 (t), U 3 (t), U 4 (t), U 5 (t)) T and A, B, C are the Jacobian matrix of (7) given by

The characteristic equation of system (7) is given by P(λ) = (λ − a 1 (R 0 e −λτ 1 − 1))(λ + η a )(λ + (µ b + r b ))(λ + (µ g + r g ))(λ + (µ c + r c )),

where

Clearly, the characteristic Equation (10) has the roots λ 1 = −η a , λ 2 = −(µ b + r b ), λ 3 = −(µ g + r g ), λ 4 = −(µ c + r c ) and the root of the equation

We suppose Re(λ) ≥ 0. From (11), we get

On the other hand, we show that (11) has a real positive root when R 0 > 1. Indeed, we put

We have that Φ(0) = −a 1 (R 0 − 1) < 0, lim λ→+∞ Φ(λ) = +∞ and function Φ is continuous on (0, +∞). Consequently, Φ has a positive root and the following result holds. Theorem 2. The disease free equilibrium of system (1), that is, (N, 0, 0, 0, 0, 0, 0, 0) , is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1.

Knowing the value of the deterministic threshold R 0 characterizes the dynamical behavior of system (1) and guarantees persistence or extinction of the disease. Similarly, now we characterize the dynamical behavior of system (4) by a sufficient condition for extinction of the disease.

Then,

Namely, I s (t) tends to zero exponentially almost surely, that is, the disease dies out with a probability of one.

Let

To simplify, we set

Then, we get

Integrating both sides of the above inequality between 0 and t, one has

where

We have

Then,

From the large number theorem for martingales [23] , we deduce that

We also have

Then,

We conclude that if β 2 2σ 2 1 − α − (1 − α)(µ s + η s ) < 0, then lim I(t) t→∞ = 0. This completes the proof.

Estimating the model parameters poses a big challenge because the COVID-19 situation changes rapidly and from one country to another. The parameters are likely to vary over time as new policies are introduced on a day-to-day basis. For this reason, in order to simulate the COVID-19 models (1) and (4), we consider some parameter values from the literature, while the remaining ones are estimated or fitted.

As the transmission rate β is unknown, we carry out the least-square method [10] to estimate this parameter, based on the actual official reported confirmed cases from 2 March to 20 March, 2020 [24] . Through this method, we estimated β as 0.4517 (95%CI, 0.4484-0.455). Since the life expectancy for symptomatic individuals is 21 days on average and the crude mortality ratio is between 3% to 4% [25] , we estimated µ s = 0.01/21 per day and η s = 0.8/21 per day. Furthermore, since the hospitals are not yet saturated and the epidemic situation is under control, we assume that mortality comes mainly from critical forms with a percentage of 40% for an average period of 13.5 days [25] . Then, we choose µ c = 0.4/13.5 per day and r c = 0.6/13.5 per day. According to [26] , the proportion of asymptomatic individuals varies from 20.6% to 39.9% and of symptomatic individuals from 60.1% and 79.4% of the infected population. The progression rates γ b , γ g and γ c , from symptomatic infected individuals to the three forms, are assumed to be 80% of diagnosed cases for benign form, 15% of diagnosed cases for severe form, and 5% of diagnosed cases for critical form, respectively [25] . The incubation period is estimated to be 5.5 days [27, 28] while the time needed before hospitalization is to be 7.5 days [29] [30] [31] . Following a clinical observation related to the situation of COVID-19 in Morocco, an evolution of symptomatic individuals is estimated towards recovery or death after 21 days without any clinical intervention. In the case when clinical intervention is applied, we estimate the evolution of the critical forms towards recovery or death after 13.3 days. The rest of the parameter values are shown in Table 1 . 

In this section, we present the forecasts of COVID-19 in Morocco related to different strategies implemented by Moroccan authorities.

Taking into account the four levels of measures attached to containment, the effectiveness level of the applied Moroccan preventive measures is estimated to be In Figure 1 , we see that the plots and the clinical data are globally homogeneous. In addition, the last daily reported cases in Morocco [32], confirm the biological tendency of our model. Thus, our models are efficient to describe the spread of COVID-19 in Morocco. However, we note that some clinical data are far from the values of the models due to certain foci that appeared in some large areas or at the level of certain industrial areas. We conclude also that the stochastic behavior of COVID-19 presents certain particularities contrary to the deterministic one, namely the magnitude of its peak is higher and the convergence to eradication is faster. On the other hand, the conditions in Theorems 2 and 3 are verified. More precisely, the basic reproduction number R 0 = 0.5230 is less than one from 12 May 2020 and σ 2 1 = 1.0609 > 1.0598 =

, which means that the eradication of disease is ensured.

To prove the biological importance of delay parameters, we give the graphical results of Figure 2 , which describe the evolution of diagnosed positive cases with and without delays. We observe in Figure 2 , a high impact of delays on the number of diagnosed positive cases, thereby the plot of model (4) without delays (τ i = 0, i = 1, 2, 3, 4) is very different to that of the clinical data. Thus, we conclude that delays play an important role in the study of the dynamic behavior of COVID-19 worldwide, especially in Morocco, and allow us to better understand the reality.

In Figure 3 , we present the forecast of susceptible, severe forms of deaths and critical forms, from which we deduce that COVID-19 will not attack the total population. In addition, the number of hospitalization beds or artificial respiration apparatus required can be estimated by the number of different clinical forms. Moreover, we see that the number of deaths given by the model is less than those declared in other countries [33] , which shows that Morocco has avoided a dramatic epidemic situation by imposing the described strategies.

Finally, we present in Figure 4 , the cumulative diagnosed cases, severe forms, deaths and critical forms 240 days from the start of the pandemic in Morocco. We summarize some important numbers in Table 2 , which gives us some information about the future epidemic situation in Morocco. 

In this study, we proposed a new deterministic model with delay and its corresponding stochastic model to describe the dynamic behavior of COVID-19 in Morocco. These models provide us with the evolution and prediction of important categories of individuals to be monitored, namely, the positive diagnosed cases, which can help to examine the efficiency of the measures implemented in Morocco, and the different developed forms, which can quantify the capacity of the public health system as well as the number of new deaths. Firstly, we have shown that our models are mathematically and biologically well posed by proving global existence and uniqueness of positive solutions. Secondly, the extinction of the disease was established. By analyzing the characteristic equation, we proved that if R 0 < 1, then the disease free equilibrium of the deterministic model is locally asymptotically stable (Theorem 2). Based on the Lyapunov analysis method, a sufficient condition for the extinction was obtained in the stochastic case (Theorem 3). Thirdly, and since there is a substantial interest in estimating the parameters, we applied the least square method to determine the confidence interval of the transmission rate β as 0.4517 (95%CI, 0.4484-0.455). In addition, the rest of the parameters were either assumed, based on some daily observations, or taken from the available literature. Finally, some numerical simulations were performed to gather information in order to be able to fight against the propagation of the new coronavirus. In 12 May 2020, the basic reproduction number was less than one (R 0 = 0.5230), which means that the epidemic was tending toward eradication, which is conditional on strict compliance with the implemented measures. Currently, the consequences of the measures taken against COVID-19 in Morocco encourage their maintenance to control the spread of the epidemic and quickly move towards extinction.

Modeling the early evolution of the COVID-19 in Brazil: Results from a susceptible-infectious-quarantined-recovered (SIQR) model, Internat. J. Modern Phys. C 2020

Modelling the effects of media coverage and quarantine on the COVID-19 infections in the UK

Prediction of confinement effects on the number of Covid-19 outbreak in Algeria

Blow-up in a parabolic-elliptic Keller-Segel system with density-dependent sublinear sensitivity and logistic source

A singular elliptic problem related to the membrane equilibrium equations

Properties of solutions to porous medium problems with different sources and boundary conditions

Analysis and explicit solvability of degenerate tensorial problems

Estimation of the Transmission Risk of the 2019-nCoV and Its Implication for Public Health Interventions

Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study

Prediction of the Epidemic Peak of Coronavirus Disease in Japan, 2020

Analysis and forecast of COVID-19 spreading in China, Italy and France

Mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan

A simple stochastic SIR model for COVID 19 infection dynamics for Karnataka: Learning from Europe

A discrete stochastic model of the COVID-19 outbreak: Forecast and control

A stochastic epidemic model of COVID-19 disease

Introduction to Functional Differential Equations

Dynamics of a stochastic viral infection model with immune response

Qualitative analysis of a stochastic epidemic model with specific functional response and temporary immunity

A stochastic model of AIDS and condom use

Stochastic Differential Equations and Applications

On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations

Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission

A stochastic differential equations SIS epidemic model

Department of Epidemiology and Disease Control

Estimating the asymptomatic proportion of coronavirus disease 2019 (COVID-19) cases on board the Diamond Princess cruise ship

COVID-19 Incubation Period: An Update

Clinical features of patients infected with 2019 novel coronavirus in Wuhan

Clinical Characteristics of 138 Hospitalized Patients with 2019 Novel Coronavirus-Infected Pneumonia in Wuhan

Ministry of Health of Morocco. The Official Portal of Corona Virus in Morocco

COVID-19 Coronavirus Pandemic, View by Country

As future work, we intend to study the regional evolution of COVID-19 in Morocco. 

The data used in this study is available from the Government of Morocco, being given in Figure 1 .

We would like to express our gratitude to the editor and the anonymous reviewers, for their constructive comments and suggestions, which helped us to enrich the paper.

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.