Stochastic Model of Dengue: Analysing the Probability of Extinction and LLN

Ragnimwendé Sawadogo, Fourtoua Victorien Konané, Wahabo Baguian

Département de Mathématiques, Unité de Recherche et de Formation en Sciences Exactes et Appliquées, Université Joseph KI-ZERBO, Ouagadougou, Burkina Faso.

DOI: 10.4236/am.2024.159036                

Abstract

In this article, we develop and analyze a continuous-time Markov chain (CTMC) model to study the resurgence of dengue. We also explore the large population asymptotic behavior of probabilistic model of dengue using the law of large numbers (LLN). Initially, we calculate and estimate the probabilities of dengue extinction and major outbreak occurrence using multi-type Galton-Watson branching processes. Subsequently, we apply the LLN to examine the convergence of the stochastic model towards the deterministic model. Finally, theoretical numerical simulations are conducted exploration to validate our findings. Under identical conditions, our numerical results demonstrate that dengue could vanish in the stochastic model while persisting in the deterministic model. The highlighting of the law of large numbers through numerical simulations indicates from what population size a deterministic model should be considered preferable.

Keywords

Dengue Fever, Continuous-Time Markov Chain, Multitype Branching Process, Probability of Disease Extinction, Law of Large Numbers

1. Introduction

Dengue, a viral disease spread by mosquitoes, has spread to become endemic worldwide [1]. The dengue virus is spread by female Aedes aegypti mosquitoes [2]. The dengue virus, belonging to the flavivirus genus, is divided into four serotypes: DENV-1, DENV-2, DENV-3, and DENV-4 [3] [4]. The danger of this disease lies primarily in its ability to affect nearly all age groups, from infants to adults, and to re-emerge rapidly in nearly all human populations [5]. This disease, after malaria, is considered one of the most dangerous [6]. The disease can also lead to significant economic losses, hindering development. For instance, during the 2011 epidemic, countries in Latin America and Asia experienced estimated losses of $12 million USD and $6.75 million USD respectively. In Thailand alone, tourist revenues plummeted by $363 million USD during the tourist season due to the epidemic [7] [8]. Although measures have been implemented to combat dengue fever [9]-[12], this disease remains without an effective treatment [13]. Therefore, understanding the transmission dynamics of this disease is of paramount importance, prompting researchers to focus on studying it. Proposals have been put forward for deterministic models of dengue fever [14]-[16], which, while providing informative results and forecasts, have limitations in predicting future dynamics. To enhance realism, research has turned to stochastic models associated with dengue fever [17]-[21]. Introducing random variations is a commonly used method to depict phenomena where a quantity is subject to constant but varying slight fluctuations. However, it’s important to note that Brownian motion cannot account for the effects of intense and sudden external disturbances such as climate change, floods, earthquakes, or tornadoes [22] [23]. This is why Driss Kiouach and colleagues proposed a stochastic model incorporating Lvy jumps to study the average extinction and persistence of dengue fever [5]. To the best of our knowledge, continuous-time Markov chain (CTMC) models for dengue are scarce. In this study, we develop a stochastic CTMC and multi-type branching process of Galton-Watson model. Continuous-time Markov chains (CTMC) are particularly suitable for modeling dengue, as transitions between the states of susceptible, infected, and recovered, as well as interactions between humans and mosquitoes, occur continuously and randomly over time. Additionally, the duration of the incubation and infectious phases can vary significantly from one person to another. Regarding multi-type Galton-Watson branching processes, they are used to represent successive generations of infections and to estimate the probability of epidemic extinction or the probability that the epidemic will persist and spread. This is crucial for understanding and anticipating the resurgence of dengue. Moreover, we use the Law of Large Numbers to demonstrate that stochastic models of dengue converge to their deterministic counterparts when the population is sufficiently large. This validates the use of deterministic models for large populations in the context of dengue.

The article is structured as follows: In Section 2, we review the definitions of the parameters of the base model. Section 3 focuses on the analysis of the proposed CTMC model. In Section 4, we present the probabilistic model and highlight the law of large numbers. Finally, Section 5 comprises numerical simulations of the proposed model to evaluate the results.

2. Recall of the Baseline Deterministic Model

For the convenience of the reader, we briefly recap the main results of the baseline model [24]. The compartments $\overline{S}\mathrm{,}\overline{I}$ and $\overline{R}$ represent respectively the total number of susceptible, infected, and recovered individuals during the epidemic. The acronyms ${\mathcal{N}}_H$ and ${\mathcal{N}}_m$ denote the population sizes of humans and mosquitoes respectively. The system of ordinary differential equations of the model is given by:

$\{\begin{array}{l}\frac{\text{d}{\overline{\mathcal{S}}}_H\left(t\right)}{dt}={\mu }_H{\mathcal{N}}_H-{\lambda }_H{\overline{\mathcal{S}}}_H\left(t\right)\frac{{\overline{ℐ}}_m\left(t\right)}{{\mathcal{N}}_m}-{\mu }_H{\overline{\mathcal{S}}}_H\left(t\right),\hfill \\ \frac{\text{d}{\overline{ℐ}}_H\left(t\right)}{\text{d}t}={\lambda }_H{\overline{\mathcal{S}}}_H\left(t\right)\frac{{\overline{ℐ}}_m\left(t\right)}{{\mathcal{N}}_m}-\left({\mu }_H+{\gamma }_H+{\alpha }_H\right){\overline{ℐ}}_H\left(t\right),\hfill \\ \frac{\text{d}{\overline{ℛ}}_H\left(t\right)}{\text{d}t}={\gamma }_H{\overline{ℐ}}_H\left(t\right)-{\mu }_H{\overline{ℛ}}_H\left(t\right),\hfill \\ \frac{\text{d}{\overline{\mathcal{S}}}_m\left(t\right)}{\text{d}t}={\mu }_m{\mathcal{N}}_m-{\lambda }_m{\overline{\mathcal{S}}}_m\left(t\right)\frac{{\overline{ℐ}}_H\left(t\right)}{{\mathcal{N}}_H}-{\mu }_m{\overline{\mathcal{S}}}_m\left(t\right),\hfill \\ \frac{\text{d}{\overline{ℐ}}_m\left(t\right)}{\text{d}t}={\lambda }_m{\overline{\mathcal{S}}}_m\left(t\right)\frac{{\overline{ℐ}}_H\left(t\right)}{{\mathcal{N}}_H}-{\mu }_m{\overline{ℐ}}_m\left(t\right).\hfill \end{array}$(1)

The description of parameters in model system (1) is given by the following Table 1.

Table 1. Model parameters and their interpretations.

Name	Description	Value
λH	is the actual contact rate between susceptible humans and mosquitoes	0.5
γH	the recovery rate of humans from dengue	0.4
αH	represents the death rate of humans induced by dengue	0.1
μH	is the natural mortality rate of humans	0.2
λm	is the actual contact rate between susceptible mosquitoes and humans	0.21
μm	is the natural mortality rate of mosquitoes	0.2

Let us now introduce the proportions

$\begin{array}{l}{S}_H\left(t\right)=\frac{{\overline{\mathcal{S}}}_H\left(t\right)}{{\mathcal{N}}_H},I_H\left(t\right)=\frac{{\overline{ℐ}}_H\left(t\right)}{{\mathcal{N}}_H},R_H\left(t\right)=\frac{{\overline{ℛ}}_H\left(t\right)}{{\mathcal{N}}_H},\\ S_m\left(t\right)=\frac{{\overline{\mathcal{S}}}_m\left(t\right)}{{\mathcal{N}}_m},I_m\left(t\right)=\frac{{\overline{ℐ}}_m\left(t\right)}{{\mathcal{N}}_m}.\end{array}$ (2)

Then, we obtain the following equalities $S_H\left(t\right)+I_H\left(t\right)+R_H\left(t\right)=1$ and $S_m\left(t\right)+I_m\left(t\right)=1$ .

3. Stochastic Model for Dengue

3.1. Continuous Time Markov Chain Model Formulation

In this section, we formulate the CTMC (Continuous-Time Markov Chain) model of the deterministic model developed in [24]. To simplify, we retain the same notations for the random variables and parameters as those used in the deterministic model. The variables ${\overline{S}}_H$ , ${\overline{I}}_H$ , ${\overline{R}}_H$ , ${\overline{S}}_m$ and ${\overline{I}}_m$ are now discrete, and t is in $\left[\mathrm{0,}\infty \right)$ . Set $x=\left(s_H\mathrm{,}{i}_H\mathrm{,}{r}_H\mathrm{,}{s}_m\mathrm{,}{i}_m\right)$ and $y=\left(s_H+{\ell }_1\mathrm{,}{i}_H+{\ell }_2\mathrm{,}{r}_H+{\ell }_3\mathrm{,}{s}_m+{\ell }_4\mathrm{,}{i}_m+{\ell }_5\right)$ . The transition probabilities associated with the model (1) are given by

(3)

where,

$\begin{array}{c}ℒ\left(t\right)={\lambda }_H{\overline{\mathcal{S}}}_H\left(t\right)\frac{{\overline{ℐ}}_m\left(t\right)}{N_m}+{\mu }_H\left({\overline{\mathcal{S}}}_H+{\overline{ℐ}}_H+{\overline{ℛ}}_H\right)+{\mu }_H{\mathcal{N}}_H\\ +{\lambda }_m{\overline{\mathcal{S}}}_m\left(t\right)\frac{{\overline{ℐ}}_H\left(t\right)}{N_H}+{\mu }_m\left({\overline{\mathcal{S}}}_m+{\overline{ℐ}}_m\right)+\left({\alpha }_H+{\gamma }_H\right){\overline{ℐ}}_H+{\mu }_m{\mathcal{N}}_m\mathrm{.}\end{array}$

We choose $\Delta t$ small enough so that all these quantities are between 0 and 1. The associated matrix is stochastic.

3.2. The Bienaymé-Galton-Watson Branching Process (BGWbp)

Here, we recall some concepts about branching process theory and apply it to our model to find the disease invasion and extinction probabilities. We start by defining a Bienaymé-Galton-Watson branching process as described in [25]-[27].

Definition 3.1. [27] A multitype BGWbp $I{\left(t\right)}_{t=0}^{\infty }$ is a collection of vector random variables $I\left(t\right)$ , where each vector consists of l different types, $I\left(t\right)=\left(I_1\left(t\right)\mathrm{,}{I}_2\left(t\right)\mathrm{,}\cdots \mathrm{,}{I}_l\left(t\right)\right)$ and each random variable $I_i\left(t\right)$ has l associated offspring random variables for the number of offsprings of type $j=\mathrm{1,2,}\cdots \mathrm{,}l$ from a parent of type i.

The offspring probability generating function (pgf) $g_i\mathrm{<mo>:</mo>}{\left[\mathrm{0,1}\right]}^n\to \left[\mathrm{0,1}\right]$ , for type i given $I_i\left(0\right)=1$ and $I_j\left(0\right)=0$ , $j\ne i$ , is defined as

$g_i\left({\omega }_1,\cdots ,{\omega }_n\right)={\displaystyle \sum _{l_n=0}^{\infty }\cdots }{\displaystyle \sum _{l_1=0}^{\infty }{\text{P}}_i\left(l_1,\cdots ,l_n\right){\omega }_1^{l_1}\cdots {\omega }_n^{l_n}}$ (4)

where ${\text{P}}_i\left(l_1\mathrm{,}\cdots \mathrm{,}{l}_n\right)=\text{Prob}\left\{Z_{1i}=l_1\mathrm{,}\cdots \mathrm{,}{Z}_{ni}=l_n\right\}$ is the probability that one infected individual of type i gives birth to $l_j$ individuals of type j and there is always a fixed point at $g_i\left(\mathrm{1,1,}\cdots \mathrm{,1}\right)=1$ . $g_i\left(\mathrm{0,0,}\cdots \mathrm{,0}\right)$ denotes the probability of extinction for $I_i$ given that $I_i\left(0\right)=1$ and $I_j\left(0\right)=0$ for all other types.

We define the expectation matrix $M=\left[m_{ji}\right]$ as an $n\times n$ , non negative and irreducible matrix where the entry $m_{ji}$ is the expected number of offsprings of individuals of type j produced by an infective individual of type i. The elements of matrix M are calculated from Equation (4) by differentiating $g_i$ with respect to ${\omega }_j$ and then evaluating all the $\omega$ variables at 1, that is,

$m_{ji}={\frac{\partial g_i}{\partial {\omega }_j}|}_{\omega =1}.$ (5)

See [28]-[30] for details. The following theorem summarizes the conditions for the probability of extinction or persistence of dengue.

Theorem 3.1. Let the initial sizes for each type be $I_i\left(0\right)=i_i$ , $i=1,2,\cdots ,l$ . Consider the generating functions $g_i$ for each of the l types are non-linear functions of ${\omega }_j$ with some $g_i\left(\mathrm{0,0,}\cdots \mathrm{,0}\right)>0$ and suppose that the expectation matrix $M=\left[m_{ji}\right]$ is an $ı\times n$ non negative and irreducible matrix, and $\rho \left(M\right)$ is the spectral radius of matrix M.

1) If $\rho \left(M\right)<1$ or $\rho \left(M\right)=1$ (sub critical and critical case respectively), then the probability of ultimate extinction is one:

$\underset{t\to \infty }{\lim }\text{Prob}\left\{I\left(t\right)=0\right\}=1.$ (6)

2) If $\rho \left(M\right)>1$ (supercritical case), then the probability of ultimate disease extinction is less than one:

$\underset{t\to \infty }{\lim }\text{Prob}\left\{I\left(t\right)=0\right\}={\omega }_1^{i_1}{\omega }_2^{i_2}\cdots q_l^{i_k}<1,$ (7)

where $\left({\omega }_1\mathrm{,}{\omega }_2\mathrm{,}\cdots \mathrm{,}{q}_l\right)$ is the unique fixed point of the k offspring pgf, $f_i\left({\omega }_1\mathrm{,}{\omega }_2\mathrm{,}\cdots \mathrm{,}{q}_l\right)=q_i$ and $0<q_i<1$ , $i=1,2,\cdots ,l$ [31]. The value of $q_i$ is the probability of disease extinction for infectious of type i and the probability of an outbreak is approximately

$1-{\omega }_1^{i_1}{\omega }_2^{i_2}\cdots q_l^{i_l}\mathrm{.}$

3.3. Extinction Probabilities for Dengue

In stochastic epidemic theory, it is feasible to predict the onset and cessation of a disease based on the initial number of infected individuals. According to Lloyd et al. [32], an ${ℛ}_0$ greater than one does not guarantee the persistence of the infection in a fully susceptible population. Unlike deterministic models, stochastic models demonstrate that even with a low initial number of infected individuals, disease extinction is possible. In such cases, stochastic models predict a minor epidemic, whereas deterministic models consistently predict a major epidemic. At the beginning of the epidemic, the decrease in the number of susceptible individuals is negligible, allowing for the estimation of invasion probabilities using a linear model. This approach, assuming that the entire population is susceptible [33] [34], often employs the theory of Galton-Watson branching processes with multiple types to calculate these probabilities [33]. In this theory, individuals in the population are classified into different types and act independently of each other.

The rest of this section, the probabilities of extinction and epidemic outbreaks are determined using probability generating functions, where we seek the fixed points of the system for each infectious type. We assume that the susceptible populations are at disease-free equilibrium (DFE), and the branching process is linear near this equilibrium and time-homogeneous. The stochastic threshold can be approximated using a two-type BGWbp. Infected humans are referred to as type 1, and mosquitoes as type 2. Therefore, the offspring pgf for type 1 humans is given by

$\begin{array}{c}{g}_1\left({\omega }_1\mathrm{,}{\omega }_2\right)=P_2\left(\mathrm{1,1}\right){\omega }_1{\omega }_2+P_2\left(\mathrm{0,0}\right)\\ =\frac{{\lambda }_H{\omega }_1{\omega }_2+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)}{{\lambda }_H+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)}\mathrm{,}{\omega }_1\mathrm{,}{\omega }_2\in \left[\mathrm{0,1}\right]\end{array}$

Likewise, in the case of type 2 mosquitoes, we have

$\begin{array}{c}{g}_2\left({\omega }_1\mathrm{,}{\omega }_2\right)=P_2\left(\mathrm{1,1}\right){\omega }_1{\omega }_2+P_2\left(\mathrm{0,0}\right)\\ =\frac{{\lambda }_m{\omega }_1{\omega }_2+{\mu }_m}{{\lambda }_m+{\mu }_m}\mathrm{,}{\omega }_1\mathrm{,}{\omega }_2\in \left[\mathrm{0,1}\right]\end{array}$

Now we can compute the expectation matrix using the offspring probabilities. The individual elements of the matrix, denoted as $m_{ji}$ , are given by

$\begin{array}{l}{m}_{11}={\frac{\partial g_1\left({\omega }_1,{\omega }_2\right)}{\partial {\omega }_1}|}_{{\omega }_1=1,{\omega }_2=1}=\frac{{\lambda }_H}{{\lambda }_H+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)},\\ m_{12}={\frac{\partial g_2\left({\omega }_1,{\omega }_2\right)}{\partial {\omega }_1}|}_{{\omega }_1=1,{\omega }_2=1}=\frac{{\lambda }_m}{{\lambda }_m+{\mu }_m},\\ m_{21}={\frac{\partial g_1\left({\omega }_1,{\omega }_2\right)}{\partial {\omega }_2}|}_{{\omega }_{1=1,{\omega }_2=1}}=\frac{{\lambda }_H}{{\lambda }_H+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)},\\ m_{22}={\frac{\partial g_2\left({\omega }_1,{\omega }_2\right)}{\partial {\omega }_2}|}_{{\omega }_1=1,{\omega }_2=1}=\frac{{\lambda }_m}{{\lambda }_m+{\mu }_m}\end{array}$ (8)

and therefore, the expectation matrix is given by

$M=\left[\begin{array}{cc}\mathcal{A}& ℬ\\ \mathcal{A}& ℬ\end{array}\right]$ (9)

where

$\mathcal{A}=\frac{{\lambda }_H}{{\lambda }_H+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)}\text{and}ℬ=\frac{{\lambda }_m}{{\lambda }_m+{\mu }_m}\mathrm{.}$ (10)

The eigenvalues of matrix M are the roots of the characteristic equation

$X^2-\left(\mathcal{A}+ℬ\right)X=0.$ (11)

The spectral radius $\rho \left(M\right)$ of matrix M obtained from Equation (11) is given by

$\rho \left(M\right)=\mathcal{A}+ℬ=\frac{\mathcal{N}}{\mathcal{D}}\mathrm{,}$ (12)

where

$\mathcal{N}={\lambda }_H\left({\lambda }_m{\mathcal{N}}_m+{\mu }_m\right)+{\lambda }_m\left({\lambda }_H+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)\right)\mathrm{,}$

$\mathcal{D}=\left({\lambda }_H+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)\right)\left({\lambda }_m+{\mu }_m\right)\mathrm{.}$

The probability of ultimate disease extinction is one if $\rho \left(M\right)<1$ . This means

${\lambda }_H{\lambda }_m<\left({\mu }_H+{\gamma }_H+{\alpha }_H\right){\mu }_m$ (13)

which reduces to

$\frac{{\lambda }_H{\lambda }_m}{\left({\mu }_H+{\gamma }_H+{\alpha }_H\right){\mu }_m}<1\iff {ℛ}_0^2<1.$ (14)

In the context of disease dynamics, especially for diseases where the basic reproductive number ${ℛ}_0>1$ (or $\rho \left(M\right)>1$ ), there is a positive probability of a major epidemic occurring. This probability is associated with a fixed point of the offspring probability generating functions (pgfs) on the interval ${\left(\mathrm{0,1}\right)}^2$ .

To find this fixed point, denoted as $\left({\omega }_1\mathrm{,}{\omega }_2\right)$ , we solve the following system of equations $g_1\left({\omega }_1\mathrm{,}{\omega }_2\right)={\omega }_1$ and $g_2\left({\omega }_1\mathrm{,}{\omega }_2\right)={\omega }_2$ .

Here, $g_1\left({\omega }_1\mathrm{,}{\omega }_2\right)$ and $g_2\left({\omega }_1\mathrm{,}{\omega }_2\right)$ represent the offspring probability generating functions (pgfs) corresponding to infected human and mosquito populations, respectively. These functions describe the probabilities of ultimate disease extinction for each population. The fixed point $\left({\omega }_1\mathrm{,}{\omega }_2\right)$ gives us the probabilities that the infected human and mosquito populations will eventually die out.

The trivial fixed point $\left(\mathrm{1,1}\right)$ (where both ${\omega }_1=1$ and ${\omega }_2=1$ ) always represents complete extinction of both infected populations. However, for ${ℛ}_0>1$ , there may exist another non-trivial fixed point $\left({\omega }_1\mathrm{,}{\omega }_2\right)$ where ${\omega }_1\mathrm{,}{\omega }_2\in \left(\mathrm{0,1}\right)$ . This non-trivial fixed point corresponds to a scenario where the disease persists in the population. So, let’s consider that $\left({\omega }_1\mathrm{,}{\omega }_2\right)\ne \left(\mathrm{1,1}\right)$ and solve the following system of equations.

$\{\begin{array}{l}\frac{{\lambda }_H{\omega }_1{\omega }_2+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)}{{\lambda }_H+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)}={\omega }_1,\\ \frac{{\lambda }_m{\omega }_1{\omega }_2+{\mu }_m}{{\lambda }_m+{\mu }_m}={\omega }_2.\end{array}$ (15)

Expressing ${\omega }_2$ in terms of ${\omega }_1$ in Equation (15) gives

${\omega }_2=\frac{{\mu }_m}{{\lambda }_m\left(1-{\omega }_1\right)+{\mu }_m}\mathrm{.}$ (16)

Substituting the Equation (16) into Equation (15) and simplifying, we derive the quadratic equation.

${\mathcal{A}}_1{\omega }_1^2+{\mathcal{A}}_2{\omega }_1+{\mathcal{A}}_3=0\text{,}$ (17)

where

$\begin{array}{l}{\mathcal{A}}_1={\lambda }_m\left({\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H\right)\mathrm{,}\\ {\mathcal{A}}_2={\lambda }_H\left({\mu }_m\right)-{\lambda }_m\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)\\ -\left({\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H\right)\left({\lambda }_m+{\mu }_H+{\gamma }_H+{\alpha }_H\right)\mathrm{,}\\ {\mathcal{A}}_3=\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)\left({\lambda }_m+{\mu }_m\right)\mathrm{.}\end{array}$ (18)

Solving Equation (17) for ${\omega }_1$ and then substituting its expression in Equation (15) yields the following expressions for ${\omega }_1$ and ${\omega }_2$ :

$\begin{array}{l}{\omega }_1=\frac{\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)\left({\lambda }_m+{\mu }_m\right)}{{\lambda }_m\left({\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H\right)}\mathrm{,}\\ {\omega }_2=\frac{\left({\mu }_m\right)\left({\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H\right)}{{\lambda }_H\left({\lambda }_m+{\mu }_m\right)}\mathrm{.}\end{array}$ (19)

We express ${\omega }_1$ and ${\omega }_2$ in terms of the basic reproduction number (11) to obtain

$\begin{array}{l}{\omega }_1=\frac{{\lambda }_H}{{\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H}\left(\frac{1}{{ℛ}_0^2}\right)+\frac{{\mu }_H+{\gamma }_H+{\alpha }_H}{{\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H}\mathrm{,}\\ {\omega }_2=\frac{{\lambda }_m}{{\lambda }_m+{\mu }_m}\left(\frac{1}{{ℛ}_0^2}\right)+\frac{{\mu }_m}{{\lambda }_m+{\mu }_m}\mathrm{.}\end{array}$ (20)

The probability ${\omega }_1$ can be interpreted epidemiologically as follows: an infected human will transmit the disease to a susceptible mosquito with a probability ${\lambda }_H/\left({\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H\right)$ or die or recover before transmitting the disease with probability $\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)/\left({\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H\right)$ . Similarly, the probability ${\omega }_2$ has the following interpretation: an infected mosquito will transmit the disease to a susceptible human with a probability ${\lambda }_m/\left({\lambda }_m+{\mu }_m\right)$ or die before transmitting the disease with probability ${\mu }_m/\left({\lambda }_m+{\mu }_m\right)$ . If transmission of the disease from the mosquito to the human is successful, then the probability that the infected mosquito transmits the disease to another susceptible mosquito is $\left(1/{ℛ}_0^2\right)$ .

We calculate the probability of disease extinction using ${\omega }_1$ and ${\omega }_2$ . If $Y\left(0\right)=y_0$ and $X\left(0\right)=x_0$ represent the initial sizes of infected humans and infected mosquitoes respectively, then the probability of disease extinction is approximate

${\mathbb{P}}_{ext}={\omega }_1^{y_0}{\omega }_2^{x_0}={\left[\frac{{\lambda }_H+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right){ℛ}_0^2}{\left({\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H\right){ℛ}_0^2}\right]}^{y_0}{\left[\frac{{\lambda }_m+{\mu }_m{ℛ}_0^2}{\left({\lambda }_m+{\mu }_m\right){ℛ}_0^2}\right]}^{x_0}\mathrm{.}$ (21)

Therefore, the probability of a major disease outbreak $P_{ps}$ is expressed as follows.

$P_{ps}=1-{\mathbb{P}}_{ext}=1-{\left[\frac{{\lambda }_H+\left({\mu }_H+{\gamma }_H+{\alpha }_H\right){ℛ}_0^2}{\left({\lambda }_H+{\mu }_H+{\gamma }_H+{\alpha }_H\right){ℛ}_0^2}\right]}^{y_0}{\left[\frac{{\lambda }_m+{\mu }_m{ℛ}_0^2}{\left({\lambda }_m+{\mu }_m\right){ℛ}_0^2}\right]}^{x_0}\mathrm{.}$ (22)

4. Law of Large Numbers

Following the same principles as those presented by É. Pardoux in [35], we propose a probabilistic model of dengue. We denote by $P_i$ , $i=\overline{\mathrm{2,7}}$ and $i=\overline{\mathrm{9,11}}$ standard mutually independent Poisson processes. Let’s suppose that each death coincides with a birth and set

$S_H\left(t\right)=\frac{{\overline{\mathcal{S}}}_H\left(t\right)}{{\mathcal{N}}_H}\mathrm{,}$

$I_H\left(t\right)=\frac{{\overline{ℐ}}_H\left(t\right)}{{\mathcal{N}}_H}\mathrm{,}$

$R_H\left(t\right)=\frac{{\overline{ℛ}}_H\left(t\right)}{{\mathcal{N}}_H}\mathrm{,}$

$S_m\left(t\right)=\frac{{\overline{\mathcal{S}}}_m\left(t\right)}{{\mathcal{N}}_m}\mathrm{.}$

The evolution system of $\left(S_H\left(t\right)\mathrm{,}{I}_H\left(t\right)\mathrm{,}{R}_H\left(t\right)\mathrm{,}{S}_m\left(t\right)\mathrm{,}{I}_m\left(t\right)\right)$ becomes

$\{\begin{array}{l}{S}_H\left(t\right)=S_H\left(0\right)+\frac{1}{{\mathcal{N}}_H}{P}_3\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)\text{d}r\right)+\frac{1}{{\mathcal{N}}_H}{P}_6\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right)+\frac{1}{{\mathcal{N}}_H}{P}_7\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{R}_H\left(r\right)\text{d}r\right)\\ -\frac{1}{{\mathcal{N}}_H}{P}_2\left({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_H}{P}_3\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)\text{d}r\right),\\ I_H\left(t\right)=-\frac{1}{{\mathcal{N}}_H}{P}_5\left({\alpha }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_H}{P}_6\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right)+I_H\left(0\right)+\frac{1}{{\mathcal{N}}_H}{P}_2\left({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right)\\ -\frac{1}{{\mathcal{N}}_H}{P}_4\left({\gamma }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right),\\ R_H\left(t\right)=R_H\left(0\right)+\frac{1}{{\mathcal{N}}_H}{P}_4\left({\gamma }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_H}{P}_7\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{R}_H\left(r\right)\text{d}r\right),\\ S_m\left(t\right)=S_m\left(0\right)+\frac{1}{{\mathcal{N}}_m}{P}_{10}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)\text{d}r\right)+P_{11}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{I}_m\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_m}{P}_9\left({\lambda }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)I_H\left(r\right)\text{d}r\right)\\ -\frac{1}{{\mathcal{N}}_m}{P}_{10}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)\text{d}r\right),\\ I_m\left(t\right)=I_m\left(0\right)+\frac{1}{{\mathcal{N}}_m}{P}_9\left({\lambda }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)I_H\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_m}{P}_{11}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{I}_m\left(r\right)\text{d}r\right).\end{array}$ (23)

The processes $P_i$ , $i\mathrm{<mo>=</mo>}\overline{\mathrm{2,7}}$ and $i\mathrm{<mo>=</mo>}\overline{\mathrm{9,11}}$ being standard mutually independent Poisson processes, we define $M_i\left(t\right)=P_i\left(t\right)-t$ , $i=\overline{\mathrm{2,7}}$ and $i=\overline{\mathrm{9,11}}$ .

Hence

$\{\begin{array}{l}{S}_H\left(t\right)=S_H\left(0\right)+\frac{1}{{\mathcal{N}}_H}{M}_3\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)\text{d}r\right)+\frac{1}{{\mathcal{N}}_H}{M}_6\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right)+\frac{1}{{\mathcal{N}}_H}{M}_7\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{R}_H\left(r\right)\text{d}r\right)\\ -\frac{1}{{\mathcal{N}}_H}{M}_2\left({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_H}{M}_3\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)\text{d}r\right)+{\chi }_1\left(t\right),\\ I_H\left(t\right)=I_H\left(0\right)+\frac{1}{{\mathcal{N}}_H}{M}_2\left({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_H}{M}_4\left({\gamma }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right)\\ -\frac{1}{{\mathcal{N}}_H}{M}_5\left({\alpha }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_H}{M}_6\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right)+{\chi }_2\left(t\right),\\ R_H\left(t\right)=R_H\left(0\right)+\frac{1}{{\mathcal{N}}_H}{M}_4\left({\gamma }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_H}{M}_7\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{R}_H\left(r\right)\text{d}r\right)+{\chi }_3\left(t\right),\\ S_m\left(t\right)=S_m\left(0\right)+\frac{1}{{\mathcal{N}}_m}{M}_{10}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)\text{d}r\right)+M_{11}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{I}_m\left(r\right)\text{d}r\right)\\ -\frac{1}{{\mathcal{N}}_m}{M}_9\left({\lambda }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)I_H\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_m}{M}_{10}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)\text{d}r\right)+{\chi }_4\left(t\right),\\ I_m\left(t\right)=I_m\left(0\right)+\frac{1}{{\mathcal{N}}_m}{M}_9\left({\lambda }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)I_H\left(r\right)\text{d}r\right)-\frac{1}{{\mathcal{N}}_m}{M}_{11}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{I}_m\left(r\right)\text{d}r\right)+{\chi }_5\left(t\right),\end{array}$ (24)

where

$\{\begin{array}{l}{\chi }_1\left(t\right)={\mu }_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)\text{d}r+{\mu }_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r+{\mu }_H{\displaystyle {\int }_0^t}{R}_H\left(r\right)\text{d}r\\ -{\lambda }_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r-{\mu }_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)\text{d}r,\\ {\chi }_2\left(t\right)={\lambda }_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r-\left({\mu }_H+{\gamma }_H+{\alpha }_H\right){\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r,\\ {\chi }_3\left(t\right)={\gamma }_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r-{\mu }_H{\displaystyle {\int }_0^t}{R}_H\left(r\right)\text{d}r,\\ {\chi }_4\left(t\right)={\mu }_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)\text{d}r+{\mu }_m{\displaystyle {\int }_0^t}{I}_m\left(r\right)\text{d}r-{\lambda }_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)I_H\left(r\right)\text{d}r\\ -{\mu }_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)\text{d}r,\\ {\chi }_5\left(t\right)={\lambda }_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)I_H\left(r\right)\text{d}r-{\mu }_m{\displaystyle {\int }_0^t}{I}_m\left(r\right)\text{d}r.\end{array}$ (25)

Let’s consider the processes

$\begin{array}{l}{ℳ}_2\left(t\right)=\frac{1}{{\mathcal{N}}_H}{M}_2\left({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right),\\ {ℳ}_3\left(t\right)=\frac{1}{{\mathcal{N}}_H}{M}_3\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)\text{d}r\right),\\ {ℳ}_4\left(t\right)=\frac{1}{{\mathcal{N}}_H}{M}_4\left({\gamma }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right),\\ {ℳ}_5\left(t\right)=\frac{1}{{\mathcal{N}}_H}{M}_5\left({\alpha }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right),\\ {ℳ}_6\left(t\right)=\frac{1}{{\mathcal{N}}_H}{M}_6\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right),\\ {ℳ}_7\left(t\right)=\frac{1}{{\mathcal{N}}_H}{M}_7\left({\mu }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{R}_H\left(r\right)\text{d}r\right),\\ {ℳ}_9\left(t\right)=\frac{1}{{\mathcal{N}}_m}{M}_9\left({\lambda }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)I_H\left(r\right)\text{d}r\right),\\ {ℳ}_{10}\left(t\right)=\frac{1}{{\mathcal{N}}_m}{M}_{10}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)\text{d}r\right),\\ {ℳ}_{11}\left(t\right)=\frac{1}{{\mathcal{N}}_m}{M}_{11}\left({\mu }_m{\mathcal{N}}_m{\displaystyle {\int }_0^t}{I}_m\left(r\right)\text{d}r\right).\end{array}$ (26)

Let

${ℱ}_t=\sigma \left\{S_H\left(r\right)\mathrm{,}{I}_H\left(r\right)\mathrm{,}{R}_H\left(r\right)\mathrm{,}{S}_m\left(r\right)\mathrm{,}{I}_m\left(r\right)\mathrm{,0}\le r\le t\right\}\mathrm{.}$

Lemma 4.1. The $\left\{{ℳ}_i\left(t\right)\mathrm{,}t\ge 0\right\}$ , are a ${ℱ}_t$ -martingale which satisfy $\mathbb{E}\left[{ℳ}_i\left(t\right)\right]=0$ , $i=\overline{\mathrm{2,7}}$ , $i=\overline{\mathrm{9,11}}$ and

$\begin{array}{l}\mathbb{E}\left[{\left|{ℳ}_2\left(t\right)\right|}^2\right]=\frac{1}{{\mathcal{N}}_H}\mathbb{E}\left({\lambda }_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right),\\ \mathbb{E}\left[{\left|{ℳ}_3\left(t\right)\right|}^2\right]=\frac{1}{{\mathcal{N}}_H}\mathbb{E}\left({\mu }_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)\text{d}r\right),\\ \mathbb{E}\left[{\left|{ℳ}_4\left(t\right)\right|}^2\right]=\frac{1}{{\mathcal{N}}_H}\mathbb{E}\left({\gamma }_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right),\\ \mathbb{E}\left[{\left|{ℳ}_5\left(t\right)\right|}^2\right]=\frac{1}{{\mathcal{N}}_H}\mathbb{E}\left({\alpha }_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right),\end{array}$

$\begin{array}{l}\mathbb{E}\left[{\left|{ℳ}_6\left(t\right)\right|}^2\right]=\frac{1}{{\mathcal{N}}_H}\mathbb{E}\left({\mu }_H{\displaystyle {\int }_0^t}{I}_H\left(r\right)\text{d}r\right),\\ \mathbb{E}\left[{\left|{ℳ}_7\left(t\right)\right|}^2\right]=\frac{1}{{\mathcal{N}}_H}\mathbb{E}\left({\mu }_H{\displaystyle {\int }_0^t}{R}_H\left(r\right)\text{d}r\right),\\ \mathbb{E}\left[{\left|{ℳ}_8\left(t\right)\right|}^2\right]=\frac{1}{{\mathcal{N}}_m}\mathbb{E}\left({\lambda }_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)I_H\left(r\right)\text{d}r\right),\\ \mathbb{E}\left[{\left|{ℳ}_9\left(t\right)\right|}^2\right]=\frac{1}{{\mathcal{N}}_m}\mathbb{E}\left({\mu }_m{\displaystyle {\int }_0^t}{S}_m\left(r\right)\text{d}r\right),\\ \mathbb{E}\left[{\left|{ℳ}_{11}\left(t\right)\right|}^2\right]=\frac{1}{{\mathcal{N}}_m}\mathbb{E}\left({\mu }_m{\displaystyle {\int }_0^t}{I}_m\left(r\right)\text{d}r\right).\end{array}$ (27)

Proof. Let’s demonstrate the case of $\left\{{ℳ}_2\left(t\right)\right\}$ ; the remainder follows similarly. The martingale property of $\left\{{ℳ}_2\left(t\right)\right\}$ follows from the fact that for all $0<r<t$

$\begin{array}{l}\mathbb{E}\left[P_2\left({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right)-P_2\left({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^r}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right)|{ℱ}_r\right]\\ ={\lambda }_H{\mathcal{N}}_H\mathbb{E}\left[{\displaystyle {\int }_r^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r|{ℱ}_r\right].\end{array}$ (28)

□

We now establish that identity. For $n\ge \mathrm{1,}0\le v\le t$ , let

${\left[v\right]}_n=\{\begin{array}{ll}v\mathrm{,}\hfill & \text{if}v\le r\mathrm{<mo>;</mo>}\hfill \\ r+\frac{k}{n}\left(t-r\right)\mathrm{,}\hfill & \text{if}r+\frac{k}{n}\left(t-r\right)\le v\mathrm{<mo><</mo>}r+\frac{k+1}{n}\left(t-r\right)\mathrm{.}\hfill \end{array}$

Let ${ℱ}_s=\sigma \left\{S_H\left(s\right)\mathrm{,}{I}_H\left(s\right)\mathrm{,}{R}_H\left(s\right)\mathrm{,}{S}_m\left(s\right)\mathrm{,}{I}_m\left(s\right)\mathrm{,0}\le v\le s\right\}$ , and

$A_n\left(s\right)={\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^s}{S}_H\left({\left[v\right]}_n\right)I_m\left({\left[v\right]}_n\right)\text{d}v\mathrm{.}$ (29)

For all $a\mathrm{,}b\in ℝ$ if $a<b$ , then let $P_2\left(\left(a\mathrm{,}b\right]\right)=P_2\left(b\right)-P_2\left(a\right)$ , and $v_k=r+\frac{k}{n}\left(t-r\right)$ . We have

$\begin{array}{l}{P}_2\left(A_n\left(t\right)\right)-P_2\left(A_n\left(r\right)\right)=P_2\left(A_n\left(r\right),A_n\left(t\right)\right]\\ =P_2\left({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^r}{S}_H\left({\left[v\right]}_n\right)I_m\left({\left[v\right]}_n\right)\text{d}v,{\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left({\left[v\right]}_n\right)I_m\left({\left[v\right]}_n\right)\text{d}v\right]\\ =P_2({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^r}{S}_H\left({\left[v\right]}_n\right)I_m\left({\left[v\right]}_n\right)\text{d}v,{\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^r}{S}_H\left({\left[v\right]}_n\right)I_m\left({\left[v\right]}_n\right)\text{d}v+{\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_r^t}{S}_H\left({\left[v\right]}_n\right)I_m\left({\left[v\right]}_n\right)\text{d}v]\end{array}$

By applying the Riemann sum, we get

$\begin{array}{l}{P}_2\left(A_n(t)\right)-P_2\left(A_n(r)\right)=P_2\left(A_n\left(r\right),A_n\left(t\right)\right]\\ ={\displaystyle \sum _{k=1}^n}{P}_2\left(A_n\left(v_{k-1}\right),A_n\left(v_{k-1}\right)+{\lambda }_H{\mathcal{N}}_H\frac{t-r}{n}{S}_H\left(v_{k-1}\right)I_m\left(v_{k-1}\right)\right]\\ ={\displaystyle \sum _{k=1}^n}{P}_2\left(A_n\left(v_{k-1}\right),A_n\left(v_k\right)\right].\end{array}$

Hence

$\begin{array}{l}\mathbb{E}\left[P_2\left(A_n\left(t\right)\right)-P_2\left(A_n\left(r\right)\right)|{ℱ}_r\right]\\ =\mathbb{E}\left[{\displaystyle \sum _{k=1}^n}{P}_2\left(\left(A_n\left(v_{k-1}\right),A_n\left(v_k\right)\right]\right)|{ℱ}_r\right]\end{array}$

$\begin{array}{l}=\mathbb{E}\left[{\displaystyle \sum _{k=1}^n}\mathbb{E}\left(P_2\left(\left(A_n\left(v_{k-1}\right),A_n\left(v_{k-1}\right)+{\lambda }_H{\mathcal{N}}_H{S}_H\left(v_{k-1}\right)I_m\left(v_{k-1}\right)\frac{t-r}{n}\right]\right)|{ℱ}_{v_{k-1}}^n\right)|{ℱ}_r\right]\\ =\mathbb{E}\left[{\displaystyle \sum _{k=1}^n}\left\{{\lambda }_H{\mathcal{N}}_H{S}_H\left(v_{k-1}\right)I_m\left(v_{k-1}\right)\frac{t-r}{n}\right\}|{ℱ}_r\right]\\ =\mathbb{E}\left[A_n\left(t\right)-A_n\left(r\right)|{ℱ}_r\right].\end{array}$

Therefore,

$\mathbb{E}\left(P_2\left(A_n\left(t\right)\right)-P_2\left(A_n\left(r\right)\right)|{ℱ}_r\right)=\mathbb{E}\left(A\left(t\right)-A\left(r\right)|{ℱ}_r\right)$ (30)

As n tends to infinity, we obtain (28). The martingale property implies that $\mathbb{E}{ℳ}_2\left(t\right)=0$ . Now we are going to focus on the expression of $\mathbb{E}\left({\left|{ℳ}_2\left(t\right)\right|}^2\right)$ .

For all $t>0$ , for all natural number $n\ge t$ , and for all $0\le i\le n$ , we define $t_i=\frac{it}{n}$ , $\forall i=\overline{\mathrm{1,}n}$ we have

$0=t_0<t_1<\cdots <t_n=t$

a discretion of the interval $\left[\mathrm{0<mo>;</mo>}t\right]$ .

$\begin{array}{l}{\left|{ℳ}_2\left(t\right)\right|}^2={\displaystyle \sum _{i=1}^n}{\left|{ℳ}_2\left(t_i\right)-{ℳ}_2\left(t_{i-1}\right)\right|}^2\\ +2{\displaystyle \sum _{1\le i<j\le n}}\left[{ℳ}_2\left(t_i\right)-{ℳ}_2\left(t_{i-1}\right)\right]\left[{ℳ}_2\left(t_j\right)-{ℳ}_2\left(t_{j-1}\right)\right]\\ \mathbb{E}\left({\left|{ℳ}_2\left(t\right)\right|}^2\right)={\displaystyle \sum _{i=1}^n}\mathbb{E}\left({\left|{ℳ}_2\left(t_i\right)-{ℳ}_2\left(t_{i-1}\right)\right|}^2\right)\end{array}$

because

$\mathbb{E}\left({ℳ}_2\left(t_i\right)-{ℳ}_2\left(t_{i-1}\right)\right)=\mathbb{E}\left({ℳ}_2\left(t_j\right)-{ℳ}_2\left(t_{j-1}\right)\right)=0,$ (31)

and

${ℳ}_2\left(t_i\right)-{ℳ}_2\left(t_{i-1}\right)\perp {ℳ}_2\left(t_j\right)-{ℳ}_2\left(t_{j-1}\right)\text{for}i\ne j\mathrm{.}$

$\mathbb{E}\left[{\left|{ℳ}_2\left(t\right)\right|}^2\right]=\mathbb{E}{\displaystyle \sum _{i=0}^n}{\left|{ℳ}_2\left(t_i\right)-{ℳ}_2\left(t_{i-1}\right)\right|}^2.$ (32)

As $n\to \infty$ ,

${\displaystyle \sum _{i=1}^n}{\left|{ℳ}_2\left(t_i\right)-{ℳ}_2\left(t_{i-1}\right)\right|}^2\to {\displaystyle \sum _{0\le r\le t}}{\left|\Delta {ℳ}_2\left(r\right)\right|}^2\text{a}\text{.s}\mathrm{,}$

where $\Delta {ℳ}_2\left(r\right)$ denotes the jump of the process ${ℳ}_2$ at time r. So

$\mathbb{E}\left[{ℳ}_2{\left(t\right)}^2\right]=\mathbb{E}\left({\displaystyle \sum _{0<r\le t}}{\left|\Delta {ℳ}_2\left(r\right)\right|}^2\right).$

Indeed, as soon as $t_i-t_{i-1}\le 1$ , since $0\le S_H\left(r\right)\mathrm{,}{I}_m\left(r\right)\le 1$ ,

$\begin{array}{l}{\left|{ℳ}_2\left(t_i\right)-{ℳ}_2\left(t_{i-1}\right)\right|}^2\le \frac{2}{{\mathcal{N}}_H^2}{\left|P_2\left(t_i\right)-P_2\left(t_{i-1}\right)\right|}^2+{\lambda }_H^2\left({\displaystyle {\int }_{t_{i-1}}^{t_i}}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right)\\ \le \frac{2}{{\mathcal{N}}_H^2}{\left|P_2\left(t_i\right)-P_2\left(t_{i-1}\right)\right|}^2+{\lambda }_H^2\left(t_i-t_{i-1}\right)\\ {\displaystyle \sum _{i=0}^{n-1}}{\left|{ℳ}_2\left(t_i\right)-{ℳ}_2\left(t_{i-1}\right)\right|}^2\le \frac{2}{{\mathcal{N}}_H^2}{\left|P_2\left(T\right)\right|}^2+{\lambda }_H^{2}T<\infty .\end{array}$

Hence ${\left({ℳ}_2\left(t\right)\right)}_{t\ge 0}$ is a bounded martingale in $L^2$ , then ${\left|{ℳ}_N^{a_0}\left(t\right)\right|}^2$ uniformly integrable. But

${\displaystyle \sum _{0<r\le t}}{\left|\Delta {ℳ}_2\left(r\right)\right|}^2=\frac{1}{{\mathcal{N}}_H^2}{P}_2\left({\lambda }_H{\mathcal{N}}_H{\displaystyle {\int }_0^t}{S}_H\left(r\right)I_m\left(r\right)\text{d}r\right).$

This equality is justified by the properties of the martingale ${ℳ}_2\left(r\right)$ .

$\begin{array}{c}\mathbb{E}\left({\left|{ℳ}_2\left(t\right)\right|}^2\right)=\mathbb{E}\left({\displaystyle \sum _{0\le r\le t}}{\left|\Delta {ℳ}_2\left(r\right)\right|}^2\right)\\ =\mathbb{E}\left(\frac{1}{{\mathcal{N}}_H^2}{P}_2\left({\mathcal{N}}_H{\displaystyle {\int }_0^t}{\lambda }_H{S}_H\left(s\right)I_m\left(s\right)\text{d}s\right)\right)\\ =\frac{1}{{\mathcal{N}}_H}\mathbb{E}\left({\displaystyle {\int }_0^t}{\lambda }_H{S}_H\left(s\right)I_m\left(s\right)\text{d}s\right).\end{array}$

This completes the demonstration.

Now we are going to state a proposition that will be used subsequently.

Proposition 4.1. Let $\left\{N\left(t\right)\mathrm{,}t\ge 0\right\}$ be a rate $\lambda$ Poisson process. Then $t^{-1}N\left(t\right)\to \lambda$ a.s. as $t\to \infty$ .

Proof. $\forall n\in \mathbb{N}$ , $n\ge 1$ ,

$\frac{N\left(n\right)}{n}=\frac{1}{n}{\displaystyle \sum _{i=1}^n}\left(N\left(i\right)-N\left(i-1\right)\right)\to \lambda ,n\to +\infty .$

Since, $\forall i=\overline{\mathrm{1,}n}$ , $N\left(i\right)-N\left(i-1\right)$ follows a Poisson distribution with parameter $\lambda$ .

For all $t>0$

$\frac{N\left(t\right)}{t}=\frac{\left[t\right]{\left[t\right]}^{-1}}{t}N\left(\left[t\right]\right)+\frac{1}{t}\left(N\left(t\right)-N\left(\left[t\right]\right)\right).$

$\left|\frac{N\left(t\right)}{t}-\lambda \right|\le \left|\frac{\left[t\right]{\left[t\right]}^{-1}}{t}N\left(\left[t\right]\right)-\lambda \right|+\frac{1}{t}\left(N\left(\left[t\right]+1\right)-N\left(\left[t\right]\right)\right).$

$\frac{1}{t}\left(N\left(\left[t\right]+1\right)-N\left(\left[t\right]\right)\right)=\frac{1}{t}N\left(\left[t\right]+1\right)-\frac{1}{t}N\left(\left[t\right]\right)\to \lambda -\lambda =\mathrm{0,}t\to +\infty \mathrm{.}$

Likewise

$\left|\frac{\left[t\right]{\left[t\right]}^{-1}}{t}N\left(\left[t\right]\right)-\lambda \right|\to \mathrm{0,}t\to +\infty \mathrm{.}$

Hence the result. □

Corollary 4.1. As ${\mathcal{N}}_H\to \infty$ , ${\mathcal{N}}_m\to \infty$ , $\forall T>0$ ,

$\underset{0\le t\le T}{\sup }\left\{{\displaystyle \sum _{k=2}^7}\left|{ℳ}_k\left(t\right)\right|+{\displaystyle \sum _{k=9}^{11}}\left|{ℳ}_k\left(t\right)\right|\right\}\to 0$

in probability.

Proposition 4.2. As ${\mathcal{N}}_H\to \infty$ , ${\mathcal{N}}_m\to \infty$ , $\forall T>0$ ,

$\underset{0\le t\le T}{\sup }\left\{{\displaystyle \sum _{k=2}^7}\left|{ℳ}_k\left(t\right)\right|+{\displaystyle \sum _{k=9}^{11}}\left|{ℳ}_k\left(t\right)\right|\right\}\to 0$ a.s.

Proof. Proof We consider the term ${ℳ}_2$ . Since the proportions $S_H\left(t\right)$ and $I_m\left(t\right)$ take values in the interval $\left[\mathrm{0,1}\right]$ , we have

$r={\mathcal{N}}_H{\displaystyle {\int }_0^t}{\lambda }_H{S}_H\left(s\right)I_m\left(s\right)\text{d}s\le {\displaystyle {\int }_0^t}{\mathcal{N}}_H{\lambda }_H\text{d}s$

then

$\begin{array}{l}\left|{ℳ}_2\left(t\right)\right|=\frac{1}{{\mathcal{N}}_H}\left|M_2\left(r\right)\right|\\ \underset{0\le t\le T}{\sup }\left|{ℳ}_2\left(t\right)\right|\le \frac{1}{N}\underset{0\le r\le {\mathcal{N}}_H{\lambda }_{H}T}{\sup }\left|M_2\left(r\right)\right|.\end{array}$

The Law of Large Numbers for Poisson processes (see below) tell us that for all $t>0$ ,

$\frac{P_2\left(Nt\right)}{N}\to t\text{a}\text{.s}\text{.}\text{as}N\to \infty \mathrm{.}$

We have a sequence of increasing functions that converges to a continuous function. Therefore, according to the second Dini’s theorem, this convergence is uniform over any compact interval in $ℝ$ .

$\forall T>0$ , $\frac{1}{{\mathcal{N}}_H}\underset{0\le r\le {\lambda }_H{\mathcal{N}}_{H}T}{\sup }\left|M_2\left(r\right)\right|\to 0$ a.s.

The other cases are demonstrated similarly. This completes the proof. □

Theorem 4.1. Law of Large Numbers.

If $\left(S_H\left(0\right)\mathrm{,}{I}_H\left(0\right)\mathrm{,}{R}_H\left(0\right)\mathrm{,}{S}_m\left(0\right)\mathrm{,}{I}_m\left(0\right)\right)\to \left(s_{H\mathrm{,0}}\mathrm{,}{i}_{H\mathrm{,0}}\mathrm{,}{r}_{H\mathrm{,0}}\mathrm{,}{s}_{m\mathrm{,0}}\mathrm{,}{i}_{m\mathrm{,0}}\right)$ a.s ${\mathcal{N}}_H\to \infty$ and ${\mathcal{N}}_m\to \infty$ , then

$\begin{array}{l}\underset{0\le t\le T}{\sup }\{\left|S_H\left(t\right)-s_H\left(t\right)\right|+\left|I_H\left(t\right)-i_H\left(t\right)\right|+\left|R_H\left(t\right)-r_H\left(t\right)\right|\\ +\left|S_m\left(t\right)-s_m\left(t\right)\right|+\left|I_m\left(t\right)-i_m\left(t\right)\right|\}\to 0\text{a}\text{.s}\text{.,}\end{array}$

where $\left\{\left(s_H\left(t\right)\mathrm{,}{i}_H\left(t\right)\mathrm{,}{r}_H\left(t\right)\mathrm{,}{s}_m\left(t\right)\mathrm{,}{i}_m\left(t\right)\right)\mathrm{,}t\ge 0\right\}$ is the unique solution of the ODE

$\{\begin{array}{l}\frac{\text{d}{s}_H\left(t\right)}{\text{d}t}={\mu }_H-{\lambda }_H{s}_H\left(t\right)i_m\left(t\right)-{\mu }_H{s}_H\left(t\right),t>0\hfill \\ \frac{\text{d}{i}_H\left(t\right)}{dt}={\lambda }_H{s}_H\left(t\right)i_m\left(t\right)-\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)i_H\left(t\right),t>0\hfill \\ \frac{\text{d}{r}_H\left(t\right)}{\text{d}t}={\gamma }_H{i}_H\left(t\right)-{\mu }_H{r}_H\left(t\right),t>0\hfill \\ \frac{\text{d}{s}_m\left(t\right)}{\text{d}t}={\mu }_m-{\lambda }_m{s}_m\left(t\right)i_H\left(t\right)-{\mu }_m{s}_m\left(t\right),t>0\hfill \\ \frac{\text{d}{i}_m\left(t\right)}{\text{d}t}={\lambda }_m{s}_m\left(t\right)i_H\left(t\right)-{\mu }_m{i}_m\left(t\right),t>0\hfill \end{array}$ (33)

where

$\{\begin{array}{l}{s}_H\left(0\right)=s_{H,0},i_H\left(0\right)=i_{H,0},\hfill \\ s_m\left(0\right)=s_{m,0},i_m\left(0\right)=i_{m,0},\hfill \\ r_H\left(0\right)=r_{H,0}.\hfill \end{array}$

Proof. Define

$X\left(t\right)=\left(\begin{array}{c}{s}_H\left(t\right)\\ i_H\left(t\right)\\ r_H\left(t\right)\\ s_m\left(t\right)\\ i_m\left(t\right)\end{array}\right),X_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)=\left(\begin{array}{c}{S}_H\left(t\right)\\ I_H\left(t\right)\\ R_H\left(t\right)\\ S_m\left(t\right)\\ I_m\left(t\right)\end{array}\right),$ (34)

$Y_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)=\left(\begin{array}{c}{ℳ}_6\left(t\right)+{ℳ}_7\left(t\right)-{ℳ}_2\left(t\right)\\ {ℳ}_2\left(t\right)-{ℳ}_4\left(t\right)-{ℳ}_5\left(t\right)-{ℳ}_6\left(t\right)\\ {ℳ}_4\left(t\right)-{ℳ}_7\left(t\right)\\ {ℳ}_{11}\left(t\right)-{ℳ}_9\left(t\right)\\ {ℳ}_9\left(t\right)-{ℳ}_{11}\left(t\right)\end{array}\right),$

${\overline{X}}_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)=X\left(t\right)-X_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)$ (35)

and finally

$F\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\end{array}\right)=\left(\begin{array}{c}{\mu }_H-{\lambda }_H{x}_1{x}_5-{\mu }_H{x}_1\\ {\lambda }_H{x}_1{x}_5-\left({\mu }_H+{\gamma }_H+{\alpha }_H\right)x_2\\ -{\mu }_H{x}_3+{\gamma }_H{x}_2\\ {\mu }_m-{\lambda }_m{x}_4{x}_2-{\mu }_m{x}_4\\ {\lambda }_m{x}_4{x}_2-{\mu }_m{x}_5\end{array}\right)\mathrm{.}$

For $0\le x_i\le 1$ and $0\le {x^{\prime }}_i\le 1$ , $i=\overline{\mathrm{1,5}}$ we have

$\Vert F\left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\end{array}\right)-F\left(\begin{array}{c}{x^{\prime }}_1\\ {x^{\prime }}_2\\ {x^{\prime }}_3\\ {x^{\prime }}_4\\ {x^{\prime }}_5\end{array}\right)\Vert \le \mathcal{C}\Vert \left(\begin{array}{c}{x}_1\\ x_2\\ x_3\\ x_4\\ x_5\end{array}\right)-\left(\begin{array}{c}{x^{\prime }}_1\\ {x^{\prime }}_2\\ {x^{\prime }}_3\\ {x^{\prime }}_4\\ {x^{\prime }}_5\end{array}\right)\Vert \mathrm{,}$

where $\mathcal{C}$ is a constant depending on ${\lambda }_H$ , ${\mu }_H$ , ${\alpha }_H$ , ${\gamma }_H$ , ${\lambda }_m$ and ${\mu }_m$ . We have

${\overline{X}}_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)={\overline{X}}_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(0\right)+{\displaystyle {\int }_0^t}\left[F\left(X\left(r\right)\right)-F\left(X_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(r\right)\right)\right]\text{d}r+Y_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right).$

From Proposition 1, for all $T>0$ , ${\sup }_{0\le t\le T}\Vert Y_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)\Vert \to 0$ a.s. as ${\mathcal{N}}_H\to \infty$ , ${\mathcal{N}}_m\to \infty$ . Let ${\epsilon }_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)={\sup }_{0\le r\le t}\Vert Y_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(r\right)\Vert$ . We have

$\Vert {\overline{X}}_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)\Vert \le \Vert {\overline{X}}_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(0\right)\Vert +\mathcal{C}{\displaystyle {\int }_0^t}\Vert {\overline{X}}_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(r\right)\Vert \text{d}r+{\epsilon }_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right).$

It then follows from Gronwall’s Lemma that

$\underset{0\le r\le t}{\sup }\Vert {\overline{X}}_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(r\right)\Vert \le \left(\Vert {\overline{X}}_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(0\right)\Vert +{\epsilon }_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)\right)\exp \left(\mathcal{C}t\right).$

The result then follows from the assumption $\Vert {\overline{X}}_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(0\right)\Vert \to 0$ , plus the fact that ${\epsilon }_{{\mathcal{N}}_H}^{{\mathcal{N}}_m}\left(t\right)\to 0$ a.s. ${\mathcal{N}}_H\to \infty$ , ${\mathcal{N}}_m\to \infty$ . □

5. Numerical Simulations of the Stochastic Model

5.1. Numeric Simulation of CTMC Model

In this section, we examine the disease dynamics using the stochastic model, utilizing the parameter values listed in Table 1. The multi-type branching process assumes that the susceptible populations are sufficiently large and are at disease-free equilibrium. Thus, the initial conditions for susceptible human hosts and mosquitoes are as follows: $S_H\left(0\right)=0.9$ and $S_m\left(0\right)=0.9$ respectively and the initial conditions for the infectious are respectively $I_H\left(0\right)=0.1$ and $I_m\left(0\right)=0.1$ .

In Figure 1, the basic reproduction number is ${ℛ}_0=1.1447$ , and the probability of a major outbreak is given by $P_{ps}=1-{\mathbb{P}}_{ext}=0.1315$ .

Figure 1. CTMC trajectories for dengue extinction.

Remark 5.1. The stochastic model allows for determining not only the extinction of a disease or the emergence of an epidemic but also the probability of these events. This is achieved by applying the theory of multitype branching processes, particularly when the epidemic starts with a small number of infected individuals, a scenario that deterministic models cannot address (Allen and van den Driessche 2013). The stochastic model suggests that it’s possible for the disease to die out, as depicted in Figure 1, while the deterministic model guarantees that the dengue epidemic will occur.

5.2. Highlighting the Law of Large Numbers

In this section, we highlight the law of large numbers by varying the population size in an increasing manner and observing the behaviour of each curve for the two models (deterministic model and stochastic model).

Figure 2 and Figure 3 demonstrate the Law of Large Numbers for ${\mathcal{N}}_H={10}^4$ , ${\mathcal{N}}_m={10}^4$ , ${\mathcal{N}}_H={10}^5$ , and ${\mathcal{N}}_m={10}^5$ respectively.

Figure 2. Illustrating the law of large numbers for ${\mathcal{N}}_H={10}^4$ and ${\mathcal{N}}_m={10}^4$ .

Figure 3. Demonstration of the law of large numbers for ${\mathcal{N}}_H={10}^5$ and ${\mathcal{N}}_m={10}^5$ .

Remark 5.2. As ${\mathcal{N}}_H$ and ${\mathcal{N}}_m$ become larger and larger, the amplitude of variations in the curves of the stochastic model depicted in Figure 2 and Figure 3 becomes almost negligible. In other words, the stochastic model converges towards the deterministic model. Thus, the Law of Large Numbers is satisfied.

6. Conclusion

In this article, we formulated and analysed a continuous-time Markov chain (CTMC) model, exploring the asymptotic behaviour of the stochastic model for large populations using the Law of Large Numbers (LLN). Our main finding is that disease extinction is possible when ${ℛ}_0$ exceeds one, and we determined the population threshold beyond which a deterministic model becomes suitable. In the context of a dengue epidemic, the corresponding stochastic model is applicable for small population sizes. However, for very large populations, transitioning to the deterministic model simplifies the study complexities. Nevertheless, several intriguing topics remain unexplored. As dengue fever spreads, individuals accumulate infection-related knowledge. Therefore, our future research aims to investigate the impact of memory on our model dynamics using new generalized and fractal fractional derivatives, as well as a novel numerical method for solving EDFs with the GHF derivative, introduced by Khalid Hattaf et al. in [36] and [37].

Acknowledgements

Sincere thanks to the members of JAMP for their professional performance, and special thanks to managing editor Hellen XU for a rare attitude of high quality.

Conflicts of Interest

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

References

[1]	Racloz, V., Ramsey, R., Tong, S. and Hu, W. (2012) Surveillance of Dengue Fever Virus: A Review of Epidemiological Models and Early Warning Systems. PLOS Neglected Tropical Diseases, 6, e1648. https://doi.org/10.1371/journal.pntd.0001648
[2]	Gibbons, R.V. (2002) Dengue: An Escalating Problem. British Medical Journal, 324, 1563-1566. https://doi.org/10.1136/bmj.324.7353.1563
[3]	Gubler, D.J. (2002) Epidemic Dengue/Dengue Hemorrhagic Fever as a Public Health, Social and Economic Problem in the 21st Century. Trends in Microbiology, 10, 100-103. https://doi.org/10.1016/s0966-842x(01)02288-0
[4]	WHO, Alert, P., et al. (2009) Dengue: Guidelines for Diagnosis, Treatment, Prevention and Control.
[5]	Kiouach, D., El-idrissi, S.E.A. and Sabbar, Y. (2023) A Novel Mathematical Analysis and Threshold Reinforcement of a Stochastic Dengue Epidemic Model with Lévy Jumps. Communications in Nonlinear Science and Numerical Simulation, 119, Article 107092. https://doi.org/10.1016/j.cnsns.2023.107092
[6]	Bhatt, S., Gething, P.W., Brady, O.J., Messina, J.P., Farlow, A.W., Moyes, C.L., et al. (2013) The Global Distribution and Burden of Dengue. Nature, 496, 504-507. https://doi.org/10.1038/nature12060
[7]	Lee, J., Mogasale, V., Lim, J.K., Ly, S., Lee, K.S., Sorn, S., et al. (2019) A Multi-Country Study of the Economic Burden of Dengue Fever Based on Patient-Specific Field Surveys in Burkina Faso, Kenya, and Cambodia. PLOS Neglected Tropical Diseases, 13, e0007164. https://doi.org/10.1371/journal.pntd.0007164
[8]	Wichmann, O., Yoon, I., Vong, S., Limkittikul, K., Gibbons, R.V., Mammen, M.P., et al. (2011) Dengue in Thailand and Cambodia: An Assessment of the Degree of Underrecognized Disease Burden Based on Reported Cases. PLOS Neglected Tropical Diseases, 5, e996. https://doi.org/10.1371/journal.pntd.0000996
[9]	Cai, L., Guo, S., Li, X. and Ghosh, M. (2009) Global Dynamics of a Dengue Epidemic Mathematical Model. Chaos, Solitons & Fractals, 42, 2297-2304. https://doi.org/10.1016/j.chaos.2009.03.130
[10]	Monath, T.P. (1994) Dengue: The Risk to Developed and Developing Countries. Proceedings of the National Academy of Sciences, 91, 2395-2400. https://doi.org/10.1073/pnas.91.7.2395
[11]	Stanaway, J.D., Shepard, D.S., Undurraga, E.A., Halasa, Y.A., Coffeng, L.E., Brady, O.J., et al. (2016) The Global Burden of Dengue: An Analysis from the Global Burden of Disease Study 2013. The Lancet Infectious Diseases, 16, 712-723. https://doi.org/10.1016/s1473-3099(16)00026-8
[12]	Swaminathan, S. and Khanna, N. (2019) Dengue Vaccine Development: Global and Indian Scenarios. International Journal of Infectious Diseases, 84, S80-S86. https://doi.org/10.1016/j.ijid.2019.01.029
[13]	Khan, M.A. and Fatmawati (2021) Dengue Infection Modeling and Its Optimal Control Analysis in East Java, Indonesia. Heliyon, 7, e06023. https://doi.org/10.1016/j.heliyon.2021.e06023
[14]	Agusto, F.B. and Khan, M.A. (2018) Optimal Control Strategies for Dengue Transmission in Pakistan. Mathematical Biosciences, 305, 102-121. https://doi.org/10.1016/j.mbs.2018.09.007
[15]	Focks, D.A., Haile, D.G., Daniels, E. and Mount, G.A. (1993) Dynamic Life Table Model for Aedes Aegypti (Diptera: Culicidae): Analysis of the Literature and Model Development. Journal of Medical Entomology, 30, 1003-1017. https://doi.org/10.1093/jmedent/30.6.1003
[16]	Wang, L. and Zhao, H. (2019) Dynamics Analysis of a Zika-Dengue Co-Infection Model with Dengue Vaccine and Antibody-Dependent Enhancement. Physica A: Statistical Mechanics and its Applications, 522, 248-273. https://doi.org/10.1016/j.physa.2019.01.099
[17]	Champagne, C. and Cazelles, B. (2019) Comparison of Stochastic and Deterministic Frameworks in Dengue Modelling. Mathematical Biosciences, 310, 1-12. https://doi.org/10.1016/j.mbs.2019.01.010
[18]	Din, A., Khan, T., Li, Y., Tahir, H., Khan, A. and Khan, W.A. (2021) Mathematical Analysis of Dengue Stochastic Epidemic Model. Results in Physics, 20, Article 103719. https://doi.org/10.1016/j.rinp.2020.103719
[19]	Guo, M., Hu, L. and Nie, L. (2021) Stochastic Dynamics of the Transmission of Dengue Fever Virus between Mosquitoes and Humans. International Journal of Biomathematics, 14, Article 2150062. https://doi.org/10.1142/s1793524521500625
[20]	Liu, Q. (2023) Dynamical Behavior of a Stochastic Dengue Model with Ornstein–Uhlenbeck Process. Journal of Mathematical Physics, 64, Article 092705. https://doi.org/10.1063/5.0147354
[21]	Shatanawi, W., Raza, A., Shoaib Arif, M., Rafiq, M., Bibi, M. and Mohsin, M. (2021) Essential Features Preserving Dynamics of Stochastic Dengue Model. Computer Modeling in Engineering & Sciences, 126, 201-215. https://doi.org/10.32604/cmes.2021.012111
[22]	Kiouach, D., Sabbar, Y. and El Azami El-Idrissi, S. (2021) New Results on the Asymptotic Behavior of an SIS Epidemiological Model with Quarantine Strategy, Stochastic Transmission, and Lévy Disturbance. Mathematical Methods in the Applied Sciences, 44, 13468-13492. https://doi.org/10.1002/mma.7638
[23]	Zhou, Y. and Zhang, W. (2016) Threshold of a Stochastic SIR Epidemic Model with Lévy Jumps. Physica A: Statistical Mechanics and its Applications, 446, 204-216. https://doi.org/10.1016/j.physa.2015.11.023
[24]	Konané, F.V. and Sawadogo, R. (2023) Mathematical Analysis of a Deterministic and a Stochastic Epidemic Models of Dengue. European Journal of Applied Mathematics, 3, 1-28.
[25]	Allen, L.J.S. and van den Driessche, P. (2013) Relations between Deterministic and Stochastic Thresholds for Disease Extinction in Continuous and Discrete-Time Infectious Disease Models. Mathematical Biosciences, 243, 99-108. https://doi.org/10.1016/j.mbs.2013.02.006
[26]	Harris, T.E., et al. (1963) The Theory of Branching Processes. Springer.
[27]	Hastings, A. and Gross, L.J. (2012) Encyclopedia of Theoretical Ecology. University of California Press.
[28]	Allen, L.J. (2015) Stochastic Population and Epidemic Models. Mathematical Biosciences Lecture Series, Stochastics in Biological Systems, 1, 120-128.
[29]	Maliyoni, M., Chirove, F., Gaff, H.D. and Govinder, K.S. (2017) A Stochastic Tick-Borne Disease Model: Exploring the Probability of Pathogen Persistence. Bulletin of Mathematical Biology, 79, 1999-2021. https://doi.org/10.1007/s11538-017-0317-y
[30]	Nandi, A. and Allen, L.J.S. (2018) Stochastic Two-Group Models with Transmission Dependent on Host Infectivity or Susceptibility. Journal of Biological Dynamics, 13, 201-224. https://doi.org/10.1080/17513758.2018.1538462
[31]	Lahodny, G.E. and Allen, L.J.S. (2013) Probability of a Disease Outbreak in Stochastic Multipatch Epidemic Models. Bulletin of Mathematical Biology, 75, 1157-1180. https://doi.org/10.1007/s11538-013-9848-z
[32]	Lloyd, A.L., Zhang, J. and Root, A.M. (2007) Stochasticity and Heterogeneity in Host–vector Models. Journal of The Royal Society Interface, 4, 851-863. https://doi.org/10.1098/rsif.2007.1064
[33]	Ball, F. (1983) The Threshold Behaviour of Epidemic Models. Journal of Applied Probability, 20, 227-241. https://doi.org/10.2307/3213797
[34]	Bartlett, M.S. (1964) The Relevance of Stochastic Models for Large-Scale Epidemiological Phenomena. Applied Statistics, 13, 2-8. https://doi.org/10.2307/2985217
[35]	Pardoux, É. (2015) Mathematical Models of Epidemics Cimpa School Ziguinchor.
[36]	Hattaf, K. (2022) On the Stability and Numerical Scheme of Fractional Differential Equations with Application to Biology. Computation, 10, Article 97. https://doi.org/10.3390/computation10060097
[37]	Hattaf, K. (2023) A New Class of Generalized Fractal and Fractal-Fractional Derivatives with Non-Singular Kernels. Fractal and Fractional, 7, Article 395. https://doi.org/10.3390/fractalfract7050395.