Modeling mosquito-borne disease dynamics via stochastic differential equations and generalized tempered stable distribution

1.   Introduction

Historically, animal pathogens have wielded significant global impact. To illustrate this, in the 14th Century, the Black Death, a disease transmitted from rodents to humans through infected fleas, led to the deaths of approximately fifty million individuals worldwide, cementing its status as one of the deadliest diseases in human history ^[1]. Subsequently, in the 20th century, diseases like monkeypox, rabies, and avian flu emerged, straining healthcare systems in numerous countries across the globe ^[2]. Zoonotic diseases, unlike those transmitted among humans, initially spread from animals to humans, presenting a substantial challenge due to the diverse nature of the affected populations ^[3]. Consequently, comprehending the potential pathways of disease transmission between different sub-populations necessitates a comprehensive examination of interactions between animals and humans, resulting in a higher level of complexity in the analysis of their dynamics ^[4,5]. Among all animals, mosquitoes stand out as the most formidable vectors in the rapid propagation of major zoonotic diseases, including zika virus, west Nile virus, chikungunya virus, dengue fever, yellow fever, and malaria. These diseases, known as mosquito-borne diseases, infect millions of people worldwide and contribute to a significant number of fatalities annually ^[6,7].

To understand the complexity of diseases transmitted by mosquitoes and identify the elements that favor their spread within populations, mathematicians have proposed numerous models. These models find their origins mainly in the fundamental contributions of Kermack and McKendrick ^[8,9] to the mathematical theory of diseses. The basic concept is to subdivide the study population into distinct groups, based on potential clinical states induced by the infection. In its simplest form, this model is represented by a system composed of interconnected nonlinear ordinary differential equations that explain how human and mosquito populations change over time, taking into account all possible interactions. In this framework, various models have been used, giving revolutionary results, as evidenced by references such as ^{[10,11,12,13,14,15,16,17]}. Nevertheless, as Mao demonstrated in ^[18], stochastic environmental fluctuations have the potential to alter the dynamics of epidemic and ecological patterns. Therefore, an extension of deterministic results to the stochastic domain becomes imperative.

Mosquito-borne diseases, such as malaria or dengue fever, are influenced by various random factors in the environment, including mosquito behavior, climate conditions, and human movement ^[19,20]. Ignoring these fluctuations can result in overly simplistic models that do not accurately reflect the complex dynamics of disease transmission. Incorporating stochastic elements into models allows for a more realistic prediction of disease outbreaks. By accounting for random fluctuations, models can provide a range of possible outcomes, including worst-case scenarios. This is particularly important for public health planning and resource allocation. To highlight the main innovations of our research, we begin by reviewing the existing body of literature regarding models of mosquito-borne epidemics. For example, Witbooi et al. introduced a stochastic mosquito-transmitted model in ^[21], incorporating Gaussian perturbations in disease transmission rates. This model characterizes the interaction between susceptible, infected, and recovered human individuals, as well as susceptible and infected mosquitoes. The authors not only demonstrated the existence of a unique and almost certainly positive solution to their model, but also outlined the conditions under which the disease is effectively eradicated within the population studied. In the context of modeling malaria dynamics, Wang et al. proposed a stochastic model in ^[22]. Similar to the Witbooi model, it divides human and mosquito populations accordingly. However, instead of Gaussian disturbances, this model uses multiplicative Gaussian noise. Furthermore, the authors extended previous results by establishing conditions guaranteeing the presence of a stationary distribution. For the spread of dengue, Liu et al. adopted the same model in ^[23], but with a modified saturated incidence rate to prevent the bilinear contact rate from being unlimited. Their work builds on the results of ^[22]. Expanding the scope to a broader stochastic model for dengue, Wei et al. introduced a model that incorporates the human population exposed in ^[24]. In the context of Chikungunya disease modeling, Gokila and Sambath, the authors of ^[25] established a well-structured stochastic model. They derived the conditions for disease extinction and the presence of a stationary distribution.

In both stochastic and deterministic scenarios, the mosquito-borne epidemic models mentioned above display three significant constraints, as delineated below:

● Restriction 1: Given that vaccines for most mosquito-borne diseases, such as Zika virus and Chikungunya virus, are not yet proven to be effective and remain in the preclinical trial phase, coupled with the absence of treatment options ^[26], governments are left with only one viable control measure for safeguarding the human population: implementing quarantine measures for the infected individuals. It is worth noting that previous models have typically overlooked the inclusion of quarantine strategies. However, it is essential to highlight that various studies have demonstrated the beneficial impact of quarantine as a control measure for different mosquito-borne diseases. Notable examples include research presented in ^[27,28,29]. Therefore, neglecting to incorporate quarantine measures in the model formulation could diminish the model's relevance.

● Restriction 2: Gaussian noise is appropriate for modeling minor fluctuations, but it is not well-suited for simulating significant and abrupt changes, such as those caused by natural events like volcanoes and earthquakes, impacting populations ^[30,31]. These types of natural disasters can play a pivotal role in the transmission of mosquito-borne diseases by forcing people to relocate and gather in unsanitary conditions with inadequate sanitation facilities, potentially leading to the resurgence of mosquito-borne diseases ^[32]. Given the sudden and discontinuous nature of these disaster events, the Gaussian noise employed in the cited literature is inadequate for capturing this effect ^[33].

Mosquito-borne diseases like malaria, dengue, Zika, and Chikungunya pose persistent global health challenges due to their complex transmission dynamics influenced by vector behavior, environmental factors, and human interactions. Accurate modeling of these diseases is crucial for understanding their epidemiology and developing effective control strategies. This study introduces an advanced stochastic model that integrates quarantine measures and employs Lévy noise ^[34,35,36]. Unlike traditional deterministic models, our approach captures the inherent variability and unpredictability of disease transmission, crucial for predicting sporadic outbreaks and assessing the impact of interventions. By incorporating stochastic elements, the model enhances realism in depicting disease spread patterns, facilitates scenario analysis for outbreak prediction under various conditions, and informs policy development to mitigate mosquito-borne disease burden globally.

● Restriction 3: Taking into account epidemic models with jumps associated with the standard measure shows certain restrictions in the modeling of some phenomena with heavy tails (see ^[37,38]). So, incorporating a generalized tempered stable (GTS) distribution into an epidemic model with jumps can offer a more realistic representation of the stochastic nature of disease transmission. The GTS distribution is a flexible probability distribution that allows for modeling various degrees of tail heaviness, which can be particularly valuable when modeling epidemic dynamics with occasional extreme events or jumps ^[39]. Mosquito-borne models often face rare but significant events that can lead to sudden and substantial changes in disease dynamics. These events may include super-spreader gatherings, sudden policy interventions, or unexpected changes in population behavior. The GTS distribution can be used to model the distribution of the sizes of these jumps, allowing for more accurate simulations of their impact on the epidemic ^[40]. Furthermore, the GTS distribution can capture the heterogeneity in disease transmission rates among individuals. Some individuals may play a more substantial role in spreading the disease than others. The GTS distribution can be used to model the variability in transmission rates, which may follow a power-law distribution ^[41].

Considering the preceding discussion, the primary innovations in this paper involve the extension of previously established findings to a novel, well-defined stochastic model that is both mathematically rigorous and biologically plausible. This model addresses the constraints denoted as Restrictions 1–3 across a broad spectrum of mosquito-borne diseases. To the best of our knowledge, this marks the inaugural attempt to present an epidemic model that tackles both of these limitations simultaneously.

In technical terms, our newly proposed stochastic system lacks both endemic and disease-free states. Consequently, the traditional approach of examining the disease's persistence or extinction by studying its asymptotic behavior around these states is not applicable. Therefore, we must adopt a novel approach rooted in stochastic analysis. It's worth noting that in the deterministic case, most information regarding disease eradication and continuation can be similarly obtained by setting the stochastic intensities to zero.

The subsequent sections of this article are organized in accordance with the following structure: In Section 2, we present our improved deterministic and stochastic models by illustrating their components and parameters. In Section 3, we present findings related to the mathematical soundness, biological viability, and the population's long-term dynamics in the absence of disease transmission. In Section 4, we establish the conditions that lead to the extinction of both infected human and mosquito populations. Subsequently, we determine the criteria for the infection to persist within both human and mosquito populations. In Section 5, we perform numerical simulations to study the sensitivity analysis and conduct experiments to validate the theoretical results. Finally, in Section 6, we deliberate on conclusions and offer additional insights.

2.   Deterministic and stochastic models

Considering the deterministic mosquito-borne epidemic models mentioned earlier and with the goal of addressing Restrictions 1–3, our initial approach involves dividing the entire human population into four distinct categories. This includes the introduction of a new category to represent individuals under quarantine. Specifically, at a given time $t\geq0$, we define the following:

● ${S}_{\mathfrak{h}}(t)$: Density of the susceptible human population.

● ${I}_{\mathfrak{h}}(t)$: Density of the infected human population.

● ${Q}_{\mathfrak{h}}(t)$: Density of the quarantined human population.

● ${R}_{\mathfrak{h}}(t)$: Density of the recovered human population.

For the mosquito population, we consider the following classes:

● ${S}_{\mathfrak{m}}(t)$: Density of the susceptible mosquito population.

● ${I}_{\mathfrak{m}}(t)$: Density of the infected mosquito population.

In the absence of stochastic noise, our proposed model is governed by the subsequent system of interconnected ordinary differential equations:

provided with the following initial value problems:

All parameters introduced in model (2.1) are supposed to be taking positive values, and their biological significations are listed as follows:

● $a$: Standard rate of recruitment associated with susceptible persons.

● ${a}_m$: Constant influx of new susceptible mosquitoes.

● ${b}$: Dissemination rate of the infection from mosquitoes to humans.

● ${b}_m$: Dissemination rate of the infection from humans to mosquitoes.

● ${c}$: Death rate of the human population.

● ${c}_0$: Disease-caused death rate of the human population.

● ${c}_m$: Natural death rate of the mosquito population.

● ${f_1}$: Rate in which the human individuals leave quarantine.

● ${f_2}$: Rate in which the infected human individuals put into quarantine recover from the disease.

● ${g_1}$: Death rate of the infected human population in quarantine.

● ${g_2}$: Disease-caused death rate of the infected human population.

● ${q_1}$: Natural recovery rate.

● ${q_2}$: Rate in which infected individuals receive medication (when available).

● $\varphi$: Rate in which the infected human population is put into quarantine.

In the realm of deterministic epidemic models, a critical parameter that holds sway over the dynamics is known as the basic reproduction number ($\mathcal{R}_{\circ}$) and is calculated using the next generation method ^[42]. For system (2.1), it holds that $\mathcal{R}_{\circ} = \dfrac{b {a} {b}_m {a}_m}{{c}{c}_m^2\left({c}+{c}_0+\varphi+q_1+q_2\right)}.$ In the deterministic framework, it is possible to delineate the endurance and disappearance of the disease by defining the circumstances under which the uninfected and endemic equilibrium points exhibit local and/or global asymptotic stability. Nevertheless, as we will illustrate, the presence of stochastic noise renders these findings no longer applicable. Indeed, the criterion ensuring the persistence or eradication of the disease in a deterministic epidemic model may not be applicable to its stochastic counterparts, influenced by either Gaussian or Lévy noise. Therefore, to establish a more resilient foundation for our analysis, we adopt a stochastic modeling approach, which enables us to address Restrictions 2 and 3. To achieve this, we introduce a stochastic perturbation in the form of Lévy noise to each equation of the deterministic model (2.1). In this context, we consider independent jump processes for each sub-population. Practically speaking, this method better mirrors real-world situations, given that the discrete factors affecting the dynamics may vary among different sub-populations. Hence, the model at hand is described by the following interconnected system of stochastic differential equations:

Here, $\mathcal{A}_{L}\; ({L} = 1, \cdots, 6)$, denote six mutually independent Brownian motions (BMs) of strengths ${a}_{L} > 0\; ({L} = 1, \cdots, 6)$, respectively. All these BMs are essentially defined on a filtered probability triple (stochastic basis) $\left(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t\geqslant 0}, \mathbb{P}\right)$ equipped with a filtration satisfying the usual criteria. ${S}_{\mathfrak{h}}(t^-), {I}_{\mathfrak{h}}(t^-), {Q}_{\mathfrak{h}}(t^-), {R}_{\mathfrak{h}}(t^-), {S}_{\mathfrak{m}}(t^-)$ and ${I}_{\mathfrak{m}}(t^-)$ are denoting the left limits of ${S}_{\mathfrak{h}}(t), {I}_{\mathfrak{h}}(t), {Q}_{\mathfrak{h}}(t), {R}_{\mathfrak{h}}(t), {S}_{\mathfrak{m}}(t)$ and ${I}_{\mathfrak{m}}(t)$. $\mathcal{Z}_{L}\; ({L} = 1, \cdots, 6)$ are six independent Poisson counting associated with six finite characteristic Lévy measures ${Q}_{L}\; ({L} = 1, \cdots, 6)$ defined on a measurable set ${\mathbb R}^6\setminus\{0\}$ as follows:

where ${1}$ denotes the indicator function, $\alpha\in(0, 2)$, and $R_L$ is the Rosiński measure defined on ${\mathbb R}^6\setminus\{0\}$ such that $R_L(0) = 0$,

$\widetilde{\mathcal{Z}}_{L}\; ({L} = 1, \cdots, 6)$ are six different compensated random measures such that

Last, ${z}_{L}: {\mathbb R}^6\setminus\{0\}\to\mathbb{R}$ are the jumps size functions which are postulated to be continuous on ${\mathbb R}^6\setminus\{0\}$.

In this paper, we consider a generalized tempered stable (GTS) distribution by taking a new Lévy measure ${Q}_{0}$ defined as follows:

Let $\alpha_{-}\in(0, 2)$, $\alpha_{+}\in(0, 2)$, $\beta_{-} > 0$, $\beta_{+} > 0$, $\sigma_{-} > 0$, and $\sigma_{+} > 0$. Then

The tempered stable distribution associated with the measure defined in (2.4) is a general framework of some well-known special cases presented in the literature:

● By picking out $\alpha_{+} = \alpha_{-} = 0$ and $\beta_{+} = \beta_{-}$, we get the variance gamma distribution exhibited in ^[38].

● By selecting $\alpha_{+} = \alpha_{-}$, we get the KoBoL distribution discussed in ^[39].

● By picking $\alpha_{+} = \alpha_{-}$ and $\sigma_{+} = \sigma_{-}$, we get the infinitely divisible distribution linked to a truncated Lévy flight introduced in ^[40].

● By choosing $\alpha_{+} = \alpha_{-} = 0$, we get the bilateral gamma distribution explained in ^[41].

● By choosing $\alpha_{+} = \alpha_{-}$ and $\beta_{+} = \beta_{-}$, we get the CGMY-distribution presented in ^[43].

3.   Mathematical well-posedness, biological feasibility, and long-time behavior in the absence of the disease transmission

The initial stage of examining model (2.2) involves confirming its well-posedness. This entails demonstrating that, for any initial conditions and constant parameters, there exists a sole solution that delineates the population's evolution over all considered time spans. Furthermore, as we are dealing with the modeling of a biological phenomenon, it is imperative to ensure that the resulting solution is biologically plausible, meaning that it maintains non-negativity throughout its trajectory.

To streamline our mathematical calculations and ensure brevity, we will consistently employ the following notations and definitions throughout the paper.

● ${ } \chi_1: = \underset{L \in (1, \cdots, 6)}{\max}\left(\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}{z}_{L}^2(\xi){Q}_L({\rm{d}}\xi)\right)$.

● ${ } \chi_2: = \underset{L \in (1, \cdots, 6)}{\max} \left(\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}\Big({z}_{L}(\xi)-\ln\left(1+{z}_{L}(\xi)\right)\Big){Q}_L({\rm{d}}\xi)\right)$.

● ${ } \chi_3: = \underset{L \in (1, \cdots, 6)}{\max}\left(\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}\Big(\Big(1+{z}_{L}(\xi)\Big)^2-1\Big)^2{Q}_L({\rm{d}}\xi)\right)$.

● ${ } \chi_4: = \underset{L \in (1, \cdots, 6)}{\max}\left(\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}\Big(\ln\Big(1+{z}_{L}(\xi)\Big)\Big)^2{Q}_L({\rm{d}}\xi)\right)$.

● ${ } \chi_5: = \underset{{L}\in(1, \cdots, 6)}{\max}\big\lbrace {a}_{L}^2 \big\rbrace$.

● ${ } \chi_6(\xi): = \underset{{L}\in(1, \cdots, 6)}{\max}\big\lbrace {z}_{L}(\xi) \big\rbrace = {z}_{L^*}(\xi)$, where $L^*$ denotes the index for which the maximum is attained.

● ${ } \chi_7(\xi): = \underset{{L}\in(1, \cdots, 6)}{\min}\big\lbrace {z}_{L}(\xi) \big\rbrace = {z}_{\overline{L}}(\xi)$, where $\overline{L}$ denotes the index for which the minimum is attained.

● ${ } \chi_8(\xi): = \left(1+\chi_6(\xi)\right)^v-v\times\chi_6(\xi)-1$.

● ${ } \chi_9(\xi): = \left(1+\chi_7(\xi)\right)^v-v\times\chi_7(\xi)-1$.

● ${ } \chi_{10}(\xi): = \max\left\lbrace\chi_8(\xi), \chi_9(\xi)\right\rbrace$.

● ${ } \chi_{11}: = {\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}} \chi_{10}(\xi) {1}_{\left(\chi_{8}(\xi)\geq \chi_{9}(\xi)\right)} {Q}_{L^*}({\rm{d}}\xi)+ {\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}} \chi_{10}(\xi) {1}_{\left(\chi_{9}(\xi) > \chi_{8}(\xi)\right)} {Q}_{\overline{L}}({\rm{d}}\xi).$

Moreover, to establish both the mathematical and biological coherence of the proposed model, we introduce the following technical assumptions

● Presumption A: ${z}_{L}(\xi) > -1, \quad \forall L \in (1, \cdots, 6)$.

● Presumption B: $\underset{{{L}\in (1, \cdots, 6)}}{\max} \big(\chi_L\big) < \infty$.

● Presumption C: $\exists v > 2$ such that ${ }\chi_{12} = \min\big({c}, {c}_{\mathfrak{m}}\big)-\frac{(v-1)}{2}\chi_5-\dfrac{1}{v}\chi_{11} > 0$.

Remark 3.1. To underscore the significance of Presumption C, we direct the reader's attention to [45,Lemma 2.5]. In this lemma, the authors provided fundamental results that serve as the foundation for establishing the main outcomes in our study.

Theorem 3.1. Assuming that Presumptions A and B are satisfied, the perturbed model (2.2) exhibits strong mathematical and biological well-posedness, as it possesses a unique solution that is almost sure global and positive.

Proof. The proof follows the Lyapunov approach inspired by Mao's work ^[44], a technique that has been applied in several studies involving diverse stochastic epidemic models. For example, authors can refer to ^{[45,46,47,48,49]}. However, the distinctive aspect in our approach is the specific form of the selected Lyapunov function, which, in our case, is defined as follows:

Since the coefficients of system (2.2) satisfy the local-Lipschitz property, there exists a unique solution $({S}_{\mathfrak{h}}, {I}_{\mathfrak{h}}, {Q}_{\mathfrak{h}}, {R}_{\mathfrak{h}}, {S}_{\mathfrak{m}}, {I}_{\mathfrak{m}}),$ defined on $(0, \tau_e),$ where $\tau_e$ denotes the explosion time.

Now, for $N \in \mathbb{N}^*$ large enough such that ${S}^0_{\mathfrak{h}}, {I}^0_{\mathfrak{h}}, {Q}^0_{\mathfrak{h}}, {R}^0_{\mathfrak{h}}, {S}^0_{\mathfrak{m}}, {I}^0_{\mathfrak{m}} \in \left(\dfrac{1}{N}, N\right),$ we consider the following stopping time

with the usual convention $\inf(\emptyset) = +\infty,$ and such that $\emptyset$ denotes the empty set.

Clearly, the sequence $(\tau_n)_{n \geq 0}$ is increasing and bounded, from whence there exists $\tau_\infty,$ such that

To conclude the proof, is suffices to prove that $\tau_\infty = +\infty.$ To this end, we reason by contradiction and assume the opposite holds. That is, there exist $\epsilon \in (0, 1)$ and $\hat{T} > 0$ such that $\mathbb{P}\left(\tau_\infty \leq \hat{T} \right) \geq \epsilon.$ By utilizing Itô's formula, for $t \in (0, \hat{T} \wedge \tau_n),$ where $a \wedge b : = \min(a, b), \forall a, b \in \mathbb{R},$ we obtain

where

By rearranging the terms and taking the positiveness of $({S}_{\mathfrak{h}}, {I}_{\mathfrak{h}}, {Q}_{\mathfrak{h}}, {R}_{\mathfrak{h}}, {S}_{\mathfrak{m}}, {I}_{\mathfrak{m}})$ in $(0, \hat{T} \wedge \tau_n)$ into account, it is straightforward that

Therefore,

Consequently,

An integration of (3.2) from 0 to $\hat{T} \wedge \tau_n$ and an evaluation of the expectation on both of its sides yield

Define $\mathcal{U}_n: = \left(\tau_n \leq \hat{T} \right).$ Then remark that for $\omega \in \Omega$ one of the quantities: ${S}_{\mathfrak{h}}(\tau_n, \omega),$ ${I}_{\mathfrak{h}}(\tau_n, \omega),$ ${Q}_{\mathfrak{h}}(\tau_n, \omega),$ ${R}_{\mathfrak{h}}(\tau_n, \omega),$ ${S}_{\mathfrak{m}}(\tau_n, \omega),$ ${I}_{\mathfrak{m}}(\tau_n, \omega)$ is equal to either $n$ or $\dfrac{1}{n}.$ With this in mind, we acquire that

Letting $n \rightarrow + \infty$ leads to the contradiction $+\infty > +\infty.$ This concludes the proof. □

We will now present findings concerning the dynamic patterns of the human and mosquito populations when there is no disease transmission occurring. From a biological perspective, this information is pivotal as it sheds light on the potential impact of various hypothetical control measures aimed at significantly reducing or completely halting disease transmission within the studied population.

To present the aforementioned findings, we examine the behavior of our model when disease transmission is absent. We will provide a concise overview of some asymptotic properties pertaining to the boundary equations incorporated into model (2.2). To do so, we employ the following two-dimensional auxiliary system with Lévy noise and GTS distribution:

Remark 3.2. Due to the structure of our model, it is insufficient to exclusively examine the first equation (pertaining to ${S}_{\mathfrak{h}}(t)$) and apply the comparison theorem, as it contains a positive term (${f_1} {Q}_{\mathfrak{h}}(t)$). Therefore, it is imperative to consider the entire human population to ensure that ${K}_{\mathfrak{h}}(t)\leq {S}_{\mathfrak{h}}(t)+{I}_{\mathfrak{h}}(t)+{Q}_{\mathfrak{h}}(t)+{R}_{\mathfrak{h}}(t)$ holds almost surely. In the case of the second equation, we can straightforwardly establish ${K}_{\mathfrak{m}}(t)\leq {S}_{\mathfrak{m}}(t)$ with certainty.

Lemma 3.1. Consider two Markov processes $\left({K}_{\mathfrak{h}}(t), {K}_{\mathfrak{m}}(t)\right)$ that conform to the bi-dimensional auxiliary system (3.3). Given the fulfillment of Presumptions A, B, and C, the following properties are observed:

● ${ }\lim\limits_{t\to\infty}t^{-1}\; \displaystyle {\int}_0^t{K}_{\mathfrak{h}}(s)\: {\rm{d}}s = \dfrac{{a}}{{c}}\; \; {\rm{a.s}}$ and ${ }\lim\limits_{t\to\infty}t^{-1}\; \displaystyle {\int}_0^t {K}_{\mathfrak{m}}(s)\: {\rm{d}}s = \dfrac{{a}_{\mathfrak{m}}}{{c}_{\mathfrak{m}}}\; \; {\rm{a.s}}$

● ${ }\lim\limits_{t\to\infty}t^{-1}\; \displaystyle {\int}_0^t{K}_{\mathfrak{h}}^2(s)\: {\rm{d}}s\leq\dfrac{2{a}^2}{{c}\chi_{13}}\; \; {\rm{a.s}}$ and ${ }\lim\limits_{t\to\infty}t^{-1}\; \displaystyle {\int}_0^t {K}_{\mathfrak{m}}^2(s)\: {\rm{d}}s\leq\dfrac{2{a}_{\mathfrak{m}}^2}{{c}_{\mathfrak{m}}\chi_{14}}\; \; {\rm{a.s}}.$

where

and

The proof of this outcome closely resembles that of Lemma 2.11 presented in ^[45], and as such, it is not reiterated here.

4.   Analytical analysis of disease extinction and survival

In this section, our primary objective is to systematically establish the necessary conditions for disease extinction and persistence, taking into account the biological parameters and the influence of Lévy noise.

4.1.   Exponential extinction

In an epidemic model, the concept of stochastic extinction refers to the random and unpredictable disappearance of a mosquito-borne disease within a population. Stochasticity in epidemiological models typically arises due to factors like chance events in disease transmission and the finite size of populations. To streamline the mathematical calculations and maintain conciseness, we will employ the following definitions:

● $\chi_{15}(\xi): = -{z}_2(\xi) {1}_{\left({z}_2(\xi)\leq0\right)}+\ln\left({z}_2(\xi)+1\right){1}_{\left({z}_2(\xi) > 0\right)}.$

● $\chi_{16}(\xi): = -{z}_6(\xi) {1}_{\left({z}_6(\xi)\leq0\right)}+\ln\left({z}_6(\xi)+1\right){1}_{\left({z}_6(\xi) > 0\right)}.$

● $\chi_{17}: = {\max}\left({\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}\chi_{15}(\xi){Q}_2({\rm{d}}\xi)}, \displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}\chi_{16}(\xi){Q}_6({\rm{d}}\xi)\right).$

● $\chi_{18}: = \frac{1}{2}({a}_2{a}_6)^2\Big({a}_2^2+{a}_6^2\Big)^{-1}.$

● $\chi_{19}: = \max\big({c}+{c}_0+\varphi+{q_1}+{q_2}, {c}_{\mathfrak{m}}\big) \times \Big(-1+\mathcal{R}_{\circ}^{\frac{1}{2}}\Big)^{0, +}$, where $\mathfrak{r}^{0, +}: = 0.5\Big(\left|\mathfrak{r}\right|+\mathfrak{r}\Big)$, $\forall \mathfrak{r}\in\mathbb{R}$.

● $\chi_{20}: = \min\big({c}+{c}_0+\varphi+{q_1}+{q_2}, {c}_{\mathfrak{m}}\big)\times \Big(1-\mathcal{R}_{\circ}^{\frac{1}{2}}\Big)^{0, +}.$

● $\chi_{21}: = \chi_{19} -\chi_{20}$.

● $\chi_{22}: = \dfrac{{b}_{\mathfrak{m}}{a}_{\mathfrak{m}}}{{c}_{\mathfrak{m}}^2\left({c}+{c}_0+\varphi+{q_1}+{q_2}\right)}$.

● $\chi_{23}: = { }\mathcal{R}_{\circ}^{\frac{1}{2}}{c}_{\mathfrak{m}}^{-1}$.

● $\chi_{24}: = 0.5{d}_{\mathfrak{m}}\mathcal{R}_{\circ}^{\frac{1}{2}}\left(\dfrac{2{d}}{\chi_{13}}-1\right)^{\frac{1}{2}}+0.5\left({c}+{c}_0+\varphi+{q_1}+{q_2}\right)\mathcal{R}_{\circ}^{\frac{1}{2}}\left(\dfrac{2{d}_{\mathfrak{m}}}{\chi_{14}}-1\right)^{\frac{1}{2}}$.

● $\chi_{25}: = \chi_{17}-\chi_{18}+\chi_{21}+\chi_{24}$.

Theorem 4.1. The distinctive feature of the solution to the perturbed system (2.2) is as follows:

As a result, the infection's stochastic extinction will almost certainly transpire when $\chi_{25} < 0$.

Proof. Initially, we introduce the subsequent function:

Utilizing Itô's rule for a two-dimensional stochastic process, we derive

where

Evidently, it is apparent that

Thus,

Moreover, we can demonstrate that

We combine (4.2) and (4.3) with (4.1) to ultimately derive the following:

Alternatively, this can be reformulated as

By reordering the terms, we reach

Now, remark that $\mathcal{R}_{\circ} = \chi_{22} {b} \dfrac{{a}}{{c}}$, $\mathcal{R}_{\circ}^{\frac{1}{2}} = \chi_{23}{c}_{\mathfrak{m}}$ and $\chi_{22}\left({c}+{c}_0+\varphi+{q_1}+{q_2}\right) = {b}_{\mathfrak{m}} \dfrac{{a}_{\mathfrak{m}}}{{c}_{\mathfrak{m}}^2}$. Then we have

As a result, we deduce that

Performing a straightforward integration from 0 to $t$ and subsequently dividing by $t$ results in

where

Utilizing Hölder's inequality, we obtain the following:

Thereby

Following the same analytical approach, we obtain

Now, by using Lemma 2.5 in ^[45], we have

Ultimately, we get that

Referring to the definition of stochastic extinction of infection presented in ^[44], we can infer that under the criterion: ${\chi_{25}} < 0,$ the infection will disappear within the human and mosquito population. □

4.2.   Asymptotic survival

A vital aspect of disease eradication is its enduring presence. The primary objective of this section is to define a set of conditions that determine whether the infection can persist in both human and mosquito populations. To maintain mathematical brevity, we will introduce the following supplementary definitions:

● $\chi_{26} = {c}+\dfrac{{a}_1^2}{2}+ {\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}\Big({z}_1(\xi)-\ln\Big(1+{z}_1(\xi)\Big)\Big){Q}_1({\rm{d}}\xi)}.$

● ${\chi_{27} = \left({c}+{c}_0+\varphi+{q_1}+{q_2}\right)+\dfrac{{a}_2^2}{2}+\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}\Big({z}_2(\xi)-\ln\Big(1+{z}_2(\xi)\Big)\Big){Q}_2({\rm{d}}\xi)}.$

● ${\chi_{28} = {c}_{\mathfrak{m}}+\dfrac{{a}_5^2}{2}+\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}\Big({z}_5(\xi)-\ln\Big(1+{z}_5(\xi)\Big)\Big){Q}_5({\rm{d}}\xi)}.$

● ${\chi_{29} = {c}_{\mathfrak{m}}+\dfrac{{a}_6^2}{2}+\displaystyle {\int}_{{\mathbb R}^6\setminus\{0\}}\Big({z}_6(\xi)-\ln\Big(1+{z}_6(\xi)\Big)\Big){Q}_6({\rm{d}}\xi)}.$

● $\chi_{30} = {\dfrac{{b}{b}_{\mathfrak{m}}{a}{a}_{\mathfrak{m}}}{\chi_{26}\chi_{27}\chi_{28}\chi_{29}}}.$

Theorem 4.2. When the condition $\chi_{30} > 1$ is met, the mosquito-borne infection persists over time in both human and mosquito populations. To be more specific,

Proof. We define the following function

where $h_1, h_3, h_4 > 0$ are to be set later and $h_2 = 1$. Utilizing Itô's rule, in the context of a four-dimensional stochastic process, we obtain

where $\mathcal{L}_{\rm{s}}\mathbb{F}_s\Big({S}_{\mathfrak{h}}(t), {I}_{\mathfrak{h}}(t), {S}_{\mathfrak{m}}(t), {I}_{\mathfrak{m}}(t)\Big)$ is expressed as follows:

Then, we obtain

Thus, through the application of the arithmetic and geometric mean inequality, we can conclude that

Now, we make the following selection: $h_1 = \dfrac{{b}{b}_{\mathfrak{m}}{a}{a}_{\mathfrak{m}}}{\chi_{26}^2\chi_{28}\chi_{29}}$, $h_3 = \dfrac{{b}{b}_{\mathfrak{m}}{a}{a}_{\mathfrak{m}}}{\chi_{26}\chi_{28}^2\chi_{29}}$ and $h_4 = \dfrac{{b}{b}_{\mathfrak{m}}{a}{a}_{\mathfrak{m}}}{\chi_{26}\chi_{28}\chi_{29}^2}$. Then

Hence,

Integrating from $0$ to $t > 0$ and dividing by $t$ lead to

Hence,

Consequently,

By evaluating the limit inferior, we acquire that

This concludes the proof. □

Remark 4.1. Persistence in the mean of a mosquito-borne disease refers to the long-term or average ability of the disease to persist within a population. It typically involves mathematical models and analysis to determine whether the disease is likely to continue to exist over time, even in the presence of stochastic or random factors that may cause fluctuations. This concept helps assess the overall impact and endurance of the disease within a given population.

Remark 4.2. Theorems 4.1 and 4.2 provide thresholds for the stochastic outcomes of extinction and persistence, on average, in the model described by Eq (2.2). In essence, these theorems furnish both the necessary and nearly sufficient conditions for the disease to either die out or persist, on average.

5.   Numerical simulation

5.1.   Sensitivity analysis of the deterministic model

Sensitivity analysis assumes a crucial role in the study of dynamic systems, particularly in fields such as ecology and epidemiology ^[50]. The normalized forward sensitivity index quantifies the relative change in the basic reproduction number $\mathcal{R}_{0}$ concerning variations in parameter values. It provides comprehensive insights into the model's resilience to such alterations. Moreover, this index serves as a tool for identifying parameters with a substantial impact on the basic reproduction numbers, guiding the formulation of targeted epidemiological intervention strategies. To be more specific, the normalized forward sensitivity index is expressed as the ratio of the relative change in the basic reproduction numbers ($\mathcal{R}_{0}$) to the relative change in the parameter (ß). Assuming that $\mathcal{R} _{0}$ is differentiable with respect to the parameter, the formula is given by:

where, $\mathcal{R}_{0}$ is considered a function of the parameter ß. Given that $\mathcal{R}_{0}$ is a rational function of model parameters, the normalized forward sensitivity index is applicable to all model parameters present in the explicit formula defining $\mathcal{R}_{0}$.

Utilizing the formula (5.1) across all parameters of the model, the observed outcomes are as follows:

(ⅰ) $a, a_{m}, b,$ and $b_{m}$ have positive index values with $a = a_{m} = b = b_{m} = 1$. This indicates that any alteration in the values of these parameters directly affects $\mathcal{R}_{0}$, leading to either an increase or decrease in its value.

(ⅱ) $c, c_{0}, c_{m}, q_{1}, q_{2},$ and $\varphi$ have negative index values with:

indicating that an increase in their values leads to a decrease in $\mathcal{R}_{0}$. Moreover, the values of $\mathcal{R}_{0}$ remain unaffected by the parameters $f_{1}, f_{2}, g_{1},$ and $g_{2}$.

To compute the sensitivity indices, we utilize the data presented in the second column of Table 1. Examination of Table 2 and Figure 1 reveal that a $10\%$ increase (or decrease) in the values of $a, a_{m}, b,$ and $b_{m}$ results in a corresponding $10\%$ increase (or decrease) in $\mathcal{R}_{0}$, for each parameter. Conversely, a 10% increase in the values of $c, c_{0}, c_{m}, q_{1}, q_{2},$ and $\varphi$ leads to a reduction in $\mathcal{R}_{0}$ by $10.7849\%, 0.0471\%, 20\%, 0.6436\%, 0.675\%,$ and $7.8493\%$, respectively.

It is crucial to note that the correlation between the rate in which the infected human population is put into quarantine $\varphi$ and $\mathcal{R} _{0}$ is negative. This signifies that as the rate of progression increases, $\mathcal{R}_{0}$ tends to decrease. This underscores the significance of the importance of quarantining in reducing the infection rate

5.2.   Simulation techniques and verification of the theoretical findings

Creating stochastic processes using computer methods requires two distinct discretization techniques. Initially, we must address the discretization of the time parameter, followed by an approximation of random variables using artificially generated finite time series datasets. In the case of a Lévy process, which exhibits stationary and independent increments, the most straightforward approach to tackling the challenge of simulating it solely for discrete time points is analogous to generating random numbers from an infinitely divisible distribution.

In this section, we will explore a technique for simulating GTS distributions and tempered stable processes. While there are various approaches available for simulating Lévy processes, many of them are not well-suited for simulating tempered stable processes because of the intricate nature of their Lévy measure.

Let $\{S_{1, j}\}_{j\geqslant 1}$ represent a sequence of independent and identically distributed (i.i.d.) random variables in the real numbers, following the distribution (2.3). Additionally, consider $\{S_{2, j}\}_{j\geqslant 1}$ and $\{S_{3, j}\}_{j\geqslant 1}$ as i.i.d. sequences of uniform random variables within the intervals $(0, 1)$ and $(0, T)$, respectively. Also, let $\{S_{4, j}\}_{j\geqslant 1}$ and $\{S_{5, j}\}_{j\geqslant 1}$ be i.i.d. sequences of random variables following exponential distribution with a rate coefficient of $1$. we assume that all mentioned random variables are mutually independent. Now, we let

Noticeably, $\{S_{6, j}\}$ can be regarded as a Poisson point process on the interval ${\mathbb R}_+$ with random intensity measure. In reference to the theory outlined in ^[37], if $\alpha\in(0, 2)$, then

converges almost surely and uniformly for t within the interval $[0, T]$ to a Lévy process, where

Ultimately, we can formulate a method for generating a GTS process with specified parameters at discrete time instances $t_i$, where $\{t_i\}_{i\in[0, I]}$ represents a partition of the interval $[0, T]$ with uniformly sized subintervals and mesh $\Delta t = T/I$, $I\in \mathbb{N}$. Then, we use the following algorithm:

(1) Select a specific time duration T and create a division of the time interval $[0, T]$ into $I$ equally sized segments.

(2) Select and fix a number $N$.

(3) Numerically replicate or emulate independent quantities $\{S_{i, j}\}$, $i\in\{1, \cdots, 6\}$ of range $N$.

(4) Determine the value of $G_t$.

Using the aforementioned algorithm, we can generate the complete path of a GTS process associated with system (2.2).

Remark 5.1. The parameter $\alpha$, known as the stability index, provides valuable insights into the tail behavior of an tempered $\alpha$-stable distribution. When $\alpha$ is less than 1, it results in heavy tails, signifying an increased likelihood of extreme events. Conversely, when $\alpha$ is greater than 1, the tails become lighter, resembling a distribution that is closer to normal. The clarity of this can be demonstrated in the visual representations presented in Figure 2.

In this following, our objective is to assess the accuracy of the results outlined in Theorems 4.1 and 4.2 with a specific focus on the impact of GTS distribution on the dynamics of mosquito-borne infections. In the ensuing examples, we initialize our compartmental model (2.2) with the following initial data: ${S}_{\mathfrak{h}}(0) = 500$, ${I}_{\mathfrak{h}}(0) = 15$, ${Q}_{\mathfrak{h}}(0) = 2$, ${R}_{\mathfrak{h}}(0) = 2$, ${S}_{\mathfrak{m}}(0) = 150$, and ${I}_{\mathfrak{m}}(0) = 10$. We assume that each time unit represents one day, and the population size is measured in millions of individuals. For stability index, we take $\alpha = 1.5$. To facilitate comparison, we illustrate three distinct types of trajectories: deterministic (without any noises), solutions incorporating jumps with standard measure (with standard distribution), and stochastic trajectories with jumps and GTS distribution (with GTS distribution). This allows us to explore various potential scenarios. It is important to note that in the following three experiments, our assumptions hold true.

5.2.1.   Test 1: Special case of extinction of the infection

In this particular context, as our first procedural step, we assign numerical values to our system parameters, aligning them with the data provided in the second column of Table 1. Additionally, we set $a_1 = 0.1$, $a_2 = 0.21$, $a_3 = 0.17$, $a_4 = 0.01$, $a_5 = 0.13$, and $a_6 = 0.16$. For the jump intensities, we adopt the following function: ${z}_L(\xi) = \frac{-u_L \xi}{0.5+\xi^2}$, $L\in\{1, \cdots, 6\}$ where $\xi = 0.3,$ with the specific values $u_1 = u_3 = u_4 = u_5 = 0.02$ and $u_2 = u_6 = 0.03$. In this context, we have achieved $\mathcal{R}_{\circ} = 1.0266 > 1$ and $\chi_{25} = -0.1096 < 0$. Consequently, we confirm that the essential conditions described in Theorem 4.1 are met. To empirically confirm this result, we illustrate three separate system trajectories corresponding to (2.2) in Figure 3. In this latter, we emphasize the significance of incorporating the GTS distribution. As demonstrated, in a system influenced by the standard Lévy distribution (model (2.2) with standard Lévy jumps), the disease exhibits continuous persistence within both the human and mosquito populations. However, when the GTS distribution is employed, the disease is seen to extinguish. This underscores the critical role played by our GTS distribution model, which addresses a key aspect that remains unresolved within the standard Lévy jump framework.

Lévy processes, when utilizing GTS distribution, exhibit heavy-tailed characteristics, signifying a heightened likelihood of extreme events compared to Gaussian (normal) distributions. In the context of infection transmission, this implies a significant probability of experiencing a sudden surge in infections over a short timeframe, which can lead to rapid disease propagation and, potentially, its extinction. The inherent stochastic nature of Lévy jumps with the GTS distribution introduces the possibility of a remarkable surge in the number of infected individuals, potentially culminating in a stochastic extinction event. Should this surge result in a substantial portion of the population becoming infected, and if the circumstances are unfavorable for sustained transmission, the disease may self-extinguish.

5.2.2.   Test 2: Complete extinction of infection

Building upon the results of Test 1, we have made slight adjustments to specific parameters in order to simulate a scenario of complete extinction. We assign numerical values to our system parameters by aligning them with the data provided in the third column of Table 1. In addition, we fix the values of $a_L = 0.3$, $L\in\{1, \cdots, 6\}$. Regarding the jump intensities, we employ the following function: ${z}_L(\xi) = \frac{-u_L \xi}{0.5+\xi^2}$, where $L$ takes values from $1$ to $6$. Here, we set $\xi = 0.3$ and assign specific values such that $u_1 = u_3 = u_4 = u_5 = 0.02$, and $u_2 = u_6 = 0.03$. Here, $\mathcal{R}_{\circ} = 0.7370 < 1$ and $\chi_{25} = -0.4266 < 0$. As a result, we have verified that the prerequisites outlined in Theorem 4.1 have been met. To validate this finding through numerical analysis, we have depicted three distinct types of trajectories associated with system (2.2) in Figure 4. It is readily apparent that the model consistently converges to a state devoid of infection for all three trajectory types. To be more specific, when we consider the initial data mentioned earlier, we observe that $S_h(t)$ stabilizes at a constant value of $140$ over time. Similarly, $S_m(t)$ reaches an equilibrium point at $340$. However, for $I_h(t)$, $Q_h(t)$, $R_h(t)$, and $I_m(t)$, irrespective of the trajectory type, the solutions ultimately extinguish after a certain duration. This behavior precisely illustrates the concept of stochastic extinction. In this context, we are referring to complete extinction, as all trajectories (i.e., those without noise, with standard distribution, and with GTS distribution) exhibit identical behavior.

5.2.3.   Test 3: Persistence of infection

Let's now explore the scenario of persistent infection. In this experiment, we consider the following parameter values: $a_1 = 0.011$, $a_2 = 0.021$, $a_3 = 0.017$, $a_4 = 0.0101$, $a_5 = 0.031$, and $a_6 = 0.0106$. In terms of the jump intensities, we utilize the function ${z}_L(\xi) = \frac{-u_L \xi}{0.5+\xi^2}$, with $\xi = 0.3$. Specifically, we set $u_1 = u_3 = u_4 = u_6 = 0.01$, and $u_2 = u_5 = 0.02$. By employing the numerical values from the last column of Table 1, we can readily confirm that our hypothèses remain valid, and $\chi_{30}$ exceeds one $(\chi_{30} = 1.2530 > 1)$. Therefore, in accordance with Theorem 4.2, we can confidently assert that our model persists on average, which is consistent with the patterns observed in Figure 5. Notably, the endemic equilibrium of the deterministic version no longer serves as the steady state for the stochastic model (2.2). Consequently, over an extended period, the intensity of noise influences the degree to which the solution fluctuates around the deterministic equilibrium states. To provide more precision, the temporal average closely aligns with the endemic equilibrium, especially for low noise intensities. Above all, this observation underscores the importance of incorporating environmental fluctuations into the biological dynamics. To draw a meaningful comparison between standard jumps fluctuations and jumps with GTS distribution, we remark in 5 that in the case of jumps with GTS distribution, the inherent volatility undergoes temporal variations and exhibits heavier-than-normal tails. Furthermore, stochastic volatility occasionally experiences substantial upward jumps and clusters at high levels, contrasting with the behavior observed in the standard scenario.

6.   Conclusions and further remarks

The significance of an epidemic model lies in its ability to comprehensively depict the characteristics of the modeled disease, encompassing a wide array of its biological aspects. With the goal of enhancing our understanding of mosquito-borne diseases, we have introduced a novel epidemic model. This model takes into consideration two key factors: The impact of quarantine measures on human populations and the influence of random and severe environmental fluctuations deriven by GTS distribution. We have employed a compartmental modeling approach to construct this model, representing it as a system of interconnected stochastic differential equations driven by Lévy noise. For this resulting model, we have rigorously established its mathematical soundness, biological plausibility, and its behavior over extended periods in the absence of the disease. Moreover, by carefully selecting appropriate stochastic parameters, we have identified conditions under which the infection can persist or be eradicated within both mosquito and human populations. Our model encompasses a variety of transitions between compartments, making it relevant for understanding and studying mosquito-borne diseases such as zika virus, west Nile virus, chikungunya virus, dengue fever, and more.

Nevertheless, it is essential to acknowledge that our research has certain constraints, prompting intriguing questions that merit deeper exploration. These limitations are summarized as follows:

● Existence of a stationary sistribution for model (2.2): As previously mentioned, it's important to note that the stochastic system described in (2.2) lacks an endemic state. Therefore, we need to employ an alternative concept of stochastic stability. Specifically, it becomes imperative to verify the presence of a stationary distribution for (2.2). From a biological standpoint, the existence of such a distribution implies the persistence of the disease within both mosquito and human populations. Consequently, it becomes intriguing to identify the conditions under which such a distribution can be established. To the best of our knowledge, in the context of mosquito-borne epidemic models driven by Lévy noise and GTS distribution, this remains an open question.

● Parameter identification for model (2.2): It is important to note that the results presented in this paper pertain to the solution of the direct problem. In simpler terms, we have operated under the assumption that all parameters are known a priori. However, when it comes to adapting the model for a specific mosquito-borne disease, it becomes essential to address the identification problem associated with model (2.2). Specifically, this involves determining the appropriate values for the model's parameters and the stochastic noise, given observations of the total population over a defined time period, in order to achieve the best possible fit to the observed data. To the best of our knowledge, within the context of epidemic models driven by Lévy noise and GTS distribution, this remains an unresolved challenge.

Due to the complexity of these questions and their merits to be treated independently, we will leave them for our next future works.

Author contributions

Yassine Sabbar: Conceptualization, Writing Original Draft, Software, Formal Analysis; Aeshah A. Raezah: Conceptualization, Software, Validation. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors extend their appreciation to the Deanship of Research and Graduate Studies at King Khalid University for funding this work through Large Research Project under grant number RGP2/174/45.

Conflict of interest

The corresponding author states that there is no conflict of interest.