Modeling and Forecasting of COVID-19 Spreading by Delayed Stochastic Differential Equations
Marouane Mahrouf, Adnane Boukhouima, Houssine Zine, El Mehdi Lotfi, Delfim F. M. Torres, Noura Yousfi
AArticle
Modeling and Forecasting of COVID-19 Spreading by DelayedStochastic Differential Equations
Marouane Mahrouf , Adnane Boukhouima , Houssine Zine , El Mehdi Lotfi , Delfim F. M. Torres * andNoura Yousfi Submitted
Axioms : 2 Dec 2020;Revised: 20 and 31 Jan 2021;This version: 4 Feb 2021.
Note:
This is a preprint of a paper whosefinal and definite form is published, openaccess, by
Axioms (ISSN: 2075-1680). Laboratory of Analysis, Modeling and Simulation (LAMS), Faculty of Sciences Ben M’sik, Hassan II University ofCasablanca, Sidi Othman, P.B. 7955 Casablanca, Morocco; [email protected] (M.M.);[email protected] (A.B.); lotfi[email protected] (E.M.L.); nourayousfi[email protected] (N.Y.) Center for Research and Development in Mathematics and Applications (CIDMA), Department of Mathematics,University of Aveiro, 3810-193 Aveiro, Portugal; [email protected] * Correspondence: delfi[email protected]
Abstract:
The novel coronavirus disease (COVID-19) pneumonia has posed a great threat to the worldrecent months by causing many deaths and enormous economic damage worldwide. The first case ofCOVID-19 in Morocco was reported on 2 March 2020, and the number of reported cases has increasedday by day. In this work, we extend the well-known SIR compartmental model to deterministic andstochastic time-delayed models in order to predict the epidemiological trend of COVID-19 in Moroccoand to assess the potential role of multiple preventive measures and strategies imposed by Moroccanauthorities. The main features of the work include the well-posedness of the models and conditionsunder which the COVID-19 may become extinct or persist in the population. Parameter values havebeen estimated from real data and numerical simulations are presented for forecasting the COVID-19spreading as well as verification of theoretical results.
Keywords:
COVID-19; coronaviruses; mathematical modeling; delayed stochastic differential equations(DSDEs)
1. Introduction
Coronaviruses are a large family of viruses that cause illnesses, ranging from the common cold to more serious illnessessuch as Middle Eastern Respiratory Syndrome (MERS-CoV) and Severe Acute Respiratory Syndrome (SARS-CoV). The newcoronavirus 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 fromcountry to country, causing enormous economic damage and many deaths around the world, prompting governments toissue 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, numerouspreventive measures and strategies to control the spread of diseases have been imposed by the Moroccan authorities. Inaddition, Morocco declared a health emergency during the period from 20 March to 20 April 2020 and gradually extended ituntil 10 June 2020 in order to control the spread of the disease. In this paper, we report the assessment of the evolution ofCOVID-19 outbreak in Morocco. Besides shedding light on the dynamics of the pandemic, the practical intent of our analysisis to provide officials with the tendency of COVID-19 spreading, as well as gauge the effects of preventives measures usingmathematical tools. Several other papers developed mathematical models for COVID-19 for particular regions in the globeand particular intervals of time, e.g., in [1] a Susceptible–Infectious–Quarantined–Recovered (SIQR) model to the analysisof data from the Brazilian Department of Health, obtained from 26 February 2020 to 25 March 2020 is proposed to betterunderstand the early evolution of COVID-19 in Brazil; in [2], a new COVID-19 epidemic model with media coverage andquarantine is constructed on the basis of the total confirmed new cases in the UK from 1 February 2020 to 23 March 2020; a r X i v : . [ q - b i o . P E ] F e b of 15 while in [3] SEIR modelling to forecast the COVID-19 outbreak in Algeria is carried out by using available data from 1 Marchto 10 April, 2020.Mathematical modeling, particularly in terms of differential equations, is a strong tool that attracts the attention of manyscientists to study various problems arising from mechanics, biology, physics, and so on. For instance, in [4], a system ofdifferential equations with density-dependent sublinear sensitivity and logistic source is proposed and blow up propertiesof solutions are investigated; paper [5] presents a mathematical model with application in civil engineering related to theequilibrium analysis of a membrane with rigid and cable boundaries; [6] studies nonnegative and classical solutions toporous 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 confirmedcases for the disease in mainland China. Wu et al. [9] provided an estimate of the size of the epidemic in Wuhan on the basisof the number of cases exported from Wuhan to cities outside mainland China by using a SEIR model. In [10], Kuniya appliedthe SEIR compartmental model for the prediction of the epidemic peak for COVID-19 in Japan, using real-time data from 15January to 29 February, 2020. Fanelli and Piazza [11] analyzed and forecasted the COVID-19 spreading in China, Italy andFrance, by using a simple Susceptible–Infected–Recovered–Deaths (SIRD) model. A more elaborate model, which includesthe transmissibility of super-spreader individuals, is proposed in Ndaïrou et al. [12]. The model we propose here is new andhas 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 beingthe inclusion of super-spreader individuals; in contrast, here we consider susceptible individuals, symptomatic infectedindividuals, which have not yet been treated, the asymptomatic infected individuals who are infected but do not transmitthe disease, patients diagnosed and under quarantine and subdivided into three categories—benign, severe and criticalforms—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 andneglect the effect of stochastic noises derived from environmental fluctuations. To the best of our knowledge, research worksthat predict the COVID-19 outbreak taking into account a stochastic component, are a rarity [13–15]. The novelty of our workis twofold: the extension of the models cited above to a more accurate model with time delay, suggested biologically in thefirst place; secondly, to combine between the deterministic and the stochastic approaches in order to well-describe reality. Todo this, Section 2 deals with the formulation and the well-posedness of the models. Section 3 is devoted to the qualitativeanalysis of the proposed models. Parameters estimation and forecast of COVID-19 spreading in Morocco is presented inSections 4 and 5, respectively. The paper ends with discussion and conclusions, in Section 6.
2. Models Formulation and Well-Posedness
Based on the epidemiological feature of COVID-19 and the several strategies imposed by the government, with differentdegrees, to fight against this pandemic, we extend the classical SIR model to describe the transmission of COVID-19 in theKingdom 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 patientsdiagnosed, 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 thefollowing 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. of 15 According to the above assumptions and the actual strategies imposed by the Moroccan authorities, the spread ofCOVID-19 in the population is modeled by the following system of delayed differential equations (DDEs): dS ( t ) dt = − β ( − u ) S ( t ) I s ( t ) N , dI s ( t ) dt = β(cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N − α I s ( t ) − ( − α )( µ s + η s ) I s ( t ) , dI a ( t ) dt = β ( − (cid:101) )( − u ) S ( t − τ ) I s ( t − τ ) N − η a I a ( t ) , dF b ( t ) dt = αγ b I s ( t − τ ) − (cid:0) µ b + r b (cid:1) F b ( t ) , dF g ( t ) dt = αγ g I s ( t − τ ) − (cid:0) µ g + r g (cid:1) F g ( t ) , dF c ( t ) dt = αγ c I s ( t − τ ) − (cid:0) µ c + r c (cid:1) F c ( t ) , dR ( t ) dt = η s ( − α ) I s ( t − τ ) + η a I a ( t − τ ) + r b F b ( t − τ ) + r g F g ( t − τ ) + r c F c ( t − τ ) , dM ( t ) dt = µ s ( − α ) I s ( t − τ ) + µ b F b ( t − τ ) + µ g F g ( t − τ ) + µ c F c ( t − τ ) , (1)where t ∈ R + , N represents the total population size and u ∈ [
0, 1 ] denotes the level of the preventive strategies on thesusceptible population. The parameter β indicates the transmission rate and (cid:101) ∈ [
0, 1 ] is the proportion for the symptomaticindividuals. The parameter α denotes the proportion of the diagnosed symptomatic infected population that moves to thethree 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 by1/ r b , 1/ r g and 1/ r c , respectively. The latter forms die also with the rates µ b , µ g and µ c , respectively. Asymptomatic infectedpopulation, 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 τ , τ , τ and τ denote the incubation period, the period of time needed before thecharge 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, D ( t ) = : µ s ( − α ) I s ( t − τ ) + µ b F b ( t − τ ) + µ g F g ( t − τ ) + µ c F c ( t − τ ) = dM ( t ) dt (2)gives the number of new death due to the disease (cf. [12]). Remark 1.
In system (1) , delays occur at the entrances, when the actions of infection take charge or the actions by the health systembegin, and not at exits. Let us see an example. A susceptible individual, after contact with an infected person at instant t, becomes himselfinfected at instant t + τ . Suddenly, the compartment of the infected is fed at the instant t by the susceptible infected at the instant t − τ .The same operation occurs at the level of the other interactions between the compartments of the model. Remark 2.
We assume that the compartment of symptomatic infected I s does not completely empty at any time t. For this reason, onehas µ s + η s < . Note also that the diagnosed symptomatic infected population is completely distributed into one of three possible forms:F b , F g and F c , respectively by the rates γ b , γ g and γ c . Then, γ b + γ g + γ c = . Remark 3.
Biologically, τ = days and τ = days are the time periods needed before dying, deriving from I s and the three formsF b , F g , F c , respectively. That is why we inserted these delays in the last equation of system (1) . Remark 4.
We consider only a short time period in comparison to the demographic time-frame. From a biological point of view, thismeans that we can assume that there is neither entry (recruitment rate) nor exit (natural mortality rate), and vital parameters can beneglected. Note also that in our model, the individuals that die due to the disease are included in the population. Therefore, the total of 15 population is here assumed to be constant, that is, N ( t ) ≡ N during the period under study. This assumption is also reinforced by thefact that the Moroccan authorities have closed geographic borders.
The initial conditions of system (1) are S ( θ ) = ϕ ( θ ) ≥ I s ( θ ) = ϕ ( θ ) ≥ I a ( θ ) = ϕ ( θ ) ≥ F b ( θ ) = ϕ ( θ ) ≥ F g ( θ ) = ϕ ( θ ) ≥ F c ( θ ) = ϕ ( θ ) ≥ R ( θ ) = ϕ ( θ ) ≥ M ( θ ) = ϕ ( θ ) ≥ θ ∈ [ − τ , 0 ] , (3)where τ = max { τ , τ , τ , τ } . Let C = C ([ − τ , 0 ] , R ) be the Banach space of continuous functions from the interval [ − τ , 0 ] into R equipped with the uniform topology. It follows from the theory of functional differential equations [16] that system(1) with initial conditions ( ϕ , ϕ , ϕ , ϕ , ϕ , ϕ , ϕ , ϕ ) ∈ C has a unique solution. On the other hand, due to continuous fluctuation in the environment, the parameters of the systemare actually not absolute constants and always fluctuate randomly around some average value. Hence, using delayedstochastic differential equations (DSDEs) to model the epidemic provide some additional degree of realism compared totheir deterministic counterparts. The parameters β and α play an important role in controlling and preventing COVID-19spreading and they are not completely known, but subject to some random environmental effects. We introduce randomnessinto system (1) by applying the technique of parameter perturbation, which has been used by many researchers (see, e.g.,[17–19]). In agreement, we replace the parameters β and α by β → β + σ ˙ B ( t ) and α → α + σ ˙ B ( t ) , where B ( t ) and B ( t ) are independent standard Brownian motions defined on a complete probability space ( Ω , F , P ) with a filtration {F t } t ≥ satisfying the usual conditions (i.e., it is increasing and right continuous while F contains all P-null sets) and σ i representsthe intensity of B i for i =
1, 2. Therefore, we obtain the following model governed by delayed stochastic differential equations: dS ( t ) = (cid:18) − β ( − u ) S ( t ) I s ( t ) N (cid:19) dt − σ ( − u ) S ( t ) I s ( t ) N dB ( t ) , dI s ( t ) = (cid:18) β(cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N − α I s ( t ) − ( − α )( µ s + η s ) I s ( t ) (cid:19) dt + σ (cid:18) (cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N (cid:19) dB ( t ) + σ ( µ s + η s − ) I s ( t ) dB ( t ) , dI a ( t ) = (cid:18) β ( − (cid:101) )( − u ) S ( t − τ ) I s ( t − τ ) N − η a I a ( t ) (cid:19) d ( t )+ σ ( − (cid:101) )( − u ) S ( t − τ ) I s ( t − τ ) N dB ( t ) , dF b ( t ) = (cid:18) αγ b I s ( t − τ ) − (cid:0) µ b + r b (cid:1) F b ( t ) (cid:19) dt + σ γ b I s ( t − τ ) dB ( t ) , dF g ( t ) = (cid:18) αγ g I s ( t − τ ) − (cid:0) µ g + r g (cid:1) F g ( t ) (cid:19) dt + σ γ g I s ( t − τ ) dB ( t ) , dF c ( t ) = (cid:18) αγ c I s ( t − τ ) − (cid:0) µ c + r c (cid:1) F c ( t ) (cid:19) dt + σ γ c I s ( t − τ ) dB ( t ) , dR ( t ) = (cid:18) η s ( − α ) I s ( t − τ ) + η a I a ( t − τ ) + r b F b ( t − τ ) + r g F g ( t − τ ) + r c F c ( t − τ ) (cid:19) dt − σ η s I s ( t − τ ) dB ( t ) , dM ( t ) = (cid:0) µ s ( − α ) I s ( t − τ ) + µ b F b ( t − τ ) + µ g F g ( t − τ ) + µ c F c ( t − τ ) (cid:1) dt − σ µ s I s ( t − τ ) dB ( t ) , (4)where the coefficients are locally Lipshitz with respect to all the variables, for all t ∈ R + .Let us denote R + = { ( x , x , x , x , x , x , x , x ) | x i > i =
1, 2, . . . , 8 } . We have the following result. of 15 Theorem 1.
For any initial value satisfying condition (3) , there is a unique solutionx ( t ) = ( S ( t ) , I s ( t ) , I a ( t ) , F b ( t ) , F g ( t ) , F c ( t ) , R ( t ) , M ( t )) to the COVID-19 stochastic model (4) that remains in R + with a probability of one. Proof.
Since the coefficients of the stochastic differential equations with several delays (4) are locally Lipschitz continu-ous, it follows from [20] that for any square integrable initial value x ( ) ∈ R + , which is independent of the consideredstandard Brownian motion B , there exists a unique local solution x ( t ) on t ∈ [ τ e ) , where τ e is the explosion time. Forshowing that this solution is global, knowing that the linear growth condition is not verified, we need to prove that τ e = ∞ . Let k > k < x ( ) < k . For each integer k ≥ k , we define the stopping time τ k = inf (cid:26) t ∈ [ τ e ) s.t. x i ( t ) / ∈ (cid:18) k , k (cid:19) for some i =
1, 2, 3 (cid:27) , where inf ∅ = ∞ . It is clear that τ k ≤ τ e . Let T > W on R ∗ + → R + as follows: W ( x ) = ( x + x + x ) + x + x + x .By Itô’s formula, for any 0 ≤ t ≤ τ k ∧ T and k ≥
1, we have dW ( x ( t )) = LW ( x ( t )) dt + ζ ( x ( t )) dB ( t ) ,where ζ is a continuous functional defined on [ + ∞ ) × C ([ − τ , 0 ] , R × ) by ζ ( x ( t )) = − σ ( − u ) S ( t ) I s ( t ) N σ (cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N σ ( µ s + η s − ) I s ( t ) σ ( − (cid:101) )( − u ) S ( t − τ ) I s ( t − τ ) N , B ( t ) = ( B ( t ) , B ( t )) T with the superscript “ T ” representing transposition, and L is the differential operator of function W defined by LW ( x ( t )) = (cid:18) ( S ( t ) + I s ( t ) + I a ( t )) − S ( t ) (cid:19)(cid:18) − β ( − u ) S ( t ) I s ( t ) N (cid:19) + (cid:18) + S ( t ) (cid:19)(cid:18) − σ ( − u ) S ( t ) I s ( t ) N (cid:19) + (cid:18) ( S ( t ) + I s ( t ) + I a ( t )) − I s ( t ) (cid:19)(cid:20) β(cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N − α I s ( t ) − ( − α )( µ s + η s ) I s ( t ) (cid:21) + (cid:18) + I s ( t ) (cid:19)(cid:34)(cid:18) σ (cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N (cid:19) + (cid:0) σ ( µ s + η s − ) I s ( t ) (cid:1) (cid:35) + (cid:18) ( S ( t ) + I s ( t ) + I a ( t )) − I a ( t ) (cid:19)(cid:18) β ( − (cid:101) )( − u ) S ( t − τ ) I s ( t − τ ) N − η a I a ( t ) (cid:19) + (cid:18) + I a ( t ) (cid:19)(cid:18) σ ( − (cid:101) )( − u ) S ( t − τ ) I s ( t − τ ) N (cid:19) . of 15 Thus, LW ( x ( t )) ≤ β ( − u ) S ( t ) I s ( t ) NS ( t ) + (cid:18) + S ( t ) (cid:19)(cid:18) σ ( − u ) S ( t ) I s ( t ) N (cid:19) + β(cid:101) ( − u ) (cid:0) S ( t ) + I s ( t ) + I a ( t ) (cid:1) S ( t − τ ) I s ( t − τ ) N + α + ( − α )( µ s + η s ) I s ( t )+ (cid:18) + I s ( t ) (cid:19)(cid:34)(cid:18) σ (cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N (cid:19) + (cid:0) σ ( µ s + η s − ) I s ( t ) (cid:1) (cid:35) + β ( − (cid:101) )( − u ) (cid:0) S ( t ) + I s ( t ) + I a ( t ) (cid:1) S ( t − τ ) I s ( t − τ ) N + η a I a ( t ) + (cid:18) + I a ( t ) (cid:19)(cid:18) σ ( − (cid:101) )( − u ) S ( t − τ ) I s ( t − τ ) N (cid:19) . (5)We now apply the elementary inequality 2 xy ≤ x + y , valid for any x , y ∈ R , by firstly taking x = β(cid:101) ( − u ) and y = S ( t ) + I s ( t ) + I a ( t ) and, secondly, x = β ( − (cid:101) )( − u ) and y = S ( t ) + I s ( t ) + I a ( t ) . In this way, we easily increase theright-hand side of inequality (5) to obtain that LW ( x ( t )) ≤ b + ψ (cid:0) S ( t ) + I s ( t ) + I a ( t ) (cid:1) + b S ( t ) + b I s ( t ) + b I a ( t ) ≤ D ( + W ( x ( t ))) ,where ψ , b , b , b , and b are positive constants and D = max ( ψ , b , b , b , b ) . By integrating both sides of equality dW ( x ( t )) = LW ( x ( t )) dt + ζ ( x ( t )) dB ( t ) between t and t ∧ τ k and acting the expectation, which eliminates the martingale part, we get E ( W ( x ( t ∧ τ k )) = E ( W ( x )) + E (cid:90) t ∧ τ k t LW ( x ( s ))) ds ≤ E ( W ( x )) + E (cid:90) t ∧ τ k t D ( + W ( x ( s ))) ds ≤ E ( W ( x )) + DT + (cid:90) t ∧ τ k t EW ( x ( s ))) ds and Gronwall’s inequality implies that E ( W ( x ( t ∧ τ k )) ≤ ( EW ( x ) + DT ) exp ( CT ) .For ω ∈ { τ k ≤ T } , x i ( τ k ) equals k or 1 k for some i =
1, 2, 3. Hence, W ( x i ( τ k )) ≥ (cid:18) k + k (cid:19) ∧ (cid:18) k + k (cid:19) .It follows that ( EW ( x ) + DT ) exp ( CT ) ≥ E (cid:16) χ { τ k ≤ T } ( ω ) W ( x τ k ) (cid:17) ≥ (cid:18) k + k (cid:19) ∧ (cid:18) k + k (cid:19) P ( τ k ≤ T ) . of 15 Letting k → ∞ , we get P ( τ e ≤ T ) =
0. Since T is arbitrary, we obtain P ( τ e = ∞ ) =
1. By defining the stoppingtime ˜ τ k = inf (cid:26) t ∈ [ τ e ) s.t. x i ( t ) / ∈ (cid:18) k , k (cid:19) for some i =
4, . . . , 8 (cid:27) , and considering the twice differentiable function ˜ W on R ∗ + → R + as ˜ W ( x ) = (cid:32) ∑ i = x i (cid:33) + ∑ i = x i ,we deduce, with the same technique, that all the variables of the system are positive on [ ∞ ) .
3. Qualitative Analysis of the Models
The basic reproduction number, as a measure for disease spread in a population, plays an important role in the courseand control of an ongoing outbreak [21]. This number is defined as the expected number of secondary cases produced, ina completely susceptible population, by a typical infective individual. Note that the calculation of the basic reproductionnumber R does not depend on the variables of the system but depends on its parameters. In addition, the R of our modeldoes not depend on the time delays. For this reason, we use the next-generation matrix approach outlined in [22] to compute R . Precisely, the basic reproduction number R of system (1) is given by R = ρ ( FV − ) = β(cid:101) ( − u )( − α )( η s + µ s ) + α , (6)where ρ is the spectral radius of the next-generation matrix FV − with F = (cid:18) β(cid:101) ( − u )
00 0 (cid:19) and V = (cid:18) ( − α )( η s + µ s ) + α η a (cid:19) .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 thelocal stability of system (1) to the local stability of dI s ( t ) dt = β(cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N − α I s ( t ) − ( − α )( µ s + η s ) I s ( t ) , dI a ( t ) dt = β ( − (cid:101) )( − u ) S ( t − τ ) I s ( t − τ ) N − η a I a ( t ) , dF b ( t ) dt = αγ b I s ( t − τ ) − (cid:0) µ b + r b (cid:1) F b ( t ) , dF g ( t ) dt = αγ g I s ( t − τ ) − (cid:0) µ g + r g (cid:1) F g ( t ) , dF c ( t ) dt = αγ c I s ( t − τ ) − (cid:0) µ c + r c (cid:1) F c ( t ) . (7)The other classes are uncoupled to the equations of system (1) and the total population size N is constant. Then, we can easilyobtain the following analytical results: S ( t ) = N − (cid:0) I s ( t ) + I a ( t ) + F b ( t ) + F g ( t ) + F c ( t ) + R ( t ) + M ( t ) (cid:1) , R ( t ) = (cid:90) t (cid:2) η s ( − α ) I s ( δ − τ ) + η a I a ( δ − τ ) + r b F b ( δ − τ ) + r g F g ( δ − τ ) + r c F c ( δ − τ ) (cid:3) d δ , M ( t ) = (cid:90) t (cid:2) µ s ( − α ) I s ( δ − τ ) + µ a I a ( δ − τ ) + µ b F b ( δ − τ ) + µ g F g ( δ − τ ) + µ c F c ( δ − τ ) (cid:3) d δ . (8)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: U ( t ) = I s ( t ) − I s , U ( t ) = I a ( t ) − I a , U ( t ) = F b ( t ) − F b , U ( t ) = F g ( t ) − F g and U ( t ) = F c ( t ) − F c . of 15 By substituting U i ( t ) , i =
1, 2, . . . , 5, into system (7) and linearizing around the free equilibrium, we get a new system that isequivalent to dX ( t ) dt = AX ( t ) + BX ( t − τ ) + CX ( t − τ ) , (9)where X ( t ) = ( U ( t ) , U ( t ) , U ( t ) , U ( t ) , U ( t )) T and A , B , C are the Jacobian matrix of (7) given by A = − α − ( − α )( µ s + η s ) − η a − ( µ b + r b ) − ( µ g + r g )
00 0 0 0 − ( µ c + r c ) , B = β(cid:101) ( − u ) β ( − (cid:101) )( − u ) ,and C = αγ b αγ g αγ c .The characteristic equation of system (7) is given by P ( λ ) = ( λ − a ( R e − λτ − ))( λ + η a )( λ + ( µ b + r b ))( λ + ( µ g + r g ))( λ + ( µ c + r c )) , (10)where a = α + ( − α )( µ s + η s ) .Clearly, the characteristic Equation (10) has the roots λ = − η a , λ = − ( µ b + r b ) , λ = − ( µ g + r g ) , λ = − ( µ c + r c ) and theroot of the equation λ − a ( R e − λτ − ) =
0. (11)We suppose Re ( λ ) ≥
0. From (11), we get Re ( λ ) = a ( R e − Re ( λ ) τ cos ( Im λ τ ) − ) < R <
1, which contradicts Re ( λ ) ≥
0. On the other hand, we show that (11) has a real positive root when R >
1. Indeed,we put Φ ( λ ) = λ − a ( R e − λτ − ) .We have that Φ ( ) = − a ( R − ) <
0, lim λ → + ∞ Φ ( λ ) = + ∞ and function Φ is continuous on ( + ∞ ) . 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 < andunstable if R > . Knowing the value of the deterministic threshold R characterizes the dynamical behavior of system (1) and guaranteespersistence or extinction of the disease. Similarly, now we characterize the dynamical behavior of system (4) by a sufficientcondition for extinction of the disease. of 15 Theorem 3.
Let x ( t ) = (cid:0) S ( t ) , I s ( t ) , I a ( t ) , F b ( t ) , F g ( t ) , F c ( t ) , R ( t ) , M ( t ) (cid:1) be the solution of the COVID-19 stochastic model (4) withinitial value x ( ) defined in (3). Assume that σ > β ( α + ( − α )( µ s + η s )) . Then, lim sup t → + ∞ ln I s ( t ) t <
0. (12)
Namely, I s ( t ) tends to zero exponentially almost surely, that is, the disease dies out with a probability of one. Proof.
Let d ln I s ( t ) = (cid:20) I s ( t ) (cid:18) β(cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N − α I s ( t ) − ( − α )( µ s + η s ) I s ( t ) (cid:19) − I s ( t ) (cid:32)(cid:18) σ β(cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N (cid:19) + (cid:0) σ ( µ s + η s − ) I s ( t ) (cid:1) (cid:33)(cid:35) dt + σ β(cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N I s ( t ) dB ( t ) + σ ( µ s + η s − ) dB ( t ) .To simplify, we set G ( t ) = (cid:101) ( − u ) S ( t − τ ) I s ( t − τ ) N , R ( t ) = σ β G ( t ) I s ( t ) , R = σ ( µ s + η s − ) , H = − α − ( − α )( µ s + η s ) .Then, we get d ln I s ( t ) = (cid:34) β G ( t ) I s ( t ) + H − (cid:32)(cid:18) σ G ( t ) I s ( t ) (cid:19) + R (cid:33)(cid:35) dt + R ( t ) dB ( t ) + R dB ( t )= (cid:34) − σ (cid:34)(cid:18) G ( t ) I s ( t ) (cid:19) − βσ G ( t ) I s ( t ) (cid:35) + H − R (cid:35) dt + R ( t ) dB ( t ) + R dB ( t )= − σ (cid:32) G ( t ) I s ( t ) − βσ (cid:33) − β σ + H − R dt + R ( t ) dB ( t ) + R dB ( t ) ≤ (cid:34) β σ + H (cid:35) dt + R ( t ) dB ( t ) + R dB ( t ) .Integrating both sides of the above inequality between 0 and t , one hasln I s ( t ) t ≤ ln I s ( ) t + β σ + H + M ( t ) t + M ( t ) t ,where M ( t ) = (cid:90) t R ( s ) dB ( s ) and M ( t ) = (cid:90) t R dB ( s ) .We have < M , M > t = (cid:90) t σ (cid:101) ( − u ) S ( s − τ ) I s ( s − τ ) N I s ( s ) ds ≤ (cid:90) t σ (cid:101) ( − u ) N N I s ( s ) ds ≤ (cid:90) t σ (cid:101) ( − u ) ds .Then, lim sup t → ∞ < M , M > t t ≤ σ (cid:101) ( − u ) < + ∞ .From the large number theorem for martingales [23], we deduce thatlim t → ∞ M ( t ) t = < M , M > t = (cid:90) t σ ( µ s + η s − ) ds = σ ( µ s + η s − ) t .Then, lim sup t → ∞ < M , M > t t ≤ σ ( µ s + η s − ) < + ∞ and lim t → ∞ M ( t ) t = t → + ∞ ln I s ( t ) t ≤ β σ − α − ( − α )( µ s + η s ) .We conclude that if β σ − α − ( − α )( µ s + η s ) <
0, then lim I ( t ) t → ∞ =
0. This completes the proof.
4. Assessment of Parameters
Estimating the model parameters poses a big challenge because the COVID-19 situation changes rapidly and from onecountry to another. The parameters are likely to vary over time as new policies are introduced on a day-to-day basis. For thisreason, in order to simulate the COVID-19 models (1) and (4), we consider some parameter values from the literature, whilethe 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 onthe actual official reported confirmed cases from 2 March to 20 March, 2020 [24]. Through this method, we estimated β as0.4517 (95%CI, 0.4484–0.455). Since the life expectancy for symptomatic individuals is 21 days on average and the crudemortality ratio is between 3% to 4% [25], we estimated µ s = η s = µ c = r c = γ b , γ g and γ c , fromsymptomatic infected individuals to the three forms, are assumed to be 80% of diagnosed cases for benign form, 15% ofdiagnosed cases for severe form, and 5% of diagnosed cases for critical form, respectively [25]. The incubation period isestimated to be 5.5 days [27,28] while the time needed before hospitalization is to be 7.5 days [29–31]. Following a clinicalobservation related to the situation of COVID-19 in Morocco, an evolution of symptomatic individuals is estimated towardsrecovery or death after 21 days without any clinical intervention. In the case when clinical intervention is applied, weestimate the evolution of the critical forms towards recovery or death after 13.3 days. The rest of the parameter values areshown in Table 1. Table 1.
Parameter values of models (1) and (4).
Parameter Value Source Parameter Value Source β u [0–1] Varied (cid:101) γ b γ g γ c α η a η s µ s µ b µ g µ c r b r g r c τ τ τ
21 Assumed τ σ σ
5. Numerical Simulation of Moroccan COVID-19 Evolution
In this section, we present the forecasts of COVID-19 in Morocco related to different strategies implemented by Moroccanauthorities.Taking into account the four levels of measures attached to containment, the effectiveness level of the applied Moroccanpreventive measures is estimated to be u = ( (
10 March, 20 March];0.4, on (
20 March, 6 April];0.8, after 6 April.In Figure 1, we see that the plots and the clinical data are globally homogeneous.
Dates
March 02 May 01 July 30 September 28 November 27 N u m b e r o f d i a gno se d c on f i r m e d cases Confirmed casesStochastic ModelDeterministic Model
Figure 1.
Comparison of the deterministic and the stochastic dynamical behavior with the daily reported cases of COVID-19 in Morocco.
In addition, the last daily reported cases in Morocco [32], confirm the biological tendency of our model. Thus, ourmodels are efficient to describe the spread of COVID-19 in Morocco. However, we note that some clinical data are far fromthe values of the models due to certain foci that appeared in some large areas or at the level of certain industrial areas. Weconclude 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 = σ = > = β ( α + ( − α )( µ s + η s )) , 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 theevolution of diagnosed positive cases with and without delays. Dates
March 02 May 01 July 30 September 28 November 27 N u m b e r o f d i a gno se d c on f i r m e d cases Confirmed casesModel with delaysModel without delays
Figure 2.
Effect of delays on the diagnosed confirmed cases.
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 = i =
1, 2, 3, 4) is very different to that of the clinical data. Thus, we conclude that delays play animportant role in the study of the dynamic behavior of COVID-19 worldwide, especially in Morocco, and allow us to betterunderstand the reality.In Figure 3, we present the forecast of susceptible, severe forms of deaths and critical forms, from which we deduce thatCOVID-19 will not attack the total population.
Days S u sce p t i b l es × Deterministic modelStochastic model
Days S eve r e f o r m s Days D ea t h s Days C r i t i ca l f o r m s Figure 3.
The evolution of susceptible, deaths, severe and critical forms from 2 March 2020.
In addition, the number of hospitalization beds or artificial respiration apparatus required can be estimated by thenumber 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 describedstrategies.Finally, we present in Figure 4, the cumulative diagnosed cases, severe forms, deaths and critical forms 240 days fromthe start of the pandemic in Morocco. We summarize some important numbers in Table 2, which gives us some informationabout the future epidemic situation in Morocco.
Dates
March 02 May 01 July 30 September 28 November 27 C u m u l a t i ve × Diagnosed confirmed casesSevere formsCritical formsDeaths
Figure 4.
Cumulative diagnosed cases, severe forms, critical forms and deaths 240 days from the start of the COVID-19 pandemic inMorocco.
Table 2.
Estimated peaks and cumulative of diagnosed cases, severe forms, critical forms and deaths.
Compartments Peak Cumulative
Diagnosed Around 190 18,890Severe forms Around 28 2233Critical forms Around 10 997Deaths Around 5 468
6. Conclusions
In this study, we proposed a new deterministic model with delay and its corresponding stochastic model to describethe dynamic behavior of COVID-19 in Morocco. These models provide us with the evolution and prediction of importantcategories of individuals to be monitored, namely, the positive diagnosed cases, which can help to examine the efficiency of themeasures implemented in Morocco, and the different developed forms, which can quantify the capacity of the public healthsystem as well as the number of new deaths. Firstly, we have shown that our models are mathematically and biologically wellposed 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 <
1, then the disease free equilibrium of the deterministicmodel is locally asymptotically stable (Theorem 2). Based on the Lyapunov analysis method, a sufficient condition for theextinction was obtained in the stochastic case (Theorem 3). Thirdly, and since there is a substantial interest in estimating theparameters, 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, ortaken from the available literature. Finally, some numerical simulations were performed to gather information in order to beable to fight against the propagation of the new coronavirus. In 12 May 2020, the basic reproduction number was less thanone ( R = As future work, we intend to study the regional evolution of COVID-19 in Morocco.
Author Contributions:
Conceptualization, M.M., A.B., H.Z., E.M.L., D.F.M.T. and N.Y.; Formal analysis, M.M., A.B., H.Z., E.M.L., D.F.M.T.and N.Y.; Investigation, M.M., A.B., H.Z., E.M.L., D.F.M.T. and N.Y.; Writing—original draft, M.M., A.B., H.Z., E.M.L., D.F.M.T. and N.Y.;Writing—review & editing, M.M., A.B., H.Z., E.M.L., D.F.M.T. and N.Y. All authors participated in the writing and reviewing of the paper.All authors have read and agreed to the published version of the manuscript.
Funding:
H.Z. and D.F.M.T. were supported by FCT within project UIDB/04106/2020 (CIDMA).
Institutional Review Board Statement:
Not applicable.
Informed Consent Statement:
Not applicable.
Data Availability Statement:
The data used in this study is available from the Government of Morocco, being given in Figure 1.
Acknowledgments:
We would like to express our gratitude to the editor and the anonymous reviewers, for their constructive commentsand suggestions, which helped us to enrich the paper.
Conflicts of Interest:
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.
References
1. Crokidakis, N. Modeling the early evolution of the COVID-19 in Brazil: Results from a susceptible-infectious-quarantined-recovered(SIQR) model, Internat.
J. Modern Phys. C , , 2050135.2. Feng, L.-X.; Jing, S.-L.; Hu, S.-K.; Wang, D.-F.; Huo, H.-F. Modelling the effects of media coverage and quarantine on the COVID-19infections in the UK. Math. Biosci. Eng. , , 3618–3636.3. Moussaoui, A.; Auger, P. Prediction of confinement effects on the number of Covid-19 outbreak in Algeria. Math. Model. Nat. Phenom. , , 14.4. Tanaka, Y.; Yokota, T. Blow-up in a parabolic-elliptic Keller-Segel system with density-dependent sublinear sensitivity and logisticsource. Math. Methods Appl. Sci. , , 7372–7396.5. Viglialoro, G.; Murcia, J. A singular elliptic problem related to the membrane equilibrium equations. Int. J. Comput. Math. , ,2185–2196.6. Li, T.; Pintus, N.; Viglialoro, G. Properties of solutions to porous medium problems with different sources and boundary conditions. Z. Angew. Math. Phys. , , 18.7. Li, T.; Viglialoro, G. Analysis and explicit solvability of degenerate tensorial problems. Bound. Value Probl. , , 13.8. Tang, B.; Wang, X.; Li, Q.; Bragazzi, N.L.; Tang, S.; Xiao, Y.; Wu, J. Estimation of the Transmission Risk of the 2019-nCoV and ItsImplication for Public Health Interventions. J. Clin. Med. , , 462.9. Wu, J.T.; Leung, K.; Leung, G.M. Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoVoutbreak originating in Wuhan, China: a modelling study. Lancet , , 689–697.10. Kuniya, T. Prediction of the Epidemic Peak of Coronavirus Disease in Japan, 2020. J. Clin. Med. , , 789.11. Fanelli, D.; Piazza, F. Analysis and forecast of COVID-19 spreading in China, Italy and France. Chaos Solitons Fractals , , 109761.12. Ndaïrou, F.; Area, I.; Nieto, J.J.; Torres, D.F.M. Mathematical modeling of COVID-19 transmission dynamics with a case study ofWuhan. Chaos Solitons Fractals , , 109846.13. Simha, A.; Prasad, R.V.; Narayana, S. A simple stochastic SIR model for COVID 19 infection dynamics for Karnataka: Learning fromEurope. arXiv , arXiv:2003.11920.14. He, S.; Tang, S.; Rong, L. A discrete stochastic model of the COVID-19 outbreak: Forecast and control. Math. Biosci. Eng. , ,2792–2804.15. Bardina, X.; Ferrante, M.; Rovira, C. A stochastic epidemic model of COVID-19 disease. arXiv , arXiv:2005.02859.16. Hale, J.; Lunel, S.M.V. Introduction to Functional Differential Equations ; Springer: New York, NY, USA, 1993.17. Mahrouf, M.; Hattaf, K.; Yousfi, N. Dynamics of a stochastic viral infection model with immune response.
Math. Model. Nat. Phenom. , , 15–32.18. Hattaf, K.; Mahrouf, M.; Adnani, J.; Yousfi, N. Qualitative analysis of a stochastic epidemic model with specific functional responseand temporary immunity. Phys. A Stat. Mech. Appl. , , 591–600.19. Dalal, N.; Greenhalgh, D.; Mao, X. A stochastic model of AIDS and condom use. J. Math. Anal. Appl. , , 36–53.20. Mao, X. Stochastic Differential Equations and Applications ; Elsevier: Amsterdam, The Netherlands, 2007.21. Diekmann, O.; Heesterbeek, J.A.P.; Metz, J.A.J. On the definition and the computation of the basic reproduction ratio R in models forinfectious diseases in heterogeneous populations. J. Math. Biol. , , 365–382.
22. Driessche, P.V.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of diseasetransmission.
Math. Biosci. , , 29–48.23. Grai, A.; Greenhalgh, D.; Hu, L.; Mao, X.; Pan, J. A stochastic differential equations SIS epidemic model. SIAM J. Appl. Math. , Coronavirus Disease 2019 (COVID-19) ; Situation Report 46, 6 March 2020; WHO: Geneva, Switzerland, 2020.26. Mizumoto, K.; Kagaya, K.; Zarebski, A.; Chowell, G. Estimating the asymptomatic proportion of coronavirus disease 2019 (COVID-19)cases on board the Diamond Princess cruise ship, Yokohama, Japan, 2020.
Euro Surveill. , , 2000180.27. WHO. Coronavirus Disease 2019 (COVID-19)
Lancet , , 497–506.30. Wang, D.; Hu, B.; Hu, C.; Zhu, F.; Liu, X.; Zhang, J.; Wang, B.; Xiang, H.; Cheng, Z.; Xiong, Y.; Zhao, Y. Clinical Characteristics of 138Hospitalized Patients with 2019 Novel Coronavirus-Infected Pneumonia in Wuhan, China. J. Am. Med. Assoc. ,323