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; while in [3] SEIR modelling to forecast the COVID-19 outbreak in Algeria is carried out by using available data from 1 March to 10 April 2020.

Mathematical modeling, particularly in terms of differential equations, is a strong tool that attracts the attention of many scientists to study various problems arising from mechanics, biology, physics, and so on. For instance, in [4], a system of differential equations with density-dependent sublinear sensitivity and logistic source is proposed and blow up properties of solutions are investigated; paper [5] presents a mathematical model with application in civil engineering related to the equilibrium analysis of a membrane with rigid and cable boundaries; [6] studies nonnegative and classical solutions to porous medium problems; and [7] a two-dimensional boundary value problem under proper assumptions on the data. Herein, we will focus on the dynamic of COVID-19. Tang et al. [8] used a Susceptible–Exposed–Infectious–Recovered (SEIR) compartmental model to estimate the basic reproduction number of COVID-19 transmission, based on data of confirmed cases for the disease in mainland China. Wu et al. [9] provided an estimate of the size of the epidemic in Wuhan on the basis of the number of cases exported from Wuhan to cities outside mainland China by using a SEIR model. In [10], Kuniya applied the SEIR compartmental model for the prediction of the epidemic peak for COVID-19 in Japan, using real-time data from 15 January to 29 February, 2020. Fanelli and Piazza [11] analyzed and forecasted the COVID-19 spreading in China, Italy and France, by using a simple Susceptible–Infected–Recovered–Deaths (SIRD) model. A more elaborate model, which includes the transmissibility of super-spreader individuals, is proposed in Ndaïrou et al. [12]. The model we propose here is new and has completely different compartments: in the paper [12], they model susceptible, exposed, symptomatic and infectious, super-spreaders, infectious but asymptomatic, hospitalized, recovered and the fatality class, with the main contribution being the inclusion of super-spreader individuals; in contrast, here we consider susceptible individuals, symptomatic infected individuals, which have not yet been treated, the asymptomatic infected individuals who are infected but do not transmit the disease, patients diagnosed and under quarantine and subdivided into three categories—benign, severe and critical forms—recovered and dead individuals. Moreover, our model has delays, while the previous model [12] has no delays; our model is stochastic, while the previous model [12] is deterministic. In fact, all mentioned models are deterministic and neglect the effect of stochastic noises derived from environmental fluctuations. To the best of our knowledge, research works that predict the COVID-19 outbreak taking into account a stochastic component, are a rarity [13,14,15]. The novelty of our work is twofold: the extension of the models cited above to a more accurate model with time delay, suggested biologically in the first place; secondly, to combine between the deterministic and the stochastic approaches in order to well-describe reality. To do this,Section 2 deals with the formulation and the well-posedness of the models.Section 3 is devoted to the qualitative analysis of the proposed models. Parameters estimation and forecast of COVID-19 spreading in Morocco is presented inSection 4 andSection 5, respectively. The paper ends with discussion and conclusions, inSection 6.

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 byS,$I_s$,$I_a$,$F_b$,$F_g$,$F_c$,R andM, whereS 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 andM are the recovered and fatality classes. This model satisfies the following assumptions:

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\in {\mathbb{R}}_{+}$,N represents the total population size and$u\in [0,1]$ denotes the level of the preventive strategies on the susceptible population. The parameter$\beta$ indicates the transmission rate and$ϵ\in [0,1]$ is the proportion for the symptomatic individuals. The parameter$\alpha$ 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${\gamma }_b$,${\gamma }_g$ and${\gamma }_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${\mu }_b$,${\mu }_g$ and${\mu }_c$, respectively. Asymptomatic infected population, which are not diagnosed, recover with rate${\eta }_a$ and the symptomatic infected ones recover or die with rates${\eta }_s$ and${\mu }_s$, respectively. The time delays${\tau }_1$,${\tau }_2$,${\tau }_3$ and${\tau }_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, gives the number of new death due to the disease (cf. [12]).

The initial conditions of system (1) are where$\tau =max\{{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4\}$. Let$\mathcal{C}=C([-\tau ,0],{\mathbb{R}}^8)$ be the Banach space of continuous functions from the interval$[-\tau ,0]$ into${\mathbb{R}}^8$ equipped with the uniform topology. It follows from the theory of functional differential equations [16] that system (1) with initial conditions 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$\beta$ and$\alpha$ 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$\beta$ and$\alpha$ by$\beta \to \beta +{\sigma }_1{\dot{B}}_1\left(t\right)$ and$\alpha \to \alpha +{\sigma }_2\dot{B_2}\left(t\right)$, where$B_1\left(t\right)$ and$B_2\left(t\right)$ are independent standard Brownian motions defined on a complete probability space$(\Omega ,\mathcal{F},\mathbb{P})$ with a filtration${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$ satisfying the usual conditions (i.e., it is increasing and right continuous while${\mathcal{F}}_0$ contains all P-null sets) and${\sigma }_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\in {\mathcal{R}}^{+}$.

Let us denote${\mathbb{R}}_{+}^8=\{(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8)\mid x_i>0,i=1,2,\dots ,8\}$. We have the following result.

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${\mathcal{R}}_0$ of system (1) is given by where$\rho$ is the spectral radius of the next-generation matrix$F{V}^{-1}$ with

$F=\left(\begin{array}{cc}\beta ϵ(1-u)& 0\\ 0& 0\end{array}\right)$ and$V=\left(\begin{array}{cc}(1-\alpha )({\eta }_s+{\mu }_s)+\alpha & 0\\ 0& {\eta }_a\end{array}\right)$.

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 sizeN is constant. Then, we can easily obtain the following analytical results:

Let$\overline{E}=(\overline{I_s},\overline{I_a},\overline{F_b},\overline{F_g},\overline{F_c})$ be an arbitrary equilibrium, and consider into system (7), the following change of unknowns:

By substituting$U_i\left(t\right)$,$i=1,2,\dots ,5$, into system (7) and linearizing around the free equilibrium, we get a new system that is equivalent to where$X\left(t\right)={(U_1\left(t\right),U_2\left(t\right),U_3\left(t\right),U_4\left(t\right),U_5\left(t\right))}^{\mathcal{T}}$ andA,B,C are the Jacobian matrix of (7) given by and

The characteristic equation of system (7) is given by where

Clearly, the characteristic Equation (10) has the roots${\lambda }_1=-{\eta }_a$,${\lambda }_2=-({\mu }_b+r_b)$,${\lambda }_3=-({\mu }_g+r_g)$,${\lambda }_4=-({\mu }_c+r_c)$ and the root of the equation

We suppose$Re\left(\lambda \right)\ge 0$. From (11), we get if${\mathcal{R}}_0<1$, which contradicts$Re\left(\lambda \right)\ge 0$. On the other hand, we show that (11) has a real positive root when${\mathcal{R}}_0>1$. Indeed, we put

We have that$\Phi \left(0\right)=-a_1({\mathcal{R}}_0-1)<0$,${lim}_{\lambda \to +\infty }\Phi \left(\lambda \right)=+\infty$ and function$\Phi$ is continuous on$(0,+\infty )$. Consequently,$\Phi$ has a positive root and the following result holds.

Knowing the value of the deterministic threshold${\mathcal{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.