Probabilistic Modelling of COVID-19 Dynamic in the Context of Madagascar

Angelo Raherinirina, Tsilefa Stefana Fandresena, Aimé Richard Hajalalaina, Haja Rabetafika, Rivo Andry Rakotoarivelo, Fontaine Rafamatanantsoa

University of Fianarantsoa, Madagascar.

DOI: 10.4236/ojmsi.2021.93014                

Abstract

We propose a probabilistic approach to modelling the propagation of the coronavirus disease 2019 (COVID-19) in Madagascar, with all its specificities. With the strategy of the Malagasy state, which consists of isolating all suspected cases and hospitalized confirmed case, we get an epidemic model with seven compartments: susceptible (S), Exposed (E), Infected (I), Asymptomatic (A), Hospitalized (H), Cured (C) and Death (D). In addition to the classical deterministic models used in epidemiology, the stochastic model offers a natural representation of the evolution of the COVID-19 epidemic. We inferred the models with the official data provided by the COVID-19 Command Center (CCO) of Madagascar, between March and August 2020. The basic reproduction number R₀ and the other parameters were estimated with a Bayesian approach. We developed an algorithm that allows having a temporal estimate of this number with confidence intervals. The estimated values are slightly lower than the international references. Generally, we were able to obtain a simple but effective model to describe the spread of the disease.

Keywords

Modified SEIR Model, COVID-19 Madagascar, Basic Reproduction Number, Markov Chain Continuous Time

1. Introduction

Since its emergence in Wuhan, coronavirus 2 (SARS-CoV-2) has shaken up scientific communities around the world [1]. Currently, there are more than 100,000 articles published according to the database “COVID-19 primer” (https://covid19primer.com) most of whom are interested in modeling and forecasting [2] [3] [4]. The case of less developed countries in sub-Saharan Africa such as Madagascar is still less studied because of insufficient data. For Madagascar, since the first case in March 2019, the COVID-19 has affected more than 14,696 people until 30 August 2020 according to official statistics. Of the contaminated, 13,709 are cured with 187 deaths officially declared by the ministry. The cure rate is, therefore, close to 93.28% and the mortality rate is about 1.30%.

This statistic is widely questionable given the insufficient number of tests on the big island. The absence of a massive test and its management principle make the choice of the procedure for responding to the disease quite particular in Madagascar. Socio-economic factors do not allow containment measures to be strictly applied. This specificity extends to the modeling approach of the evolution of the pandemic to represent reality [5] [6]. Studies have been conducted on the identification of the various symptoms and clinical signs of COVID-19 to anticipate the detection of confirmed cases.

From the compartmental models of epidemics introduced by Daniel Bernoulli [7]; several extensions have been proposed to explain the evolution of more diseases [8] [9] [10]. These models have been widely used by scientists around the world to model the evolution of the COVID-19 epidemic. Several compartments can be used depending on the modeling context. Faïçal et al. [8] consider an 8 compartments model and Kyrychko et al. arrived up to 10 compartments [11]. The number of parameters as well as the complexity of the model increase with this number of compartments. Our challenge in this article is to have a relatively simple, understandable but realistic model in the context of Madagascar.

As in all epidemic models, the basic reproduction number gives us information on the spread of the disease. In order to have a better estimate of this parameter, we combine several approaches starting from a deterministic model in order to study its stochastic extension. For the deterministic case, the main results on the ordinary differential equation make it possible to identify important properties on the propagation of the disease [12] [13]. The basic reproduction number is considered to be constant.

Considering the fact that people do not have the same sensitivities to contact with the virus and different factors can influence the rate of spread of the disease; the use of probabilistic models seems realistic to us. From section 4, we focus on the stochastic extension of the deterministic model, with the introduction of random noise. The techniques used are detailed in the article by Alen et al. [14].

In [15], with its six-compartment model, authors introduced a similar approach to model the evolution of the new coronavirus in the island of Mayotte (French). We have adapted this model for the case of Madagascar with our seven compartments epidemic model. We propose new methods for simulation and estimation of parameters. The models will be inferred with official data on COVID-19 in Madagascar, during the critical moment of the epidemic, i.e. between March and August 2020.

2. Preliminaries

2.1. Local Context of the COVID-19 Epidemic in Madagascar

Madagascar remains vulnerable to coronavirus attack. Despite the small number of cases detected, the impact of the new coronavirus on the Malagasy health system is considerable. The country is still at risk of an unmanageable and dangerous health crisis.

As soon as the virus appeared in Madagascar, the CCO took the initiative to introduce restrictions in order to limit the spread of the virus. Studies have shown that this initiative is early and has not provided more information on how the virus spreads.

According to official data, during the first wave, between March and August 2020, the outbreaks are in two regions: Analamanga region (Antananarivo, capital) and Antsinana region (Tamatave) (Figure 1).

Much of the data used in this article comes from these two strategic regions of Madagascar: Antananarivo the capital and Toamasina the economic capital. Figure 2 and Figure 3 represent the evolution of new cases observed in these two regions.

These graphs show that the epidemic has not had the same dynamic in these two major regions of Madagascar. It first hit the east side of the big island before entering the highlands.

For the entire territory of Madagascar, the SRAS-COV2 coronavirus epidemic has experienced several phases. The numbers of new cases and deaths were considered constant during the first trimester. These numbers followed an upward phase from June 2020. The impact of the epidemic was major in July 2020. The peak of this first wave is assumed to have reached around the beginning of August 2020, see Figure 4.

Figure 5 describes the cumulative data observed for all 22 regions of Madagascar.

2.2. Modelling approach

The effectiveness of SIR (Susceptible-infected-removed) models has already been shown in epidemiological modelling problems. Several modifications have been proposed based on hypotheses on how the disease is transmitted [7] [8].

To limit the spread of the COVID-19 virus, Madagascar’s strategy is to isolate all suspected cases and hospitalize confirmed case. Given the inadequacy or lack of mass screening (PCR test), the effects of this strategy remain questionable. Many people with the virus are asymptomatic and show no signs of infection.

In line with this strategy, we propose in this article a model with seven compartments:

• S (Susceptibles): individuals who are not infected but at some risk of becoming infected;

• E (Exposed): individuals who have been in contact with the virus but not yet infectious until after a latency period;

Figure 1. Official data observed for the 22 regions of Madagascar during the first wave of the covid epidemic 19.

Figure 2. Evolution of new cases of COVID-19 in the Antsinanana (Toamasina) region during the first wave.

Figure 3. Evolution of new cases of COVID-19 in the Analamanga (Antananarivo) region during the first wave.

Figure 4. Daily evolution of the numbers of new cases, healed and dead.

Figure 5. Cumulative numbers of confirmed, cured and dead cases for all 22 regions of Madagascar during the first wave.

• I (Infectious): individuals who can infect and are not yet hospitalized or isolated,

• H (Hospitalized): confirmed and hospitalized (isolated) cases;

• A (Asymptomatic): infected individuals but do not have symptoms (not isolated);

• C (Cured): healed confirmed cases;

• D (Death): deaths due to COVID-19.

Diagram in Figure 6 describes the evolution of our system according to the flows between the different compartments. The particularity of this model among others are first of all the separation of healthy (susceptible) individuals and individuals who are not yet sick but already in contact with the virus (Exposed). We have also put in another compartment infected individuals asymptomatic but who can transmit the virus.

Using the assumptions of Manou-Abi and Balicchi in [15] all the models’ parameters is summarized in Table 1.

The basic reproductions number R₀ corresponding to the model will be estimated using the officiad data recenced by the CCO-COVID-19. It represents the average number of secondary cases produced by a fundamental case over its infectious period. It is the initial growth rate of the epidemic.

3. Classic Deterministic Model

3.1. Deterministic SEIR Modified Model

As a reference, we will consider the deterministic case adapted to the context of.

Figure 6. COVID-19 modified SEIR model with seven compartements.

Table 1. Parameters of the COVID-19 SEIR modified model.

Madagascar. This type of model is widely used in mathematical modeling of the spread of the COVID-19 [5] [16].

At a moment t, let $S\left(t\right),E\left(t\right),I\left(t\right),A\left(t\right),H\left(t\right),C\left(t\right)$ et $D\left(t\right)$ the number of individuals in the compartments “Susceptible”, “Exposed”, “Infections”, “Asymptomatic”, “Hospitalized”, “Cured” et “Death”.

Assume that the total number of population is constant, i.e.,

$N=N\left(t\right)$, (1)

with $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+A\left(t\right)+H\left(t\right)+C\left(t\right)+D\left(t\right)$, for any momentt.

In a time interval $\left[t,t+\text{d}t\right]$, the evolution of the system is described by Equation (2):

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=-\alpha \frac{S\left(t\right)I\left(t\right)}{N},\\ \frac{\text{d}E\left(t\right)}{\text{d}t}=\alpha \frac{S\left(t\right)I\left(t\right)}{N}-\beta E\left(t\right),\\ \frac{\text{d}I\left(t\right)}{\text{d}t}=\beta E\left(t\right)-p\gamma I\left(t\right)-\left(1-p\right)\lambda I\left(t\right),\\ \frac{\text{d}A\left(t\right)}{\text{d}t}=p\gamma I\left(t\right),\\ \frac{\text{d}H\left(t\right)}{\text{d}t}=\left(1-p\right)\lambda I\left(t\right)-\left(1-q\right)\delta H\left(t\right)-q\theta H\left(t\right),\\ \frac{\text{d}C\left(t\right)}{\text{d}t}=q\theta H\left(t\right),\\ \frac{\text{d}D\left(t\right)}{\text{d}t}=\left(1-q\right)\delta H\left(t\right).\end{array}$ (2)

with:

- $\left(1-p\right)\lambda I\left(t\right)\text{d}t$ the number of confirmed cases (isolated, hospitalized),

- $p\gamma I\left(t\right)$ the number of asymptomatic cases,

- $\alpha \left(S\left(t\right)I\left(t\right)\right)/N\text{d}t$ the number of new infections in the interval $\left[t,t+\text{d}t\right]$.

By grouping asymptomatic cases and healed people into R(t) (“removed”):

$R\left(t\right)=C\left(t\right)+A\left(t\right),$ (3)

the simulation of the model with an extract from the population, in an environment where the virus circulates, gives the graphs of Figure 7.

3.2. Estimation of the Parameters and Simulation

To estimate the infection rate β and recovery rate γ in the SIR sub-process, we use the method introduced by Kermarck and McKendick in [17] using the observed data.

(The parameters β and γ of the SIR sub-process).The number of basic reproductions R₀ is obtained by

$R_0=\beta /\gamma$, (4)

after training the ODE system of Equation (3) with the data of the COVID-19 in Madagascar.

Note:If R₀ < 1, the infected individual contaminates less than one other individual on average. This is why the epidemic is disappearing from the population.

If R₀ > 1, then the virus continues to spread in the population.

For Madagascar, by running the model with the observed data and its population of about 26,026,000, we obtain:

$\beta =11.5734$ et $\gamma =11.5274$ , and $R_0=11.5734/9.5245=1.2152$

Figure 7. Simulation of the model with a population extract in an environment where the virus circulates.

The estimated number of basic reproductions with official data is quite small. This may be due to the specificity of the Malagasy on the vulnerability to the COVID-19. Figure 8 shows the simulation of the model with the epidemic parameters of Madagascar. Without intervention, the dynamics stabilize over the long term, see Figure 9.

4. Stochastic Model

In this section, we study the stochastic pattern of the propagation dynamics of COVID-19 in Madagascar with a continuous time Makov chain (CTMC) [18]. In complex and imprecise situations, the use of stochastic models is always a better alternative. These models are very useful in epidemiology [19] [20] [21].

4.1. Stochastic Evolution of the Spread of COVID-19

The spread and speeds of transmission of the new coronavirus vary widely from country to country. It depends on many factors such as demographic, socioeconomic, geographical, biologic and especially epidemic [22]. At the same time, individuals strive to find good ways to deal with the disease and adapt their behaviors to the different requirements of the riposte strategy.

We assume that factors external to the virus affect its rate of spread and turn it into a random variable. Thus, the process formed by transmission speeds $\left(\alpha \left(t\right)\right),t\ge 0$ can be considered a stochastic process in continuous time with:

$\alpha \left(t\right)\text{d}t=\alpha \text{d}t+\sigma \text{d}B\left(t\right)$ (5)

where B is a standard Brownian movement [23].

$\alpha \left(t\right)\text{d}t$ represents the number of infection transmitted by an infectious individual during a period dt. It is a normally distributed random variable [15]:

$\alpha \left(t\right)\text{d}t~\mathcal{N}\left(\alpha \text{d}t,{\sigma }^2\text{d}t\right)$.

Figure 8. Simulation of the modified SEIR model Madagascar data and the estimated R₀ = 1.2 parameter.

Figure 9. Long-term SEIAHD process dynamics.

4.2. Continuous-Time Markov Chain SEIR Modified Model

The works of Norris [18] and Durret [23] are good references to the theories on continuous time Markov chain. We use results that are demonstrated in these works.

For the SEIAHCD process of Figure 6, for each individual, consider the following events:

i) Exposition: a susceptible individual became exposed through contact with an infected individual, under a probability of p₁

ii) Infection: an exposed individual became infected under a probability of p₂

iii) Recovery: an infected individual has passed the incubation period and become an asymptomatic case, with a probability of p₃

iv) Hospitalization: an infected individual is hospitalized with a probability of p₄

v) Healing: an individual hospitalized is cured of COVID-19

vi) Death: an individual hospitalized died of COVID-19

Events i), ii), iii), iv), v) et vi) correspond to jumps in processes $S\left(t\right),E\left(t\right),I\left(t\right),A\left(t\right),H\left(t\right),C\left(t\right)$ et $D\left(t\right)$. Consider the following subset of events:

- Events after exposition:

$G_1=\left\{\text{exposition},\text{infection},\text{removed}\right\}$, where

$\text{removed}=\left\{\text{revovery},\text{hospitalization}\right\}$.

- Events after infections:

$G_2=\left\{\text{recovery},\text{hospitalization},\text{retired}\right\}$, where

$\text{retired}=\left\{\text{death},\text{healing}\right\}$, with

- Events after hospitalization:

$G_3=\left\{\text{death},\text{healing}\right\}$.

In our modeling approach, we assume that events in each subset are independent uniformly distributed. Probabilities of jumps corresponding to the different events can be written:

- For G₁:

$\left(\text{exposition}\right)~\frac{p_1}{p_1+p_2+p_3+p_4}$, $\left(\text{infection}\right)~\frac{p_2}{p_1+p_2+p_3+p_4}$,

$\left(\text{removed}\right)~\frac{p_3+p_4}{p_1+p_2+p_3+p_4}$.

- For G₂

$\left(\text{recovery}\right)~\frac{p_3}{p_3+p_4+p_5+p_6}$, $\left(\text{hospitalization}\right)~\frac{p_4}{p_3+p_4+p_5+p_6}$,

$\left(\text{retired}\right)~\frac{p_5+p_6}{p_3+p_4+p_5+p_6}$.

- For G₃:

$\left(\text{healing}\right)~\frac{p_5}{p_5+p_6}$, $\left(\text{death}\right)~\frac{p_4}{p_3+p_4+p_5+p_6}$.

with the below assumption and with an additional hypotheses that the transition time are exponentially distributed, we can consider that processes $S\left(t\right),E\left(t\right),I\left(t\right),A\left(t\right),H\left(t\right),C\left(t\right)$ et $D\left(t\right)$ can be modeled by continuous time Markov chain. Algorithm 1 can simulate this model with the parameters estimated using official data. The result is in Figure 10.

5. Extending the Model to a Stochastic Differential Equation

5.1. Stochastic Equivalent of the Deterministic SEIR Model

The introduction of random effects into Equations (2) leads to a stochastic differential equation system. The deduced model is the stochastic equivalent of the modified SEIR model [24].

In a general way, a classic differential equation in explicit form

$\frac{\text{d}{x}_t}{\text{d}t}=b_t(\; x\; t\; )$

transforms into a stochastic differential equation by adding random noise.

This noise can be taken into account by a process of Winner W_t [25]:

$\frac{\text{d}{x}_t}{\text{d}t}=\mu \left(x_t\right)+\sigma \left(x_t\right)\text{d}{W}_t$, (6)

where $\mu \left(X,\theta \right)$ the average et ${\sigma }^2\left(X,\theta \right)$ the covariance matrix.

Algorithm 1. Simulation of the continuous time Markov chain SEIAHCD model.

Figure 10. Simulation of a COVID-19 propagation model with Markov chains in continuous time.

1) First, consider the sub-model SEI (Susceptible, Exposed, Infected). Let $\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right)$, associated random variables. For a small interval of time ∆t; transition probabilities are written:

$P\left(S+k,E+j,I+m\right)=\{\begin{array}{l}{p}_1=\frac{\alpha SI}{N}\Delta t,\left(k,j,m\right)=\left(-1,1,0\right)\\ p_2=\beta E\Delta t,\left(k,j,m\right)=\left(0,-1,1\right)\\ p_3=\left(p\gamma I+\left(1-p\right)\lambda I\right)\Delta t,\left(k,j,m\right)=\left(0,0,-1\right)\\ 1-\left({\displaystyle \underset{i=1}{\overset{3}{\sum }}{p}_i}\right)\Delta t,\left(k,j,m\right)=\left(0,0,0\right)\\ 0,\text{otherwise}\end{array}$ (7)

$\begin{array}{c}E\left(\Delta x\right)={\displaystyle \underset{i=1}{\overset{3}{\sum }}{p}_i{\left(\Delta x\right)}_i}=p_1{\left(\Delta x\right)}_1+p_2{\left(\Delta x\right)}_2+p_3{\left(\Delta x\right)}_3\\ =p_1\left(\begin{array}{c}-1\\ 1\\ 0\end{array}\right)+p_2\left(\begin{array}{c}0\\ -1\\ 1\end{array}\right)+p_3\left(\begin{array}{c}0\\ 0\\ -1\end{array}\right)\end{array}$

$E\left(\Delta x\right)=\left(\begin{array}{c}-\frac{\alpha SI}{N}\\ \frac{\alpha SI}{N}-\beta E\\ \beta E-p\gamma I-\left(1-p\right)\lambda I\end{array}\right)\Delta t=\mu \left(x_1,x_2,x_3\right)\Delta t$

$\begin{array}{c}E\left(\left(\Delta x\right)\left(\Delta x\right)\right)={\displaystyle \underset{i=1}{\overset{3}{\sum }}{p}_i{\left(\Delta x\right)}_i{\left(\Delta x\right)}_i^{\text{T}}}\\ =p_1\left(\begin{array}{c}-1\\ 1\\ 0\end{array}\right)\left(\begin{array}{ccc}-1& 1& 0\end{array}\right)+p_2\left(\begin{array}{c}0\\ -1\\ 1\end{array}\right)\left(\begin{array}{ccc}0& -1& 1\end{array}\right)+p_3\left(\begin{array}{c}0\\ 0\\ -1\end{array}\right)\left(\begin{array}{ccc}0& 0& -1\end{array}\right)\\ =\left(\begin{array}{ccc}{p}_1& -p_1& 0\\ -p_1& p_1+p_2& -p_2\\ 0& -p_2& p_2+p_3\end{array}\right)\\ =\left(\begin{array}{ccc}\frac{\alpha SI}{N}& -\frac{\alpha SI}{N}& 0\\ -\frac{\alpha SI}{N}& \frac{\alpha SI}{N}+\beta E& -\beta E\\ 0& -\beta E& \beta E+p\gamma I+\left(1-p\right)\lambda I\end{array}\right)\Delta t\\ ={\sigma }^2\left(x_1,x_2,x_3\right)\Delta t\end{array}$

We obtain:

$\{\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=-\frac{\alpha S\left(t\right)I\left(t\right)}{N}+{\sigma }_{1,1}\frac{\text{d}W}{\text{d}t}+{\sigma }_{1,2}\frac{\text{d}W}{\text{d}t}\\ \frac{\text{d}E\left(t\right)}{\text{d}t}=\frac{\alpha S\left(t\right)I\left(t\right)}{N}-\beta E\left(t\right)+{\sigma }_{2,1}\frac{\text{d}W}{\text{d}t}+{\sigma }_{2,2}\frac{\text{d}W}{\text{d}t}+{\sigma }_{2,3}\frac{\text{d}W}{\text{d}t}\\ \frac{\text{d}I\left(t\right)}{\text{d}t}=\beta E\left(t\right)-\text{p}\gamma \text{I(t)}-\text{(1}-\text{p)}\lambda \text{I(t)}+{\sigma }_{3,2}\frac{\text{d}W}{\text{d}t}+{\sigma }_{3,3}\frac{\text{d}W}{\text{d}t}\end{array}$ (8)

where a Wienner i process.

2) For the sub-model “after infections”:

Note $\left(\begin{array}{c}{x}_4\\ x_5\end{array}\right)$ the random variables associated with sub-model AH.

For a small interval ∆t. Transition probabilities are:

$P\left(A+k,H+j\right)=\{\begin{array}{l}{p}_4=p\gamma I\Delta t,\left(k,j\right)=\left(1,0\right)\\ p_5=\left(1-p\right)\lambda I\Delta t,\left(k,j\right)=\left(0,1\right)\\ p_6=\left(\left(1-q\right)\delta H+q\theta H\right)\Delta t,\left(k,j\right)=\left(0,-1\right)\\ 1-\left({\displaystyle \underset{i=4}{\overset{6}{\sum }}{p}_i}\right)\Delta t,\left(k,j,m\right)=\left(0,0\right)\\ 0\text{otherwise}\end{array}$

$\begin{array}{c}E\left(\Delta x\right)={\displaystyle \underset{i=4}{\overset{6}{\sum }}{p}_i{\left(\Delta x\right)}_i}=p_4{\left(\Delta x\right)}_4+p_5{\left(\Delta x\right)}_5+p_6{\left(\Delta x\right)}_6\\ =p_4\left(\begin{array}{c}1\\ 0\end{array}\right)+p_5\left(\begin{array}{c}0\\ 1\end{array}\right)+p_6\left(\begin{array}{c}0\\ -1\end{array}\right)\end{array}$

$\begin{array}{c}E\left(\Delta x\right)=\left(\begin{array}{c}p\gamma I\\ \left(1-p\right)\lambda I-q\theta H-\left(1-q\right)\delta H\end{array}\right)\Delta t\\ =\mu \left(x_1,x_2,x_3\right)\Delta t\end{array}$

$\begin{array}{c}E\left(\left(\Delta x\right)\left(\Delta x\right)\right)={\displaystyle \underset{i=4}{\overset{6}{\sum }}{p}_i{\left(\Delta x\right)}_i{\left(\Delta x\right)}_i^{\text{T}}}\\ =p_4\left(\begin{array}{c}0\\ 1\end{array}\right)\left(\begin{array}{cc}0& 1\end{array}\right)+p_5\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{cc}1& 0\end{array}\right)+p_6\left(\begin{array}{c}-1\\ 0\end{array}\right)\left(\begin{array}{cc}-1& 0\end{array}\right)\\ =\left(\begin{array}{cc}{p}_5+p_6& 0\\ 0& p_4\end{array}\right)\\ =\left(\begin{array}{cc}\left(1-p\right)\lambda I+q\theta H+\left(1-q\right)\delta H& 0\\ 0& p\gamma I\end{array}\right)\Delta t\\ ={\sigma }^2\left(x_1,x_2,x_3\right)\Delta t\end{array}$.

We obtain:

$\{\begin{array}{l}\frac{\text{d}A\left(t\right)}{\text{d}t}=p\gamma I\left(t\right)+{\sigma }_{1,1}\frac{\text{d}{W}_1}{\text{d}t}\\ \frac{\text{d}H\left(t\right)}{\text{d}t}=\left(1-p\right)\lambda I\left(t\right)-\left(1-q\right)\delta H\left(t\right)-q\theta H\left(t\right)+{\sigma }_{2,2}\frac{\text{d}W}{\text{d}t}\end{array}$ (9)

where W the independent processes of Winner.

3) For the last sub-model (compartments after hospitalization):

Note $\left(\begin{array}{c}{x}_6\\ x_7\end{array}\right)$ random variables corresponding to the CD submodel. For a small inteval ∆t, the transition probabilities are:

$P\left(c+k,d+j\right)=\{\begin{array}{l}{p}_6=q\theta H\Delta t,\left(k,j\right)=\left(1,0\right)\\ p_7=\left(1-q\right)\delta I\Delta t,\left(k,j\right)=\left(0,1\right)\\ 1-\left({\displaystyle \underset{i=6}{\overset{7}{\sum }}{p}_i}\right)\Delta t,\left(k,j,m\right)=\left(0,0\right)\\ 0\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{otherwise}\end{array}$

$\begin{array}{c}E\left(\Delta x\right)={\displaystyle \underset{i=6}{\overset{7}{\sum }}{p}_i{\left(\Delta x\right)}_i}\\ =p_6\left(\begin{array}{c}1\\ 0\end{array}\right)+p_7(\; 0\; 1\; )\end{array}$

$\begin{array}{c}E\left(\Delta x\right)=\left(\begin{array}{c}q\theta H\\ 0\end{array}\right)\Delta t+\left(\begin{array}{c}0\\ \left(1-q\right)\delta H\end{array}\right)\Delta t\\ =\mu \left(x_1,x_2,x_3\right)\Delta t\end{array}$

$\begin{array}{c}E\left(\left(\Delta x\right)\left(\Delta x\right)\right)={\displaystyle \underset{i=6}{\overset{7}{\sum }}{p}_i{\left(\Delta x\right)}_i{\left(\Delta x\right)}_i^{\text{T}}}\\ =p_6\left(\begin{array}{c}1\\ 0\end{array}\right)\left(\begin{array}{cc}1& 0\end{array}\right)+p_7\left(\begin{array}{c}0\\ 1\end{array}\right)\left(\begin{array}{cc}0& 1\end{array}\right)\\ =\left(\begin{array}{cc}{p}_6& 0\\ 0& p_7\end{array}\right)\\ =\left(\begin{array}{cc}q\theta H& 0\\ 0& \left(1-q\right)\delta H\end{array}\right)\Delta t\\ ={\sigma }^2\left(x_1,x_2,x_3\right)\Delta t\end{array}$

The system becomes:

$\{\begin{array}{l}\frac{\text{d}C\left(t\right)}{\text{d}t}=q\theta H+{\sigma }_{1,1}\frac{\text{d}W}{\text{d}t}\\ \frac{\text{d}D\left(t\right)}{\text{d}t}=\left(1-q\right)\delta H+{\sigma }_{2,2}\frac{\text{d}W}{\text{d}t}\end{array}$ (10)

where W a Wienner process.

Figure 11 represents a simulation of the model through a numerical resolution of stochastic differential Equations (8)-(10).

5.2. Estimation of the Basic Reproduction Number R₀

1) Estimate by the method of maximum likelihood:

In [26], the authors propose maximum likelihood method to estimate the number of basic reproductions of an epidemic. This technique is based on the assumption that the number of secondary cases caused by an index case is a Poisson distribution, with the expected value R₀.

Let $\left(I_0,I_1,\cdots ,I_T\right)$ the observation on consecutive time and ω generation time distribution. The number of basic reproductions R₀ is estimated by a maximization of the log-likelihood:

$LL\left(R\right)={\displaystyle {\sum }_{t=1}^T\log \left(\frac{\exp \left(-{\mu }_t\right){\mu }_t^{I_t}}{I_t}\right)}$ (11)

where ${\mu }_t=R{\displaystyle {\sum }_{t-i}^t{I}_{t-i}{\omega }_i}$.

The likelihood function must be calculated over a period of exponential growth, and the R-square deviance can be used to select the best period. For Madagascar, using the observed data, the maximum likelihood estimate (MLE) of the basic reproduction number is:

Figure 11. Simulation of the stochastic equivalent of the modifed model SEIR.

$R_0=1.1027$

with 95% of confidence interval.

2) Bayesian estimate:

Introduced by Bettencourt and Ribeiro [27], this technique allows sequential estimate of the number of reproduction R. It is based on the assumption that the number of infections $I\left(t\right)$ is a discrete random variable and $I\left(t+1\right)$ is a Poisson random variable: $\mathcal{P}\left(\lambda \right)$, where

$\lambda =I\left(t\right)\exp \left(\gamma \left(R-1\right)\right)\text{et}\frac{1}{\gamma }$

the average duration of the infectious period.

With a non informative prior distribution of the basic number of reproduction R₀, we have implemented an algorithm that can calculate the posterior distribution:

$P\left(R_0|I_0,I_1,\cdots ,I_{t+1}\right)=\frac{P\left(I_{t+1}|R_0,I_0,I_1,\cdots ,I_t\right)P\left(R_0|I_0,I_1,\cdots ,I_t\right)}{P\left(I_0,\cdots ,I_{t+1}\right)}$ (12)

Knowledge about R₀ is updated as new data is observed. On each new day, the posterior distribution after the previous day is considered to be the prior distribition. At any given moment, the Bayesian estimate can be calculated with the highest probability density interval.

As before, the method requires that the epidemic be in a period of exponential growth, i.e., it does not take into account the depletion of sensitivity. It implicitly uses an exponential distribution for generation time; and assumes a random mix in the population [28].

For our case, on the graph of Figure 12, the grayed area corresponds to the

Figure 12. Estimated number of temporal reproductions until 2020-08-22 with a Bayesian approach.

confidence interval and the black curve to the average number of reproductions in time function $R\left(t\right)$. According to these results, the basic reproduction number in the Malagasy territory is lower than the references in neighboring countries [29] [30].

6. Conclusions

We have proposed a mathematical modelling framework for the spread of the COVID-19 epidemic in Madagascar taking into account local contexts. We developed a seven compartments model corresponding to the pandemic management strategies adopted by Malagasy authorities.

The implementation of containment and the imposition of various social distancing measures from the first outbreaks of the pandemic give a particular dynamic to the evolution of coronavirus in Madagascar. Due to insufficient screening tests, the detection rate of positive cases remains very low. There has been a sufficiently slow dynamic compared to other countries.

For the first deterministic model, the estimated basic reproduction number is close to 1.2. With the stochastic extensions of this model, we estimated values slightly below (near 1.1). These parameters suggest a low rate of spread of the epidemic during the first wave of coronavirus between March and August 2020 in Madagascar. According to official data, the epidemic appears to be under control. This deduction leaves much discussion given the lack of a massive test and the low rate of positive case detection.

Probabilistic tools provide a natural approach to modelling the evolution of the epidemic. Results already established in probability theory allow for more in-depth studies. The first stochastic model derived from a continuous-time Markov chain is reliable for a small-scale population. Managing transition time and events is a big challenge for a large number of populations. The stochastic differential equation (SDE) model is more interesting in this situation. Yet analyses of model characteristics, such as stability and convergence, are meriting further study.

For the proposed models, we were able to estimate the key parameters of epidemiological models such as the basic reproduction number. It can be interpreted as the measure of contagiousness and severity of the epidemic. Based on our results the Malagasy people are less vulnerable to the attack of the COVID-19. The epidemic is a little mild if we look at the official observed data.

With a Bayesian approach, we calculated the real-time estimation of the reproduction number. This time function of the $R0$ can be used to introduce control over the spread of the disease. This is also one of the perspectives of this work.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Anna O., Sandro G., David S., Carlo S. and the University Vita-Salute San Raffaele COVID-19 Literature Monitoring Working Group (2020) The First 10,000 COVID-19 Papers in Perspective: Are We Publishing What We Should Be Publishing? European Journal of Public Health, 30, 849-850. https://doi.org/10.1093/eurpub/ckaa170
[2]	Malik, V.S., Anton, K., Anne, P. and Boris, K. (2020) A New Modeling of the COVID 19 Pandemic. Chaos, Solitons & Fractals, 139, Article ID: 110039. https://doi.org/10.1016/j.chaos.2020.110039
[3]	Abdul, M.I. and Osei, A. (2020) Modelling Covid-19 Deaths in Ghana as a Discrete State Process in Continuous Time. American Journal of Applied Mathematics, 8, 344-355.
[4]	Raberahona, M., Rakotomalala, R., Rakotomijoro, E., Rahaingoalidera, T., Andry, C. E., Mamilaza, N., Razafindrabekoto, L., Rafanomezantsoa, E., Andriananja, V., Andrianasolo, R.L., Razafimahefa, S.H., Rakotoarivelo, R.A. and Randria, M. (2020) Clinical and Epidemiological Features Discriminating Confirmed COVID-19 Patients from SARS-CoV-2 Negative Patients at Screening Centers in Madagascar. International Journal of Infectious Diseases, 103, 6-8. https://doi.org/10.1016/j.ijid.2020.11.151
[5]	Faïçal, N., Iván, A., Juan, J.N. and Delfim, F.M. T. (2020) Mathematical Modeling of COVID-19 Transmission Dynamics with a Case Study of Wuhan. Chaos, Solitons & Fractals, 135, Article ID: 109846. https://doi.org/10.1016/j.chaos.2020.109846
[6]	Itai, C., Mathias, D., Mbuzeleni, H. and Tafadzwa, D. (2020) COVID-19: Unpacking the Low Number of Cases in Africa. Public Health in Practice, 1, Article ID: 100038. https://doi.org/10.1016/j.puhip.2020.100038
[7]	Marius, G., Giulia, P., Francesco, P., Eugenio, V., Chiara, P., Pierre-Yves, B., Eric, D., Yazdan, Y., Serge, P.E., Mathias, A., Bernardo, G., Moritz, U.G.K. and Vittoria, C. (2020) Preparedness and Vulnerability of African Countries against Importations of COVID-19: A Modelling Study. The Lancet, 395, 871-877. https://doi.org/10.1016/S0140-6736(20)30411-6
[8]	Klaus, D. and Heesterbeek, J.A.P. (2002) Daniel Bernoulli’s Epidemiological Model Revisited. Mathematical Biosciences, 180, 1-21. https://doi.org/10.1016/S0025-5564(02)00122-0
[9]	Mandal, S., Sarkar, R.R. and Sinha, S. (2011). Mathematical Models of Malaria—A Review. Malaria Journal, 10, Article NO. 202. https://doi.org/10.1186/1475-2875-10-202
[10]	Fred, B. (2017) Mathematical Epidemiology: Past, Present, and Future. Infectious Disease Modelling, 2, 113-127. https://doi.org/10.1016/j.idm.2017.02.001
[11]	Kyrychko, Y.N., Blyuss, K.B. and Brovchenko I. (2020) Mathematical Modelling of the Dynamics and Containment of COVID-19 in Ukraine. Scientific Reports, 10, Article No. 19662. https://doi.org/10.1038/s41598-020-76710-1
[12]	Piu, S., Jayanta, M. and Subhas, K. (2020) A Mathematical Model for COVID-19 Transmission Dynamics with a Case Study of India. Chaos, Solitons & Fractals, 140, Article ID: 110173. https://doi.org/10.1016/j.chaos.2020.110173
[13]	Nguyen, H.T., Hakimeh, M. and Shahram, R. (2020) A Mathematical Model for COVID-19 Transmission by Using the Caputo Fractional Derivative. Chaos, Solitons & Fractals, 140, Article ID: 110107. https://doi.org/10.1016/j.chaos.2020.110107
[14]	Allen, L.J.S. (2008) An Introduction to Stochastic Epidemic Models. In: Brauer, F., van den Driessche, P. and Wu, J., Eds., Mathematical Epidemiology, Lecture Notes in Mathematics, Vol. 1945, Springer, Berlin, Heidelberg, 81-130. https://doi.org/10.1007/978-3-540-78911-6_3
[15]	Manou-Abi, S. and Balicchi, J. (2020) Analysis of the COVID-19 Epidemic in French Overseas Department Mayotte Based on a Modified Deterministic and Stochastic SEIR Model. Cold Spring Harbor Laboratory Press.
[16]	Carcione, J.M., Santos, J.E., Bagaini, C. and Ba, J. (2020) A Simulation of a COVID-19 Epidemic Based on a Deterministic SEIR Mode. Frontiers in Public Health, 8, Article No. 230. https://doi.org/10.3389/fpubh.2020.00230
[17]	Kermack, W.O. and McKendrick, A.G. (1927) A Contribution to the Mathematical Theory of Epidemics. Proceedings of the Royal Society A, 115, 700-721. https://doi.org/10.1098/rspa.1927.0118
[18]	Norris, J. (1997) Markov Chains. Cambridge Series in Statistical and Probabilistic Mathematics, Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511810633
[19]	Witbooi, P.J. (2013) Stability of an SEIR Epidemic Model with Independent Stochastic Perturbations. Physica A: Statistical Mechanics and its Applications, 392, 4928-4936. https://doi.org/10.1016/j.physa.2013.06.025
[20]	Charles, J.M., Herman, E.G. and Nira, H. (1988) A Methodological Study of a Stochastic Model of an AIDS Epidemic. Mathematical Biosciences, 92, 201-229. https://doi.org/10.1016/0025-5564(88)90032-6
[21]	Cai, Y., Kang, Y. and Wang, W. (2017) A Stochastic SIRS Epidemic Model with Nonlinear Incidence Rate. Applied Mathematics and Computation, 305, 221-240. https://doi.org/10.1016/j.amc.2017.02.003
[22]	Byambasuren, O., Cardona, M., Bell, K., Clark, J., McLaws, M.L. and Glasziou, P. (2020) Estimating the Extent of Asymptomatic COVID-19 and Its Potential for Community Transmission: Systematic Review and Meta-Analysis. Cold Spring Harbor Laboratory Press.
[23]	Durrett, R. (2010) Brownian Motion. In: Durrett, R., Ed., Probability: Theory and Examples, Cambridge University Press, Cambridge, 353-400. https://doi.org/10.1017/CBO9780511779398.009
[24]	Ndanguza, D., Mbalawata, I. and Nsabimana, J. (2016) Analysis of SDEs Applied to SEIR Epidemic Models by Extended Kalman Filter Method. Applied Mathematics, 7, 2195-2211. http://dx.doi.org/10.4236/am.2016.717175
[25]	Allen, E.J., Allen, L.J.S., Armando, A. and Greenwood, P.E. (2008) Construction of Equivalent Stochastic Differential Equation Models. Stochastic Analysis and Applications, 26, 274-297. https://doi.org/10.1080/07362990701857129
[26]	White, L.F. and Pagano, M. (2008) A Likelihood-Based Method for Real-Time Estimation of the Serial Interval and Reproductive Number of an Epidemic. Statistics in Medicine, 27, 2999-3016. https://doi.org/10.1002/sim.3136
[27]	Bettencourt, L.M. and Ribeiro, R.M. (2008) Real Time Bayesian Estimation of the Epidemic Potential of Emerging Infectious Diseases. PLoS ONE, 3, e2185. https://doi.org/10.1371/journal.pone.0002185
[28]	Wang, K., Zhao, S., Li, H., Song, Y., Wang, L., Wang, M.H., Peng, Z., Li, H. and He, D. (2020) Real-Time Estimation of the Reproduction Number of the Novel Coronavirus Disease (COVID-19) in China in 2020 Based on Incidence Data. Annals of Translational Medicine, 8, Article No. 689. https://doi.org/10.21037/atm-20-1944
[29]	Evans, M.V., Garchitorena, A., Rakotonanahary, R.J., Drake, J.M., Andriamihaja, B., Rajaonarifara, E., Ngonghala, C.N., Roche, B., Bonds, M.H. and Rakotonirina, J. (2020) Reconciling Model Predictions with Low Reported Cases of COVID-19 in Sub-Saharan Africa: Insights from Madagascar. Global Health Action, 13, Article No. 1816044. https://doi.org/10.1080/16549716.2020.1816044
[30]	Musa, S.S., Zhao, S., Wang, M.H., et al. (2020) Estimation of Exponential Growth Rate and Basic Reproduction Number of the Coronavirus Disease 2019 (COVID-19) in Africa. Infectious Diseases of Poverty, 9, Article No. 96. https://doi.org/10.1186/s40249-020-00718-y