Optimizing microbe-infected mosquito release: a stochastic model for malaria prevention

Malaria remains a critical public health challenge in Africa, demanding innovative control strategies. This study introduces a novel approach using Microsporidia MB-infected mosquitoes and stochastic optimal control within a Lévy process framework to regulate mosquito release strategies. The primary goal is to optimize Microsporidia MB prevalence within mosquito populations to disrupt Plasmodium transmission to humans. By incorporating Lévy noise into the modeling process, we capture the inherent randomness of mosquito dynamics, improving intervention accuracy. The model, guided by the Hamilton–Jacobi–Bellman (HJB) equation, optimizes release protocols while accounting for key environmental factors like seasonality and temperature fluctuations. Results show that intervention success depends on local climatic conditions, underscoring the need for flexible, region-specific strategies in malaria-endemic areas. Focus regions include Kenya, Ghana, Niger, and Benin, where Microsporidia MB has been confirmed. Findings suggest that targeted mosquito releases could significantly reduce malaria transmission, offering valuable insights for public health efforts.

1 Introduction

Malaria is a mosquito-borne disease that results from an individual being infected by any of the five Plasmodium species known as P. falciparum, P. vivax, P.ovale, P. malariae and P.knowlesi [1–3]. Approximately 50 mosquito species from the genus, Anopheles have been identified as major vectors of the disease [2]. According to the World Health Organization [4], malaria deaths have been estimated at more than 600,000 during the year 2021. The 2022 World Malaria Report states that the disease continues to persist in the African region more than any other region of the world as the region accounts for the highest cases and deaths, especially among children.

To control the spread of malaria, control and intervention measures have been applied over the years and these efforts have helped in lowering the prevalence of malaria in most regions across the globe. As an effort to lower the disease infection and transmission, insecticide-treated nets (ITNs) and indoor residual spraying (IRS) are primarily used in order to control malaria vectors [5]. However, the rising issues of resistance threaten the progress and efforts achieved [6].

The World Health Organization has put forward its support for the development of improved, innovative, and sustainable methods to reduce and eliminate malaria. One of the aforementioned involves the use of endosymbionts as plasmodium transmission blockers in mosquitoes. In literature, the bacteria Wolbachia has been proposed and explored as an endosymbiont for malaria control [7, 8]. Another possible endosymbiont is a microsporidian denoted Microsporidia MB that has been discovered to impair plasmodium transmission [9]. This Microsporidia MB species has been identified in association with Anopheles arabiensis, Anopheles gambiae, and Anopheles coluzzii, as evidenced by research conducted in Kenya, Ghana, and, more recently, in Niger and the Republic of Benin [9–12]. In the pursuit of effective malaria control strategies, researchers have explored the intriguing role of microsporidia-infected mosquitoes. Fox and Weiser's investigation into Anopheles gambiae infected with Nosema stegomyiae shed light on natural infections in Liberia [13]. Meanwhile, laboratory studies by other scientists focused on Plistophora culicis infections in Culex pipiens fatigans mosquitoes, revealing altered mosquito biology and reduced Plasmodium development [14]. Additionally, a recent study examined the effects of Nosema algerae on Plasmodium yoelii nigeriensis in Anopheles stephensi, emphasizing the intricate interactions within the mosquito host [15]. Notably, Bano's work highlighted the partial inhibitory effect of P. culicis on the sporogonic cycle of P. cynomolgi [16]. While these studies don't directly prescribe a solution, they contribute valuable knowledge for advancing biocontrol strategies against malaria transmission.

Although information on Microsporidia MB is limited, past literature suggests that microsporidia could serve as a promising malaria vector control strategy, warranting further investigation. The potential for these advancements is both challenging and exciting; however, practical application requires careful consideration of deploying endosymbiont-infected mosquitoes into natural environments. Various mathematical models have been employed to study vector-borne diseases, providing valuable insights into effective control strategies. For instance, models depicting the release of Wolbachia-infected mosquitoes have demonstrated practical guidance for implementing release strategies to achieve desirable mosquito populations in the field [17–20]. Additionally, optimal control strategies using sterile mosquitoes [21–26] and genetically modified mosquitoes have been explored [27, 28]. Furthermore, recent studies have employed fractional differential equations to model vector-borne diseases such as dengue fever, capturing complex behaviors and memory effects in disease transmission dynamics [29–33]. While these approaches provide valuable insights, they often do not fully account for the inherent variability and uncertainty driven by environmental factors. This is where stochastic methods become essential. By incorporating stochastic elements, such as Stochastic Differential Equations (SDEs), we can better capture the unpredictable dynamics of mosquito populations, influenced by environmental fluctuations. SDEs have been widely used across various fields to model systems under uncertainty. In epidemiology, they are pivotal for understanding disease transmission dynamics and informing control strategies, making them highly relevant for developing effective frameworks for malaria vector control. Beyond epidemiology, SDEs have proven their utility in finance, engineering, ecology, and even agriculture, where they are used to model complex systems influenced by random factors [34–38]. By integrating stochastic elements with the dynamics of Microsporidia MB, our approach aims to create a more realistic and effective framework for malaria vector control, better equipped to handle the variability of real-world environments.

Building upon previous work in vector control, this study introduces a novel approach by utilizing Microsporidia MB-infected mosquitoes, rather than transgenic mosquito populations, to block malaria transmission. The primary objective is to optimize this biological intervention by incorporating stochastic elements into the modeling of mosquito population dynamics, with the ultimate goal of enhancing microbe prevalence in mosquito populations to disrupt Plasmodium transmission to humans. A key innovation of our methodology is the use of Lévy noise, which captures the inherent randomness and unpredictability in mosquito population fluctuations. This stochastic framework is crucial for accurately modeling how different release frequencies and patterns of infected mosquitoes affect overall population dynamics, overcoming the limitations of traditional deterministic models. Additionally, we expand our analysis spatially by mapping optimal infection prevalence rates across selected West African countries-namely Niger, Ghana, and Benin-where malaria is endemic, and also including Kenya. This geographical focus addresses the need for tailored vector control strategies that account for regional variability in mosquito behavior, infection prevalence, and environmental conditions. By using advanced simulations and mapping techniques, our research highlights the complexities of controlling mosquito populations, particularly under varying environmental and seasonal conditions. We examine stochastic growth rates, infection prevalence ratios, and uninfected mosquito proportions to deliver insights into the effectiveness of mosquito release strategies. This study is particularly innovative in its application of stochastic processes to epidemiological modeling and its focus on Microsporidia MB-based biological interventions, offering a promising alternative to traditional transgenic methods.

This research applies stochastic optimal control and Lévy processes to develop an optimized strategy for releasing Microsporidia MB-infected mosquitoes to block malaria transmission. By integrating mathematical modeling with ecological and biological insights, we address the complexities of mosquito population dynamics and environmental variability. Our approach formulates a stochastic control problem, solved using the Hamilton–Jacobi–Bellman (HJB) equation under jump diffusions, to determine the optimal mosquito release strategy. In Section 2, we define the theoretical framework by introducing the Lévy process and its applications, outlining stochastic control with Lévy jumps, and formulating the dynamics of mosquito populations influenced by both white and Lévy noise. This section also establishes the framework for the stochastic optimal mosquito release problem, detailing the data sources for temperature measurements and model validation. Section 3 progresses to the analytical derivation of optimal release strategies, including Python-based simulations, visualization maps based on temperature variations, and validation of these strategies. Section 4 discusses the implications of our findings within the broader context of malaria control and concludes by summarizing our research outcomes, along with suggesting future directions for vector-borne disease management.

2 Materials and methods

Our methodology applies a Lévy process-based mathematical model to analyze the population dynamics of Anopheles mosquitoes, distinguishing between Microsporidia MB-infected and uninfected groups. The model combines traditional population dynamics with stochastic elements, such as white noise and Lévy noise, to address environmental fluctuations and unpredictable changes. Central to our approach is the strategic release of infected mosquitoes to influence the wild population and enhance symbiont prevalence, thus supporting malaria control efforts. Key parameters of the model, including growth rates, mortality, and environmental factors, are essential for understanding and optimizing this strategy [39].

2.1 Area of study

This research examines the population dynamics of Anopheles mosquitoes in selected regions of West Africa-namely Niger, Ghana, and Benin Republic and Kenya, East Africa (Figure 1). These areas are home to key malaria vectors such as Anopheles gambiae, Anopheles coluzzii, Anopheles funestus, and Anopheles arabiensis. In these mosquito species, symbiotic microbes known to inhibit the malaria parasite Plasmodium have been identified. The discovery of these microbial inhibitors presents a promising strategy for controlling malaria transmission.

Figure 1

Figure 1. Selected regions of West Africa, including Niger, Ghana, and Benin, along with Kenya in East Africa, are the focus of this study, as Microsporidia MB has been identified in these areas.

Given the significant variability in environmental factors such as temperature across these geographic areas, our approach tailors the release strategies to align with local ecological conditions. This regional customization aims to optimize the conditions under which microbe-infected mosquitoes are released, enhancing the prevalence and effectiveness of the microbial symbionts in reducing malaria transmission. By examining how different environmental settings impact the success of these strategies, we aim to identify the most suitable environments for effective mosquito releases, thereby maximizing the potential for controlling malaria in diverse settings.

2.2 Data

2.2.1 Temperature data

We utilized the Bio1 bioclimatic variable dataset from WorldClim, offering annual mean temperature data spanning the period 1970–2000. This dataset, characterized by a spatial resolution of 30 s (~1 km²), served as a foundational component for our study. To tailor the dataset to our study's specific requirements, we employed QGIS to clip the raster layer to our designated study sites, ensuring a focused analysis of the relevant geographic areas.

2.2.2 Validation data

For the validation process, Microsporidia-infected mosquito data was collected in Kenya from the years 2022 to 2023 after the rainy seasons. The data specifically focus on Anopheles mosquitoes-An. Gambiae, An. Funestus, An. Arabiensis, and An. Coluzzii known for transmitting the Plasmodium parasite across diverse regions in Africa. The alignment of this information with our mapped regions of malaria transmission potential enhances the robustness of our validation process.

2.3 Preliminaries

In this section, we introduce the fundamental concepts of stochastic optimal control, focusing on the application of Lévy processes and jump-diffusion methods in the modeling of evolutionary dynamics in population biology. We also provide preliminary results that demonstrate how our proposed model integrates these stochastic elements to effectively capture the complex variabilities and interactions within biological populations.

2.3.1 Lévy Process

Definition 2.1. Let us define a Lévy process on a filtered probability space $(\Omega ,\mathcal{F},{({\mathcal{F}}_t)}_{t\ge 0},\mathbb{P})$, where Ω is the sample space, $\mathcal{F}$ is the σ-algebra, ${({\mathcal{F}}_t)}_{t\ge 0}$ is the filtration representing the evolution of information over time, and ℙ is the probability measure. A stochastic process {_Xt}t≥0, taking values in ℝ and adapted to the filtration ${\mathcal{F}}_t$, is termed a Lévy process if it satisfies the following three conditions:

1. Initialization: X₀ = 0 almost surely, establishing the process at the origin with probability one.

2. Stochastic Continuity: For any ϵ > 0 and for any s ≥ 0, $\underset{t\to s}{lim}\mathbb{P}(|X_t-X_s|>ϵ)=0$, ensuring that the process is continuous in probability, allowing for jumps but not “wild” behavior.

3. Independent and Stationary Increments: The process exhibits increments that are independent and stationary; that is, for any 0 ≤ s < t, the increment X_t − X_s is independent of ${\mathcal{F}}_s$ and the distribution of X_t − X_s depends only on t − s, not on the actual values of t or s. This property indicates that the process evolves in a way where the future increment is independent of the past, with its distribution reliant solely on the length of the time interval.

2.3.1.1 The dynamics of Itô-Lévy processes

Within the framework of stochastic calculus, the Itô-Lévy process stands out for its ability to model complex systems influenced by both continuous fluctuations and discrete jumps. Herein, we present a foundational theorem characterizing the evolution of such processes:

Theorem 2.2 (Evolution of Itô-Lévy Processes). Consider a process X(t), taking values in ℝ, defined on a probability space. The dynamics of X(t) as an Itô-Lévy process are captured by the differential equation:

$dX(t)=\lambda (t,\omega )dt+\sigma (t,\omega )dW(t)+{\displaystyle {\int }_{ℝ}}\gamma (t,y,\omega )Ñ(dt,dy),$

where Ñ(dt, dy) is a compensated Poisson random measure defined by:

$Ñ(dt,dy)=\{\begin{array}{ll}N(dt,dy)-\vartheta (dy)dt,& \text{ for }\left|y\right|<R,\\ N(dt,dy),& \text{ for }\left|y\right|\ge R,\end{array}$

with R ∈ [0, ∞] demarcating the threshold for jump sizes. Let h : ℝ² → ℝ be a twice continuously differentiable function, and define a new process Z(t) = h(t, X(t)). The differential of Z(t), embodying the Itô-Lévy dynamics, is then given by:

$\begin{array}{l}dZ(t)=\frac{\partial h}{\partial t}(t,X(t))dt+\frac{\partial h}{\partial x}(t,X(t))\left[\lambda (t,\omega )dt+\sigma (t,\omega )dW(t)\right]\\ +\frac{1}{2}{\sigma }^2(t,\omega )\frac{{\partial }^{2}h}{\partial x^2}(t,X(t))dt\\ +{\displaystyle {\int }_{\left|y\right|<R}}[h(t,X(t^{-})+\gamma (t,y,\omega ))-h(t,X(t^{-}))\\ -\frac{\partial h}{\partial x}(t,X(t^{-}))\gamma (t,y,\omega )]\vartheta (dy)dt\\ +{\displaystyle {\int }_{ℝ}}[h(t,X(t^{-})+\gamma (t,y,\omega ))-h(t,X(t^{-}))]N(dt,dy).\end{array}$

This theorem elucidates the nuanced behavior of Itô-Lévy processes, highlighting their capacity to integrate both gradual trends and abrupt shifts in state, thereby offering a robust model for systems exhibiting complex dynamics.

2.3.1.2 Control state dynamics and the generator in Lévy-driven systems

The dynamics of a controlled system influenced by Lévy processes can be characterized by its state equation and the associated infinitesimal generator, which together describe how the system evolves under a given control policy. The state equation incorporates both continuous and jump components, while the generator captures the system's local behavior under the control policy.

Definition 2.3 (Controlled Lévy Process State Equation). Consider a controlled stochastic system where the state dynamics are governed by a Lévy-driven Stochastic Differential Equation (SDE) in the form:

$\begin{array}{r}dX(t)=\lambda (t,X(t),u(t))dt+\sigma (t,X(t),u(t))dW(t)\\ +{\displaystyle {\int }_{ℝ}}\gamma (t,X(t^{-}),z,u(t))Ñ(dt,dz),& (1)\end{array}$

where X(t) ∈ ℝⁿ represents the state of the system at time t, $u(t)\in \mathcal{U}$ denotes the control action applied at time t, chosen from a set of admissible controls $\mathcal{U}$, W(t) is a standard Brownian motion, and Ñ(dt, dz) is the compensated Poisson random measure reflecting the jumps.

Definition 2.4 (Infinitesimal Generator of a Controlled Lévy Process). The infinitesimal generator ${\mathcal{A}}^u$ associated with a controlled Lévy process under policy u(t) is defined for a sufficiently smooth function h : ℝⁿ → ℝ as:

${\mathcal{A}}^{u}h(x)={\displaystyle \underset{\delta \to 0}{lim}}\frac{\text{E}\left[h(X(t+\delta ))|X(t)=x,u(t)=u\right]-h(x)}{\delta },$

which, under the Lévy process dynamics, can be explicitly written as:

$\begin{array}{l}{\mathcal{A}}^{u}h(x)={\displaystyle \sum _{i=1}^n{\lambda }_i}(t,x,u)\frac{\partial h}{\partial x_i}(x)+\frac{1}{2}{\displaystyle \sum _{i,j=1}^n{\left[\sigma {\sigma }^T\right]}_{ij}}(t,x,u)\frac{{\partial }^{2}h}{\partial x_i\partial x_j}(x)\\ +{\displaystyle {\int }_{{ℝ}^n}[h(x+\gamma (t,x,z,u))}\\ -h(x)-{\displaystyle \sum _{i=1}^n{\gamma }_i}(t,x,z,u)\frac{\partial h}{\partial x_i}(x)]\vartheta (dz).\end{array}$

Here, λ_i, σ, and γ_i are the components of the drift, diffusion, and jump intensity functions, respectively, and ϑ is the Lévy measure that characterizes the jump distribution.

The infinitesimal generator ${\mathcal{A}}^u$ plays a crucial role in the formulation of dynamic programming and the Hamilton–Jacobi–Bellman (HJB) equation for controlled Lévy processes. It encapsulates the expected rate of change of a function of the system's state under a specific control policy, providing a vital link between the stochastic dynamics and the optimization of the control policy.

2.3.1.3 Dynamic programming for controlled Lévy diffusions

Dynamic programming forms the cornerstone of solving stochastic control problems, particularly when the system's dynamics exhibit both continuous fluctuations and discrete jumps. Herein, we present a concise definition tailored to controlled Lévy diffusions:

Definition 2.5. Let us consider a controlled Lévy diffusion process Z(t), initiated at Z(0) = z ∈ ℝ^k, evolving according to:

$\begin{array}{r}dZ(t)=b(Z(t),u(t))dt+\sigma (Z(t),u(t))dW(t)\\ +{\displaystyle {\int }_{{ℝ}^k}}\gamma (Z(t^{-}),u(t^{-}),y)N_0(dt,dy),& (2)\end{array}$

where the drift b, volatility σ, and jump intensity γ are functions defining the system's dynamics under a control policy u(t) from a set of admissible controls U. N₀(dt, dy) denotes the compensated Poisson measure, capturing the jump behavior.

The objective is encapsulated in the performance criterion J^u(z), expressed as:

$J^u(z)={\text{E}}^z\left[{\displaystyle {\int }_0^{{\tau }_S}}h(Z(t),u(t))dt+i(Z({\tau }_S))I_{\left\{{\tau }_S<\infty \right\}}\right],$

aiming to optimize over the course until exit time ${\tau }_S=inf\left\{t>0;Z^u(t)\notin S\right\}$, with h and i being the cost functions.

A control u is deemed admissible if it ensures a unique strong solution of the diffusion process for any initial condition in S and satisfies the integrability condition for the negative parts of h and i.

The essence of the stochastic control problem lies in determining the optimal value function Φ(y) and identifying an optimal control policy $u^{*}\in \mathcal{U}$, fulfilling:

$\Phi (y)={\displaystyle \underset{u\in \mathcal{U}}{sup}}{J}^u(y)=J^{u^{*}}(y).$

For Markov controls, where u = u(z), the diffusion Z(t) under u transforms into a Lévy diffusion characterized by the generator A^u, defined as:

$\begin{array}{l}{A}^u\varphi (z)={\displaystyle \sum _{i=1}^k}{b}_i(z,u(z))\frac{\partial \varphi }{\partial z_i}(z)+\frac{1}{2}{\displaystyle \sum _{i,j=1}^k}{(\sigma {\sigma }^T)}_{ij}(z,u(z))\frac{{\partial }^2\varphi }{\partial z_i\partial z_j}(z)\\ +{\displaystyle {\int }_{{ℝ}^k}}\Delta \varphi (z;y)\vartheta (dy),& (3)\end{array}$

where Δϕ(z; y) = ϕ(z + γ(z, u(z), y))−ϕ(z)−∇ϕ(z)·γ(z, u(z), y), incorporating the jump effects through the Lévy measure ϑ.

Theorem 2.6. Hamilton–Jacobi–Bellman criterion for jump diffusion systems. Consider a controlled jump diffusion system within a domain S ⊂ ℝ^k, subject to the following conditions:

1. Let ϕ be a function such that ϕ ∈ C²(S) for the interior of S and continuous on $\overline{S}$, the closure of S. It must satisfy:

(a) For every control v ∈ U and state z ∈ S, the inequality

$A^v\varphi (z)+h(z,v)\le 0$

holds, where A^v denotes the infinitesimal generator under control v.

(b) The process exits S through its boundary ∂S with certainty if the exit time τ_S is finite, and

${\displaystyle \underset{t\uparrow {\tau }_S}{lim}}\varphi (Z(t))=i(Z({\tau }_S))I_{\left\{{\tau }_S<\infty \right\}},$

almost surely, for all admissible controls.

(c) The expectation of the absolute value of ϕ at the exit time and the accumulated generator and gradient terms over the interval [0, τ_S] is finite, i.e.,

$\begin{array}{l}{\text{E}}^z\left[|\varphi (Z(\tau ))|+{\displaystyle {\int }_0^{{\tau }_S}}(|A\varphi (Z(t))|)dt\right]<\infty ,\end{array}$

for all admissible controls u and stopping times τ.

(d) The negative part of ϕ(Z(τ)) is uniformly integrable up to τ_S.

Then, ϕ(z) serves as an upper bound for the value function Φ(z) for all z ∈ S, i.e,

$\varphi (z)\ge \Phi (z)\forall y\in S.$

2. Furthermore, if there exists a feedback control û such that for any y ∈ S,

(a) Applying û(y) yields A^ûϕ(y)+h(y, û(y)) = 0,

(b) The process under û maintains the uniform integrability of ϕ up to τ_S.

Then, the strategy û, when applied as u^*(t) = û(Z(t⁻)), constitutes an optimal control. Consequently, we have that ϕ(y) equals the optimal value function Φ(y), and thus satisfies:

$\varphi (y)=\Phi (y)=J^{u^{*}}(y)\forall y\in S.$

Proof. See proof in [39].

2.4 Model formulation

2.4.1 Population dynamics of anopheles mosquitoes X(t)

We use a stochastic differential equation based on ecological modeling principles to explore the dynamics of Anopheles mosquito populations. Initially adhering to a fundamental growth model with death rate λ, expressed as

$\begin{array}{ll}\frac{dX(t)}{dt}=rX\left(t\right)\left(1-\frac{X\left(t\right)}{K}\right)-\lambda X\left(t\right),& (4)\end{array}$

where r signifies the intrinsic growth rate and K denotes the carrying capacity, we extend our exploration to encapsulate more nuanced dynamics. That is integrating stochastic perturbation into Equation 4 such as :

$\begin{array}{ll}X:-\lambda \to -\lambda +\sigma \frac{dW\left(t\right)}{dt}.& (5)\end{array}$

Therefore the population dynamics of Anopheles become:

$\begin{array}{ll}dX\left(t\right)=\left(r\left(1-\frac{X\left(t\right)}{K}\right)-\lambda \right)X\left(t\right)dt+\sigma X\left(t\right)dW\left(t\right),& (6)\end{array}$

where σ captures random fluctuations driven by the infinitesimal increments of a Wiener process or simply white noise dW(t). This model provides a foundational understanding of how the population grows and is constrained by environmental capacity and natural mortality.

However, to account for additional complexities such as the release of mosquitoes and other sudden ecological perturbations, we adapt this model by incorporating a more dynamic representation of the growth rate. This adaptation is crucial to reflect more accurately the variability and unpredictability inherent in real-world ecological systems. Consequently, the logistic growth term $r\left(1-\frac{X\left(t\right)}{K}\right)$ is approximated by a time-dependent factor $\tilde{\beta }\left(t\right)$, such that,

$\begin{array}{ll}\tilde{\beta }(t)dt=\beta (t)dt+\zeta {\displaystyle {\int }_{ℝ}}yÑ(dt,dy),& (7)\end{array}$

where β(t) is modeled using the Cox-Ingersoll-Ross (CIR):

$\begin{array}{ll}d\beta (t)=a(b-\beta (t))dt+c\sqrt{\beta (t)}dV(t),\beta (0)={\beta }_0.& (8)\end{array}$

Therefore (6) becomes:

$\begin{array}{c}dX(t)=(X(t)·\beta (t)-\lambda X(t))dt+\sigma X(t)dW(t)\\ +\zeta X(t){\displaystyle {\int }_{ℝ}}yÑ(dt,dy),X(0)=X_0.& (9)\end{array}$

The stochastic nature of β(t) captures intrinsic growth adjustments and external shocks, thus refining the population dynamics model to more accurately reflect sudden changes and interventions. Additionally, the correlation between the stochastic processes dW(t) and dV(t), quantified by the correlation coefficient ρ, introduces further complexity. This enhances the model's capability to represent the interconnected ecological and environmental factors affecting mosquito populations.

This comprehensive approach stands on the shoulders of pioneering work by Allen [34], Capocelli and Ricciardi [40], Turelli [41], and Nipa et al. [42], who, in various contexts, enriched our understanding of population dynamics by incorporating stochastic elements and exploring diverse growth rate equations. Their collective contributions underscore the versatility of these modeling frameworks, transcending disciplinary boundaries and fostering a more profound comprehension of the intricate interplay of factors influencing population growth.

For clarity, the parameters used in our model are summarized in the Table 1.

Table 1

Table 1. Parameters of the anopheles mosquito population dynamics model.

2.4.2 Dynamics of infected and uninfected mosquito populations

We categorize the mosquito population into two groups: X₁(t), representing those infected with Microsporidia MB, and X₂(t), comprising the uninfected or wild mosquitoes. The total population is X(t) = X₁(t)+X₂(t).

The dynamics of these populations are described by a stochastic differential equation with Lévy noise:

$\begin{array}{l}d{X}_k\left(t\right)=\left({\alpha }_k\left(1-\frac{X\left(t\right)}{K}\right)-{\lambda }_k\right)X_k\left(t\right)dt+{\sigma }_k{X}_k\left(t\right)d{W}_k\left(t\right)\\ +{\zeta }_k{X}_k\left(t\right){\displaystyle {\int }_{ℝ}}y{Ñ}_k\left(dt,dy\right)& (10)\end{array}$

This model captures the growth α_k, mortality λ_k, and stochastic influences within each mosquito population, providing a foundation for optimizing symbiont-infected mosquito release strategies for malaria control.

2.4.3 Interaction between infected and uninfected mosquitoes

This section explores the interactions between two key mosquito populations: those infected with Microsporidia MB, a symbiont known for inhibiting malaria transmission, and the uninfected ones. Understanding the dynamics and potential coexistence of these populations is crucial for developing effective malaria control strategies. Drawing from the works of Rafikov et al. [28] and Wyse et al. [43], we model these dynamics using a stochastic differential equation with Lévy noise:

$\begin{array}{l}\frac{d{X}_k\left(t\right)}{X_k\left(t^{-}\right)}=\left[{\displaystyle \sum _{j=1}^n}{\alpha }_k{\gamma }_{k,j}{\theta }_j\left(t\right)\left(1-\frac{{\displaystyle {\sum }_{k=1}^n}{X}_k\left(t\right)}{K}\right)-{\lambda }_k\right]dt\\ +{\sigma }_{k}d{W}_k\left(t\right)+{\zeta }_k{\displaystyle {\int }_{ℝ}}y{Ñ}_k\left(dt,dy\right),& (11)\end{array}$

where

$\begin{array}{ll}{\theta }_k\left(t\right)=\frac{X_k\left(t\right)}{{\displaystyle {\sum }_{k=1}^n}{X}_k\left(t\right)}& (12)\end{array}$

signifies the proportion of class k mosquitoes within the total population.

Focusing on a specific case with n = 2 classes (infected and uninfected), we define the total mosquito population as X(t) = X₁(t)+X₂(t). The fractions of infected and uninfected mosquitoes are represented as

$\begin{array}{ll}{\theta }_1\left(t\right)=\frac{X_1\left(t\right)}{X_1\left(t\right)+X_2\left(t\right)},& (13)\end{array}$

and

$\begin{array}{ll}{\theta }_2(t)=\frac{X_2(t)}{X_1(t)+X_2(t)}.& (14)\end{array}$

Accordingly, Equation 11 can be reformulated for this two-class scenario as follows:

$\begin{array}{l}\frac{d{X}_k(t)}{X_k(t^{-})}=\left[{\displaystyle \sum _{j=1}^n}{\gamma }_{k,j}{\theta }_j(t)\beta (t)-{\lambda }_k\right]dt+{\sigma }_{k}d{W}_k(t)\\ +{\zeta }_k{\displaystyle {\int }_{ℝ}}y{Ñ}_k(dt,dy).& (15)\end{array}$

The parameters in this equation are defined as follows:

• γ_k,j: Interaction coefficients

– γ₁₁: Infected-Infected Mating Rate

– γ₁₂: Infected-Uninfected Mating Rate

– γ₂₁: Uninfected-Infected Mating Rate

– γ₂₂: Uninfected-Uninfected Mating Rate

• θ_j: Proportional rates

– θ₁: Infection prevalence ratio

– θ₂: Uninfected ratio

• λ_k: Mortality rate of the mosquito population

• σ_k: Volatility term for random population fluctuations

• dW_k(t): Wiener process for standard Brownian motion

• Ñ_k(dt, dy): Lévy measure for impulsive behavior in the population

In the following, we consider that the diffusion term of σ₂ and the scaling constant for the intensity of Lévy jumps ζ₂ for the wild (uninfected) mosquitoes are null, and that the corresponding intensity of Lévy jumps for infected mosquitoes ζ₁ is the unit.

2.4.4 Release control problem formulation

In the context of this study, the performance function is formulated as follows:

$\begin{array}{ll}{\mathcal{H}}_f^{R,{\theta }_1}(s,x_1,x_2)={\text{E}}^{s,x_1,x_2}\left[{\displaystyle {\int }_0^{{\tau }_S}}{e}^{-\alpha (s+t)}g(X(t),u(t))dt\right]& (16)\end{array}$

Here, α > 0 and γ ∈ (0, 1) are constants, and ${\text{E}}^{x_1,x_2}$ denotes the expectation with respect to the probability law ${\mathbb{P}}^{x_1,x_2}$ of the state vector (X₁, X₂) given initial conditions $(X_1(0^{-}),X_2(0^{-}))=(x_1,x_2)$. The function $g:\mathcal{S}\times \mathcal{A}\to ℝ$ is given by $\frac{R^{\gamma }}{\gamma }$, where $\mathcal{S}$ is a fixed domain in ℝ and $\mathcal{A}$ is a given subset of ℝ².

Let us consider u(t) = (R(t), θ₁(t)) ∈ [0, ∞) × [0, 1] as the control applied at time t, chosen from a set of admissible controls $\mathcal{U}$. Furthermore, biological constraints impose that the total mosquito population X(t) = X₁(t) + X₂(t) is non-negative for all t ≥ 0.

The objective is to identify an optimal release control strategy, defined as follows:

$\begin{array}{ll}\begin{array}{l}\text{Find }\Psi (s,x_1,x_2)\text{and }u\in \mathcal{U}\text{such that}\\ \text{Ψ}(s,x_1,x_2)={\displaystyle \underset{u\in \mathcal{U}}{sup}}{\mathcal{H}}_f^u(s,x_1,x_2)={\mathcal{H}}_f^{u^{*}}(s,x_1,x_2).\end{array}& (17)\end{array}$

This is subject to the dynamic equation:

$\begin{array}{l}d{X}^u(t)=((m_1(\theta )\beta (t)-m_2(\lambda ))X^u(t)-R(t))dt\\ +{\sigma }_1{\theta }_1{X}^u(t)dW(t)+{\zeta }_1{\theta }_1{X}^u(t){\displaystyle {\int }_{ℝ}}y{Ñ}_k(dt,dy)& (18)\end{array}$

Here, $X(0^{-})=x=x_1+x_2>0$, and β(t) describes a mean-reverting process. The functions m₁(θ) and m₂(λ) are defined as follows:

$\begin{array}{ll}{m}_1(\theta )={\gamma }_{11}{\theta }_1^2+({\gamma }_{12}+{\gamma }_{21}){\theta }_1{\theta }_2+{\gamma }_{22}{\theta }_2^2,& (19)\end{array}$

$\begin{array}{ll}{m}_2(\lambda )={\lambda }_1{\theta }_1+{\lambda }_2{\theta }_2.& (20)\end{array}$

2.4.4.1 Stopping time τ_S

The stopping time τ_S is defined as the first time at which the mosquito population reaches zero:

$\begin{array}{ll}{\tau }_S=inf\left\{t\ge 0:X(t)=0\right\}.& (21)\end{array}$

Assumption 2.7. The generator of the Lévy diffusion process $Z(t)=\left[\begin{array}{l}s+t\\ X(t)\end{array}\right]$ for t ≥ 0 and initial condition $Z(0^{-})=\left[\begin{array}{l}s\\ x\end{array}\right]$ is given by:

$\begin{array}{ll}\begin{array}{l}{A}^u\psi =\frac{\partial \psi }{\partial s}+\Theta (x)\frac{\partial \psi }{\partial x}+\frac{1}{2}{\Sigma }^2(x)\frac{{\partial }^2\psi }{\partial x^2}\\ +{\displaystyle {\int }_{ℝ}}\left\{\psi (s,x(1+{\theta }_{1}y))-\psi (s,x)-{\theta }_{1}xy\frac{\partial \psi }{\partial x}(s,x)\right\}\vartheta (dy).\end{array}& (22)\end{array}$

For simplicity and to avoid compromising the theoretical integrity of the control problem, we assume that the derivative of ψ with respect to β is not directly captured in the generator. Consequently, we use the simplified generator as shown above.

In this context, β(t) is a stochastic function representing the mosquito growth rate, and X(t) represents the population of infected and uninfected mosquitoes. While β(t) was originally a deterministic function, making it stochastic adds realism to the model. However, we choose to simplify the problem by not directly including β(t) in the generator. This is justified because β(t) influences the system implicitly through its effect on X(t), and its direct stochastic influence can be captured indirectly.

By focusing on the primary stochastic process X(t) and simplifying the model, we make the problem more tractable for analysis and control without significantly compromising accuracy. This assumption allows us to balance model complexity with practical applicability, especially when empirical evidence supports that the effects of β(t) can be captured indirectly through X(t).

3 Mathematical analysis

3.1 Analytical results

In this section, we rigorously analyze the solution to the stochastic optimal control problem via dynamic programming, a methodology crucial for optimizing the Hamiltonian–Jacobi–Bellman (HJB) equation. The HJB equation serves as the fundamental equation governing the release control problem, underpinning our analytical framework.

Proposition 3.1. The solution to the dynamic (9) representing the total mosquito population is given by:

$\begin{array}{l}X(t)=X_0.\exp (({\displaystyle {\int }_{\left|y\right|<R}\left(ln(1+\zeta y)-\zeta y\right)\vartheta (dy)}\\ -\left(\lambda +\frac{1}{2}{\sigma }^2\right))t+\sigma W(t))\times \exp \left({\Lambda }_1(t)+{\Lambda }_2(t)\right),& (23)\end{array}$

with :

$\begin{array}{ll}{\Lambda }_1(t)={\displaystyle {\int }_0^t}\beta (s)ds,& (24)\end{array}$

$\begin{array}{ll}{\Lambda }_2(t)={\displaystyle {\int }_0^t}{\displaystyle {\int }_{ℝ}}ln(1+\zeta y)Ñ(ds,dy).& (25)\end{array}$

In particular, if −1 < ζy ≤ 0,

$\begin{array}{r}X\left(t\right)=X_0.exp\left({\Lambda }_1\left(t\right)-\left(\lambda +\frac{1}{2}{\sigma }^2+\zeta {\displaystyle {\int }_{ℝ}}y\vartheta \left(dy\right)\right)t\right)\\ \times exp\left(\sigma W\left(t\right)+{\Lambda }_2\left(t\right)\right).& (26)\end{array}$

Proof. Consider Z(t) = h(t, X(t)) such that h(t, X(t)) = ln(X(t)) is an Ito-Lévy process. Using the Ito-Lévy formula, we get,

$\begin{array}{l}dZ(t)=\frac{\partial h}{\partial t}(t,X(t))dt+\frac{\partial h}{\partial x}(t,X(t))[\left(X(t)·\beta (t)-\lambda X(t)\right)dt\\ +\sigma X(t)dW(t)]+\frac{1}{2}{\sigma }^{2}X{(t)}^2\frac{{\partial }^{2}h}{{\partial }^2{x}^2}(t,X(t))dt\\ +{\displaystyle {\int }_{\left|y\right|<R}\{h(t,X(t^{-})+\zeta yX(t^{-}))-h(t,X(t^{-}))}\\ -\frac{\partial h}{\partial x}(t,X(t))\zeta yX(t^{-})\}\vartheta (dy)dt+{\displaystyle {\int }_{ℝ}(h(t,X(t^{-})}\\ +\zeta yX(t^{-}))-h(t,X(t^{-})))\tilde{N}(dt,dy)\\ =\frac{\partial h}{\partial t}(t,X(t))dt+\frac{\partial h}{\partial x}(t,X(t))\left[\left(\beta (t)-\lambda \right)dt+\sigma dW(t)\right]X(t)\\ +\frac{1}{2}{\sigma }^{2}X{(t)}^2\frac{{\partial }^{2}h}{{\partial }^2{x}^2}(t,X(t))dt\\ +{\displaystyle {\int }_{\left|y\right|<R}\{h(t,(1+\zeta y)X(t^{-}))-h(t,X(t^{-}))}\\ -\frac{\partial h}{\partial x}(t,X(t))\zeta yX(t^{-})\}\vartheta (dy)dt\\ +{\displaystyle {\int }_{ℝ}\left(h(t,(1+\zeta y)X(t^{-}))-h(t,X(t^{-}))\right)\tilde{N}(dt,dy)}.\end{array}$

By replacing h by ln(x), dZ(t) is expressed as follows:

$\begin{array}{l}dZ(t)=0.dt+\frac{1}{X(t)}\left[\left(\beta (t)-\lambda \right)dt+\sigma dW(t)\right]X(t)\\ +\frac{1}{2}{\sigma }^{2}X{(t)}^2\frac{1}{X^2(t)}dt+{\displaystyle {\int }_{\left|y\right|<R}\{ln((1+\zeta y)X(t^{-}))}\\ -ln(X(t^{-}))-\frac{1}{X(t)}\zeta yX(t^{-})\}\vartheta (dy)dt\\ +{\displaystyle {\int }_{ℝ}\left(ln((1+\zeta y)X(t^{-}))-ln(X(t^{-}))\right)\tilde{N}(dt,dy)}\\ =\left(\beta (t)-(\lambda +\frac{1}{2}{\sigma }^2)\right)dt+\sigma dW(t)\\ +{\displaystyle {\int }_{\left|y\right|<R}\left\{ln(1+\zeta y)-\zeta y\right\}\vartheta (dy)dt}\\ +{\displaystyle {\int }_{ℝ}\left(ln(1+\zeta y)\right)\tilde{N}(dt,dy)}.\end{array}$

Thus, by integrating both sides of the previous equality, we have :

$\begin{array}{l}Z(t)=Z(0)+{\displaystyle {\int }_0^t\left(\beta (s)-(\lambda +\frac{1}{2}{\sigma }^2)\right)}ds+\sigma W(t)\\ +{\displaystyle {\int }_0^t{\displaystyle {\int }_{\left|y\right|<R}\left(ln(1+\zeta y)-\zeta y\right)\vartheta (dy)ds}}\\ +{\displaystyle {\int }_0^t{\displaystyle {\int }_{ℝ}\left(ln(1+\zeta y)\right)\tilde{N}(ds,dy)}}\\ =Z(0)+{\displaystyle {\int }_0^t\beta }(s)ds-\left(\lambda +\frac{1}{2}{\sigma }^2\right)t+\sigma W(t)\\ +\left({\displaystyle {\int }_{\left|y\right|<R}\left(ln(1+\zeta y)-\zeta y\right)\vartheta (dy)}\right)t\\ +{\displaystyle {\int }_0^t{\displaystyle {\int }_{ℝ}\left(ln(1+\zeta y)\right)\tilde{N}(ds,dy)}}\\ =Z(0)+(\left({\displaystyle {\int }_{\left|y\right|<R}\left(ln(1+\zeta y)-\zeta y\right)\vartheta (dy)}\right)\\ -\left(\lambda +\frac{1}{2}{\sigma }^2\right))t+\sigma W(t)+{\displaystyle {\int }_0^t\beta }(s)ds\\ +{\displaystyle {\int }_0^t{\displaystyle {\int }_{ℝ}\left(ln(1+\zeta y)\right)\tilde{N}(ds,dy)}}.\end{array}$

Taking the exponential of both sides facilitates deriving the solution presented in Equation 23. For further details, the reader is referred to the comprehensive discussions in Øksendal and Sulem [39] and Nunno et al. [44].

Proposition 3.2. A solution to the release control problem (17) is given by :

$\begin{array}{ll}\{\begin{array}{l}{R}^{*}={(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}x\\ {\theta }_1^{*}=\frac{(({\gamma }_{12}+{\gamma }_{21}){\theta }_2)\beta -{\lambda }_1}{{\sigma }_1^2(1-\gamma )-2{\gamma }_{11}\beta }-{\displaystyle {\int }_{ℝ}}\frac{y\left[1-{(1+{\theta }_{1}y)}^{\gamma -1}\right]}{{\sigma }_1^2(1-\gamma )-2{\gamma }_{11}\beta }\vartheta (dy),\end{array}& (27)\end{array}$

where $\mathcal{C}$ is a constant determining the corresponding fraction over the total mosquito required for the release. An explicit formula for the constant is given in the following proposition.

Proof. Consider the following controlled process

$\begin{array}{l}Z(t)=\left[\begin{array}{l}s+t\\ X(t)\\ \beta (t)\end{array}\right];t\ge 0,Z(0^{-})=\left[\begin{array}{l}s\\ x\\ \beta \end{array}\right].\end{array}$

The corresponding generator to Z(t) denoted by 𝔸^u is given by:

$\begin{array}{l}{\text{A}}^u\psi (z)=\frac{\partial \psi }{\partial s}+\left((m_1(\theta )\beta -m_2(\lambda ))x-R\right)\frac{\partial \psi }{\partial x}+a(b-\beta )\frac{\partial \psi }{\partial \beta }\\ +\frac{1}{2}{\sigma }_1^2{\theta }_1^2{x}^2\frac{{\partial }^2\psi }{\partial x^2}+\frac{1}{2}{c}^2\beta \frac{{\partial }^2\psi }{\partial {\beta }^2}+{\sigma }_1{\theta }_{1}c\sqrt{\beta }\frac{{\partial }^2\psi }{\partial x\partial \beta }\rho (x,\beta )\\ +{\displaystyle {\int }_{ℝ}\{\psi (s,x(1+{\theta }_{1}y),\beta )-\psi (s,x,\beta )}\\ -{\theta }_{1}xy\frac{\partial \psi }{\partial x}(s,x,\beta )\}\vartheta (dy).\end{array}$

Choose

$\begin{array}{l}\psi (z)=\psi (s,x,\beta )=e^{-\alpha s}\varphi (x,\beta ).\end{array}$

Based on Assumption 2.7, we have

$\begin{array}{l}{\text{A}}^u\psi (z)=e^{-\alpha s}{\text{A}}_0^u\varphi (x),\end{array}$

with

$\begin{array}{l}{\mathbb{A}}_0^u\varphi (x)=-\alpha \varphi (x)+\left((m_1(\theta )\beta -m_2(\lambda ))x-R\right){\varphi }^{\prime }(x)\\ +\frac{1}{2}{\sigma }_1^2{\theta }_1^2{x}^2{\varphi }^{″}(x)+{\displaystyle {\int }_{ℝ}\{\varphi (x(1+{\theta }_{1}y))}\\ -\varphi (x)-{\theta }_{1}xy{\varphi }^{\prime }(x)\}\vartheta (dy).\end{array}$

Setting $\varphi (x)=\mathcal{C}{x}^{\gamma }$ we get,

$\begin{array}{l}\begin{array}{l}{\text{A}}_0^u\varphi (x)+g(x,u)=-\alpha \mathcal{C}{x}^{\gamma }\\ +((m_1(\theta )\beta -m_2(\lambda ))x-R)\mathcal{C}\gamma x^{\gamma -1}\\ +\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2{x}^2\gamma (\gamma -1)x^{\gamma -2}+\\ \mathcal{C}{x}^{\gamma }{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy)+\frac{R^{\gamma }}{\gamma }.\end{array}\end{array}$

Let us define the function H(R, θ₁) as follows:

$\begin{array}{ll}\begin{array}{l}H(R,{\theta }_1)=-\alpha \mathcal{C}{x}^{\gamma }+((m_1(\theta )\beta -m_2(\lambda ))x-R)\mathcal{C}\gamma x^{\gamma -1}\\ +\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2{x}^2\gamma (\gamma -1)x^{\gamma -2}\\ +\mathcal{C}{x}^{\gamma }{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy)+\frac{R^{\gamma }}{\gamma }.\end{array}& (28)\end{array}$

For any release value R_v ≥ R and prevalence ratio θ_v,1 ≥ θ₁ simultaneously,

$\begin{array}{l}H(R_v,{\theta }_{v,1})\ge H(R,{\theta }_1).\end{array}$

This implies that H(R, θ₁) is a non-decreasing function with respect to R and θ₁. Furthermore, H is concave in (R, θ₁) and attains its maximum, that is:

$\begin{array}{ll}\frac{\partial H}{\partial R}=-\mathcal{C}\gamma x^{\gamma -1}+R^{\gamma -1}=0,& (29)\end{array}$

$\begin{array}{l}\frac{\partial H}{\partial {\theta }_1}=\mathcal{C}\gamma x^{\gamma -1}((2{\gamma }_{11}{\theta }_1+({\gamma }_{12}+{\gamma }_{21}){\theta }_2)\beta -{\lambda }_1)x\\ +\mathcal{C}{\sigma }_1^2{\theta }_1{x}^2\gamma (\gamma -1)x^{\gamma -2}\\ +\mathcal{C}{x}^{\gamma }{\displaystyle {\int }_{ℝ}}\left[\gamma y{(1+{\theta }_{1}y)}^{\gamma -1}-\gamma y\right]\vartheta (dy)=0.& (30)\end{array}$

The relation (29) implies that,

$\begin{array}{ll}R={(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}x.& (31)\end{array}$

Moreover, from (30) we have

$\begin{array}{l}0=\mathcal{C}\gamma x^{\gamma }\left((2{\gamma }_{11}{\theta }_1+({\gamma }_{12}+{\gamma }_{21}){\theta }_2)\beta -{\lambda }_1\right)+\mathcal{C}{\sigma }_1^2{\theta }_1\gamma (\gamma -1)x^{\gamma }\\ +\mathcal{C}{x}^{\gamma }{\displaystyle {\int }_{ℝ}\left[\gamma y{(1+{\theta }_{1}y)}^{\gamma -1}-\gamma y\right]}\vartheta (dy).\\ 0=\mathcal{C}\gamma x^{\gamma }[(2{\gamma }_{11}{\theta }_1+({\gamma }_{12}+{\gamma }_{21}){\theta }_2)\beta -{\lambda }_1+{\sigma }_1^2{\theta }_1(\gamma -1)\\ +{\displaystyle {\int }_{ℝ}\left[y{(1+{\theta }_{1}y)}^{\gamma -1}-y\right]}\vartheta (dy)].\\ 0=[(2{\gamma }_{11}\beta +{\sigma }_1^2(\gamma -1)){\theta }_1+(({\gamma }_{12}+{\gamma }_{21}){\theta }_2)\beta -{\lambda }_1\\ +{\displaystyle {\int }_{ℝ}\left[y{(1+{\theta }_{1}y)}^{\gamma -1}-y\right]}\vartheta (dy)].\end{array}$

Thus

$\begin{array}{r}({\sigma }_1^2(1-\gamma )-2{\gamma }_{11}\beta ){\theta }_1=(({\gamma }_{12}+{\gamma }_{21}){\theta }_2)\beta -{\lambda }_1\\ +{\displaystyle {\int }_{ℝ}}\left[y{(1+{\theta }_{1}y)}^{\gamma -1}-y\right]\vartheta (dy).\end{array}$

Accordingly,

$\begin{array}{l}{\theta }_1=\frac{(({\gamma }_{12}+{\gamma }_{21}){\theta }_2)\beta -{\lambda }_1}{({\sigma }_1^2(1-\gamma )-2{\gamma }_{11}\beta )}\\ +\frac{1}{({\sigma }_1^2(1-\gamma )-2{\gamma }_{11}\beta )}{\displaystyle {\int }_{ℝ}}\left[y{(1+{\theta }_{1}y)}^{\gamma -1}-y\right]\vartheta (dy)\\ =\frac{(({\gamma }_{12}+{\gamma }_{21}){\theta }_2)\beta -{\lambda }_1}{({\sigma }_1^2(1-\gamma )-2{\gamma }_{11}\beta )}-{\displaystyle {\int }_{ℝ}}\frac{y\left[1-{(1+{\theta }_{1}y)}^{\gamma -1}\right]}{({\sigma }_1^2(1-\gamma )-2{\gamma }_{11}\beta )}\vartheta (dy).\end{array}$

From Equations 29, 30, and based on the HJB criterion for jumps diffusion in Theorem (2.6), the optimal release R^* and symbiont transmission rate ${\theta }_1^{*}$ can be expressed as follows:

$\begin{array}{ll}{R}^{*}={(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}x,& (32)\end{array}$

$\begin{array}{ll}{\theta }_1^{*}=\frac{(({\gamma }_{12}+{\gamma }_{21}){\theta }_2)\beta -{\lambda }_1}{{\sigma }_1^2(1-\gamma )-2{\gamma }_{11}\beta }-{\displaystyle {\int }_{ℝ}}\frac{y\left[1-{(1+{\theta }_{1}y)}^{\gamma -1}\right]}{{\sigma }_1^2(1-\gamma )-2{\gamma }_{11}\beta }\vartheta (dy).& (33)\end{array}$

Proposition 3.3. The constant $\mathcal{C}$ in Equation 32 is given by:

$\begin{array}{ll}{(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}=\{\begin{array}{ll}\frac{1}{1-\gamma }\left[\alpha +\frac{1}{2}{\sigma }_1^2{\theta }_1^2(1-\gamma )\gamma -\gamma (m_1(\theta )\beta -m_2(\lambda ))\right],& \text{ if }\vartheta =0,\\ \frac{1}{1-\gamma }\left[\alpha +\frac{1}{2}{\sigma }_1^2{\theta }_1^2(1-\gamma )\gamma -\gamma (m_1(\theta )\beta -m_2(\lambda ))-{\displaystyle {\int }_{ℝ}}f(y;{\theta }_1)\vartheta (dy)\right],& \text{ if }\vartheta \ne 0.\end{array}& (34)\end{array}$

where

$\begin{array}{l}{\displaystyle {\int }_{ℝ}}f(y;{\theta }_1)\vartheta (dy)={\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy).\end{array}$

Remark 3.4. This constant, $\mathcal{C}$, can be denoted as $\mathcal{C}(T)$ because it depends on temperature T. This dependence arises from the fact that the mosquito death rates λ₁ and λ₂ are functions of temperature. Consequently, this makes the release strategy $R^{*}={(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}x$ temperature dependent, even though it can inherit the dependency from x, which is already influenced by environmental conditions. This ensures that R is proportional to the total mosquito population at any time, and that the proportionality factor reflects the impact of temperature on the release.

As depicted in Figure 2, the release factor ${(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}$ quantifies the necessary adjustments in mosquito release quantities throughout the year 2023. This model reflects how local temperature fluctuations critically influence mosquito population management strategies. The graph highlights the relationship between temperature and the scaling of mosquito release, underlining the importance of temperature in optimizing biological control measures against mosquito-borne diseases. This adaptive strategy aids in tailoring interventions to be more effective and responsive to changing climatic conditions, enhancing the efficacy of disease control efforts.

Figure 2

Figure 2. Temperature-Dependent Release Factor from January 1, 2023, to January 1, 2024: This plot displays the computed release factor ${(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}$ as a function of daily mean temperatures from Niger, Ghana, Benin and Kenya. Each data point illustrates the factor adjustment required to modulate daily mosquito releases in response to environmental temperature variations.

Proof. [Proposition 3.3] The release strategy $R^{*}={(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}x$ of infected mosquitoes by the symbiont microbe imposes that ${\text{A}}_0^u\psi (x)+g(x,u)=0$. Therefore,

$\begin{array}{l}0=-\alpha \mathcal{C}{x}^{\gamma }+((m_1(\theta )\beta -m_2(\lambda ))x-{(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}x)\mathcal{C}\gamma x^{\gamma -1}\\ +\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2{x}^2\gamma (\gamma -1)x^{\gamma -2}+\mathcal{C}{x}^{\gamma }{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy)\\ +\frac{{({(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}x)}^{\gamma }}{\gamma }=-\alpha \mathcal{C}{x}^{\gamma }+((m_1(\theta )\beta -m_2(\lambda ))-{(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}})\mathcal{C}\gamma x^{\gamma }\\ +\frac{1}{2}{\sigma }_1^2{\theta }_1^2\mathcal{C}\gamma (\gamma -1)x^{\gamma }+\mathcal{C}{x}^{\gamma }{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy)\\ +\frac{{(\mathcal{C}\gamma )}^{\frac{\gamma }{\gamma -1}}}{\gamma }{x}^{\gamma }.\end{array}$

Since the total mosquito population is assumed to be strictly positive (x > 0), we have:

$\begin{array}{l}0=-\alpha \mathcal{C}+((m_1(\theta )\beta -m_2(\lambda ))-{(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}})\mathcal{C}\gamma \\ +\frac{1}{2}{\sigma }_1^2{\theta }_1^2\mathcal{C}\gamma (\gamma -1)+\mathcal{C}{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy)+\frac{{(\mathcal{C}\gamma )}^{\frac{\gamma }{\gamma -1}}}{\gamma }\\ =-\alpha +((m_1(\theta )\beta -m_2(\lambda ))-{(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}})\gamma +\frac{1}{2}{\sigma }_1^2{\theta }_1^2\gamma (\gamma -1)\\ +{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy)+\frac{{(\mathcal{C})}^{\frac{1}{\gamma -1}}{\gamma }^{\frac{\gamma }{\gamma -1}}}{\gamma }\\ =-\alpha +((m_1(\theta )\beta -m_2(\lambda )))\gamma +\frac{1}{2}{\sigma }_1^2{\theta }_1^2\gamma (\gamma -1)\\ +{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy)+(1-\gamma ){(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}.\end{array}$

Thus

$\begin{array}{l}(1-\gamma ){(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}=\alpha -\gamma ((m_1(\theta )\beta -m_2(\lambda )))+\frac{1}{2}{\sigma }_1^2{\theta }_1^2\gamma (1-\gamma )\\ -{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy).\end{array}$

Accordingly, if ϑ = 0,

$\begin{array}{l}{\left(\mathcal{C}\gamma \right)}^{\frac{1}{\gamma -1}}=\frac{1}{\left(1-\gamma \right)}(\alpha -\gamma \left((m_1(\theta )\beta -m_2(\lambda ))\right)\\ +\frac{1}{2}{\sigma }_1^2{\theta }_1^2\gamma (1-\gamma )).& (35)\end{array}$

Now let us consider that ϑ ≠ 0. If θ₁ ≠ 0, define the function f(y; θ₁) as:

$f(y;{\theta }_1)={(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y$

such that f(y; θ₁) ≠ 0 for all y ∈ ℝ. Therefore,

${\displaystyle {\int }_{ℝ}}f(y;{\theta }_1)\vartheta (dy)\ne 0.$

Hence, the following equation holds:

$\begin{array}{l}{\left(\mathcal{C}\gamma \right)}^{\frac{1}{\gamma -1}}=\frac{1}{1-\gamma }[\alpha +\frac{1}{2}{\sigma }_1^2{\theta }_1^2(1-\gamma )\gamma -\gamma (m_1(\theta )\beta -m_2(\lambda ))\\ -{\displaystyle {\int }_{ℝ}f}(y;{\theta }_1)\vartheta (dy)].& (36)\end{array}$

Proposition 3.5. Given that θ₁ and R are solutions to the release problem, the following inequality holds:

$\begin{array}{l}|{\text{A}}^u\psi (t,x;\beta )|\le e^{-\alpha t}{x}^{\gamma }[{\kappa }_0+\varsigma (\beta )+\mathcal{C}{\displaystyle {\int }_{ℝ}|}{(1+{\theta }_{1}y)}^{\gamma }-1\\ -\gamma {\theta }_{1}y|\vartheta (dy)],& (37)\end{array}$

where ς(β) is defined as:

$\begin{array}{ll}\varsigma (\beta )={\kappa }_1{e}^{-at}+{\kappa }_2{\displaystyle {\int }_0^t}{e}^{-a(t-s)}\sqrt{\beta (s)}dV(s),& (38)\end{array}$

and κ₀, κ₁, κ₂ are constants.

Moreover, if there exists a finite value,

$\begin{array}{ll}\Delta >{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy)& (39)\end{array}$

such that:

$\begin{array}{ll}\text{E}\left[{\displaystyle {\int }_0^T}{e}^{-2\alpha t}{X}^{2\gamma }(t)dt\right]=\Delta -{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy),& (40)\end{array}$

then the expected absolute value of ψ at the exit time, along with the accumulated generator and gradient terms over the interval [0, τ_S], is finite.

Proof.

$\begin{array}{l}{\text{A}}_0^u\varphi (x)+g(x,u)=-\alpha \mathcal{C}{x}^{\gamma }+((m_1(\theta )\beta -m_2(\lambda ))x-R)\mathcal{C}\gamma x^{\gamma -1}\\ +\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2{x}^2\gamma (\gamma -1)x^{\gamma -2}\\ +\mathcal{C}{x}^{\gamma }{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy)+\frac{R^{\gamma }}{\gamma }.\end{array}$

Therefore, we have

$\begin{array}{l}{\text{A}}_0^u\varphi (x)=-\alpha \mathcal{C}{x}^{\gamma }+((m_1(\theta )\beta -m_2(\lambda ))x-R)\mathcal{C}\gamma x^{\gamma -1}\\ +\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2{x}^2\gamma (\gamma -1)x^{\gamma -2}\\ +\mathcal{C}{x}^{\gamma }{\displaystyle {\int }_{ℝ}}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}\vartheta (dy)\end{array}$

$\begin{array}{l}=x^{\gamma }[-\alpha \mathcal{C}+(m_1(\theta )\beta -m_2(\lambda ))-{\left(\mathcal{C}\gamma \right)}^{\frac{1}{\gamma -1}}\mathcal{C}\gamma +\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2\gamma (\gamma -1)\\ +\mathcal{C}{\displaystyle {\int }_{ℝ}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}}\vartheta (dy)]\\ =x^{\gamma }[(m_1(\theta )\beta -m_2(\lambda ))-\left({\left(\mathcal{C}\gamma \right)}^{\frac{\gamma }{\gamma -1}}+\alpha \mathcal{C}\right)+\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2\gamma (\gamma -1)\\ +\mathcal{C}{\displaystyle {\int }_{ℝ}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}}\vartheta (dy)].\end{array}$

$\begin{array}{l}{\text{A}}_0^u\varphi (x)=x^{\gamma }[m_1(\theta )\beta -({\left(\mathcal{C}\gamma \right)}^{\frac{\gamma }{\gamma -1}}+\alpha \mathcal{C}+\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2\gamma (1-\gamma )\\ +m_2(\lambda ))\\ +\mathcal{C}{\displaystyle {\int }_{ℝ}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}}\vartheta (dy)]\\ =x^{\gamma }[\left({\gamma }_{11}{\theta }_1^2+({\gamma }_{12}+{\gamma }_{21}){\theta }_1{\theta }_2+{\gamma }_{22}{\theta }_2^2\right)\beta \\ -\left({\left(\mathcal{C}\gamma \right)}^{\frac{\gamma }{\gamma -1}}+\alpha \mathcal{C}+\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2\gamma (1-\gamma )+m_2(\lambda )\right)\\ +\mathcal{C}{\displaystyle {\int }_{ℝ}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}}\vartheta (dy)]\\ (obtainedbyreplacing{m}_1(\theta )).\end{array}$

$\begin{array}{l}|{\text{A}}_0^u\varphi (x)|\le x^{\gamma }|[\delta \left(b+\left({\beta }_0-b\right)e^{-at}+c{\displaystyle {\int }_0^t{e}^{-a(t-s)}}\sqrt{\beta (s)}dV(s)\right)\\ +\mathcal{C}{\displaystyle {\int }_{ℝ}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}}\vartheta (dy)]|\\ \le x^{\gamma }[\delta \left(b+{\beta }_0{e}^{-at}+c{\displaystyle {\int }_0^t{e}^{-a(t-s)}}\sqrt{\beta (s)}dV(s)\right)\\ +|\mathcal{C}{\displaystyle {\int }_{ℝ}\left\{{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y\right\}}\vartheta (dy)|]\\ \le x^{\gamma }[{\kappa }_0+{\kappa }_1{e}^{-at}+{\kappa }_2{\displaystyle {\int }_0^t{e}^{-a(t-s)}}\sqrt{\beta (s)}dV(s)\\ +\mathcal{C}{\displaystyle {\int }_{ℝ}|}{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy)].\end{array}$

since ${\text{A}}^u\psi (z)=e^{-\alpha t}{\text{A}}_0^u\varphi (x)$, we obtain the following :

$\begin{array}{ll}|{\text{A}}^u\psi (t,x)|=|e^{-\alpha t}{\text{A}}_0^u\varphi (x)|& (41)\end{array}$

$\begin{array}{ll}=e^{-\alpha t}|{\text{A}}_0^u\varphi (x)|& (42)\end{array}$

$\begin{array}{ll}\le e^{-\alpha t}{x}^{\gamma }[{\kappa }_0+{\kappa }_1{e}^{-at}+{\kappa }_2{\displaystyle {\int }_0^t{e}^{-a(t-s)}}\sqrt{\beta (s)}dV(s)& (43)\end{array}$

$\begin{array}{ll}+\mathcal{C}{\displaystyle {\int }_{ℝ}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|}\vartheta (dy)].& (44)\end{array}$

By applying Fubini's theorem, we interchange the expectation and integral, justified by the non-negativity of the integrand:

$\begin{array}{r}\text{E}\left[{\displaystyle {\int }_0^{\tau }}{e}^{-2\alpha t}{X}^{2\gamma }(t)dt\right]=\text{E}\left[{\displaystyle {\int }_0^{\infty }}{e}^{-2\alpha t}{X}^{2\gamma }(t)1_{\left\{t\le \tau \right\}}dt\right]\\ ={\displaystyle {\int }_0^{\infty }}\text{E}\left[e^{-2\alpha t}{X}^2(t)1_{\left\{t\le \tau \right\}}\right]dt.\end{array}$

Since $\psi (t)=\mathcal{C}{e}^{-\alpha t}{X}^{\gamma }(t)$, we have :

$\begin{array}{l}\text{E}\left[{\mathcal{C}}^2{e}^{-2\alpha t}{X}^{2\gamma }(t)\right]=\text{E}\left[\psi {(t)}^2\right].\end{array}$

Thus

$\begin{array}{l}{\displaystyle {\int }_0^{\infty }}\text{E}\left[e^{-2\alpha t}{X}^{2\gamma }(t)1_{\left\{t\le \tau \right\}}\right]dt=\frac{1}{{\mathcal{C}}^2}{\displaystyle {\int }_0^{\infty }}\text{E}\left[\psi {(t)}^2{1}_{\left\{t\le \tau \right\}}\right]dt\\ =\frac{1}{{\mathcal{C}}^2}{\displaystyle {\int }_0^{{\tau }_S}}\text{E}\left[\psi {(t)}^2\right]dt\\ =\Delta -{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy).\end{array}$

We get

$\begin{array}{ll}{\displaystyle {\int }_0^{{\tau }_S}}\text{E}\left[\psi {(t)}^2\right]dt={\mathcal{C}}^2(\Delta -{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy)).& (45)\end{array}$

Since ψ(t)² ≥ 0, we have $\left|\psi (t)\right|=\sqrt{\psi {(t)}^2}$. Applying Jensen's inequality for the concave function $\sqrt{x}$, we obtain:

$\begin{array}{l}\text{E}\left[|\psi (t)|1_{\left\{t\le \tau \right\}}\right]\le \sqrt{\text{E}\left[\psi {(t)}^2{1}_{\left\{t\le \tau \right\}}\right]}.\end{array}$

Therefore,

$\begin{array}{l}{\displaystyle {\int }_0^{\infty }}\text{E}\left[|\psi (t)|1_{\left\{t\le \tau \right\}}\right]dt\le {\displaystyle {\int }_0^{\infty }}\sqrt{\text{E}\left[\psi {(t)}^2{1}_{\left\{t\le \tau \right\}}\right]}dt\\ \le \sqrt{{\tau }_S}.\sqrt{{\displaystyle {\int }_0^{{\tau }_S}}\text{E}\left[\psi {(t)}^2\right]dt}\text{(Cauchy-Schwarz inequality)}\\ \le \sqrt{{\tau }_S}.\mathcal{C}\sqrt{(\Delta -{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy))}.\end{array}$

Considering the relation (44), we note that

$\begin{array}{l}{\displaystyle {\int }_0^{{\tau }_S}}\text{E}|{\text{A}}^u\psi (X(t))|dt\le \left[{\kappa }_0+{\kappa }_2+\mathcal{C}{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy)\right]\text{E}\left[{\displaystyle {\int }_0^{\tau }}{e}^{-\alpha t}{X}^{\gamma }(t)dt\right]\\ \le \left[\kappa +\mathcal{C}{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy)\right]{\displaystyle {\int }_0^{\tau }}\text{E}\left[e^{-\alpha t}{X}^{\gamma }(t)\right]dt\\ \le \left[\kappa +\mathcal{C}{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy)\right].\sqrt{{\tau }_S}.\sqrt{\text{E}\left[{\displaystyle {\int }_0^{\tau }}{e}^{-2\alpha t}{X}^{2\gamma }(t)dt\right]}.\end{array}$

Hence,

$\begin{array}{l}{\displaystyle {\int }_0^{{\tau }_S}}\text{E}\left|\frac{{\text{A}}^u\psi (X(t))}{\left[{\kappa }_{\mathcal{C}}+{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy)\right]}\right|dt\\ \le \sqrt{{\tau }_S}.\mathcal{C}\sqrt{(\Delta -{\displaystyle {\int }_{ℝ}}|{(1+{\theta }_{1}y)}^{\gamma }-1-\gamma {\theta }_{1}y|\vartheta (dy))}& (46)\end{array}$

where ${\kappa }_{\mathcal{C}}$ is a positive constant.

Remark 3.6. Since ${\theta }_1^{*}+{\theta }_2^{*}=1$, the expression (33) can be written as :

$\begin{array}{l}{\theta }_1^{*}(t)=\frac{({\gamma }_{12}+{\gamma }_{21})\beta -{\lambda }_1}{{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta }\\ -{\displaystyle {\int }_{ℝ}}\frac{y\left[1-{(1+{\theta }_{1}y)}^{\gamma -1}\right]}{{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta }\vartheta (dy).& (47)\end{array}$

Theorem 3.7. Suppose that for each mosquito population $x\in \mathcal{S}$, there exists a release strategy v such that v: = û(x), satisfying the condition:

$\begin{array}{l}{\text{A}}^v\psi (x)+G(x,v)=0.\end{array}$

and that u^*(t) = û(X(t)). Then $u^{*}(t)=(R^{*}(t),{\theta }_1^{*}(t))$ is an optimal control given by the following:

$\begin{array}{ll}{R}^{*}(t)={(\mathcal{C}\gamma )}^{\frac{1}{\gamma -1}}x(t),& (48)\end{array}$

$\begin{array}{l}{\theta }_1^{*}(t)=\frac{({\gamma }_{12}+{\gamma }_{21})\beta (t)-{\lambda }_1}{{\sigma }^2(1-\gamma )-\left[2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21})\right]\beta (t)}\\ -{\displaystyle {\int }_{ℝ}}\frac{y\left[1-{(1+{\theta }_{1}y)}^{\gamma -1}\right]}{{\sigma }^2(1-\gamma )-\left[2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21})\right]\beta (t)}\vartheta (dy).& (49)\end{array}$

The value function for the release control problem described by (17) is then expressed by

$\begin{array}{ll}\text{Ψ}(s,x_1,x_2)=\psi (s,x;\beta )=e^{-\alpha s}\mathcal{C}{x}^{\gamma }(s),& (50)\end{array}$

where $\mathcal{C}$ is a constant determined by (34), which differentiates cases where the measure ϑ(dy) is zero and non-zero. The term involving the integral over ℝ signifies contributions from jumps, with ϑ(dy) representing a measure.

Proof. Let $u\in \mathcal{U}$ and for k = 1, 2, …, set τ_k = k∧τ_S.

Since Proposition 3.5 ensures that

$\begin{array}{l}\text{E}\left[|\psi (X(\tau ))|+{\displaystyle {\int }_0^{{\tau }_S}}|{\text{A}}^u\psi (X(t))|dt\right]<\infty ,\end{array}$

we can apply Dynkin's Formula II to obtain

$\begin{array}{l}{\text{E}}^x\left[\psi (X({\tau }_k))\right]=\psi (x)+{\text{E}}^x\left[{\displaystyle {\int }_0^{{\tau }_k}}{\text{A}}^u\psi (X(t))dt\right]\\ \le \psi (x)-\text{E}\left[{\displaystyle {\int }_0^{{\tau }_k}}G(X(t),u(t))dt\right].\end{array}$

Therefore,

$\begin{array}{l}\psi (x)\ge {\text{E}}^x\left[\psi (X({\tau }_k))\right]+\text{E}\left[{\displaystyle {\int }_0^{{\tau }_k}}G(X(t),u(t))dt\right]\\ \ge {\text{E}}^x\left[\psi (X({\tau }_k))+{\displaystyle {\int }_0^{{\tau }_k}}G(X(t),u(t))dt\right]\\ \ge {\displaystyle \underset{k\to \infty }{lim inf}}{\text{E}}^x\left[\psi (X({\tau }_k))+{\displaystyle \underset{0}{\overset{{\tau }_k}{\int }}}G(X(t),u(t))dt\right]\\ \ge {\text{E}}^{x_1,x_2,\beta }\left[{\displaystyle {\int }_0^{\infty }}G(X(t),u(t))dt\right].\end{array}$

Thus, for $u\in \mathcal{U}$,

$\begin{array}{l}\psi (x)\ge {\mathcal{H}}_f^u(s,x_1,x_2).\end{array}$

Accordingly,

$\begin{array}{ll}\psi (x)\ge {\displaystyle \underset{u\in \mathcal{U}}{sup}}{\mathcal{H}}_f^u(s,x_1,x_2)=\text{Ψ}(x),\forall x\in \mathcal{S}.& (51)\end{array}$

Furthermore, since u^* = û(x) then,

$\begin{array}{ll}\psi (x)={\text{E}}^{x_1,x_2,\beta }\left[{\displaystyle {\int }_0^{\infty }}G(X(t),v(t))dt\right]& (52)\end{array}$

$\begin{array}{ll}\le {\displaystyle \underset{u\in \mathcal{U}}{sup}}{\mathcal{H}}_f^u(s,x_1,x_2,\beta )=\text{Ψ}(x),\forall x\in \mathcal{S}.& (53)\end{array}$

Therefore

$\begin{array}{l}\psi (s,x)=\text{Ψ}(s,x)\\ =\text{Ψ}(s,x;\beta )\\ =e^{-\alpha s}\mathcal{C}{x}^{\gamma }(s).\end{array}$

Using Propositions 3.2 and 3.3, the proof is completed.

Proposition 3.8. Assuming that ${\theta }_1^{*}$ and ${\theta }_2^{*}$ satisfy the conditions of the release control problem and that γ₁₁ and γ₂₂ are distinct from $\frac{1}{2}({\gamma }_{12}+{\gamma }_{21})$, then the integral

$\begin{array}{ll}{\mathcal{I}}_{\beta }={\displaystyle {\int }_{ℝ}}\frac{y\left[1-{(1+{\theta }_{1}y)}^{\gamma -1}\right]}{{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta }\vartheta (dy),& (54)\end{array}$

can take the form :

$\begin{array}{ll}{\mathcal{I}}_{\beta }=\frac{({\gamma }_{12}+{\gamma }_{21})\beta -{\lambda }_1}{{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta }-\frac{2{\gamma }_{22}\beta -{\lambda }_2}{(2{\gamma }_{22}-({\gamma }_{12}+{\gamma }_{21}))\beta }.& (55)\end{array}$

Consequently, the infection prevalence ratio during release can be articulated as follows :

$\begin{array}{ll}{\theta }_1^{*}(t)=\{\begin{array}{ll}\frac{({\gamma }_{12}+{\gamma }_{21})\beta (t)-{\lambda }_1}{{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta (t)}& if\vartheta =0\\ \frac{2{\gamma }_{22}\beta (t)-{\lambda }_2}{(2{\gamma }_{22}-({\gamma }_{12}+{\gamma }_{21}))\beta (t)}& if\vartheta \ne 0\end{array}& (56)\end{array}$

Note that the integral ${\mathcal{I}}_{\beta }$ is independent of θ₁ and θ₂.

Proof. Suppose that $\frac{\partial H}{\partial {\theta }_2}=0,$ where H is given by

$\begin{array}{l}H=-\alpha \mathcal{C}{x}^{\gamma }+(({\gamma }_{11}{\theta }_1^2+({\gamma }_{12}+{\gamma }_{21}){\theta }_1{\theta }_2+{\gamma }_{22}{\theta }_2^2)\beta \\ -({\lambda }_1{\theta }_1+{\lambda }_2{\theta }_2))x-R)\mathcal{C}\gamma x^{\gamma -1}\\ +\mathcal{C}\frac{1}{2}{\sigma }_1^2{\theta }_1^2{x}^2\gamma (\gamma -1)x^{\gamma -2}+\mathcal{C}{x}^{\gamma }{\displaystyle {\int }_{ℝ}\{{(1+{\theta }_{1}y)}^{\gamma }}\\ -1-\gamma {\theta }_{1}y\}\vartheta (dy)+\frac{R^{\gamma }}{\gamma }.\end{array}$

Thus,

$\begin{array}{ll}((({\gamma }_{12}+{\gamma }_{21}){\theta }_1+2{\gamma }_{22}{\theta }_2)\beta -{\lambda }_2)\mathcal{C}\gamma x^{\gamma -1}=0.& (57)\end{array}$

Since ${\theta }_1^{*}=1-{\theta }_2^{*}$, Equation 57 becomes:

$\begin{array}{ll}((({\gamma }_{12}+{\gamma }_{21})(1-{\theta }_2^{*})+2{\gamma }_{22}{\theta }_2^{*})\beta -{\lambda }_2)\mathcal{C}\gamma x^{\gamma -1}=0.& (58)\end{array}$

Therefore, we obtain:

$\begin{array}{ll}(({\gamma }_{12}+{\gamma }_{21})\beta +(2{\gamma }_{22}-({\gamma }_{12}+{\gamma }_{21}))\beta {\theta }_2^{*}-{\lambda }_2)\mathcal{C}\gamma x^{\gamma -1}=0.& (59)\end{array}$

Given that the total mosquito population X(t) and the constants $\mathcal{C}$ and γ are non-zero, it follows that:

$\begin{array}{ll}{\theta }_2^{*}=\frac{{\lambda }_2-({\gamma }_{12}+{\gamma }_{21})\beta }{(2{\gamma }_{22}-({\gamma }_{12}+{\gamma }_{21}))\beta }.& (60)\end{array}$

This expression for ${\theta }_2^{*}(t)$ is derived by rearranging the terms and solving for ${\theta }_2^{*}$ under the condition that the denominator is non-zero, which should be verified in the context of the problem. Therefore

$\begin{array}{ll}1-{\theta }_2^{*}(t)=1-\frac{{\lambda }_2-({\gamma }_{12}+{\gamma }_{21})\beta (t)}{(2{\gamma }_{22}-({\gamma }_{12}+{\gamma }_{21}))\beta (t)}& (61)\end{array}$

$\begin{array}{ll}=\frac{2{\gamma }_{22}\beta (t)-{\lambda }_2}{(2{\gamma }_{22}-({\gamma }_{12}+{\gamma }_{21}))\beta (t)}& (62)\end{array}$

$\begin{array}{ll}={\theta }_1^{*}(t)& (63)\end{array}$

$\begin{array}{ll}=\frac{({\gamma }_{12}+{\gamma }_{21})\beta (t)-{\lambda }_1}{{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta (t)}-{\mathcal{I}}_{\beta }.& (64)\end{array}$

Hence,

$\begin{array}{ll}{\mathcal{I}}_{\beta }=\frac{({\gamma }_{12}+{\gamma }_{21})\beta -{\lambda }_1}{{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta }-\frac{2{\gamma }_{22}\beta -{\lambda }_2}{(2{\gamma }_{22}-({\gamma }_{12}+{\gamma }_{21}))\beta }.& (65)\end{array}$

Furthermore the integral ${\mathcal{I}}_{\beta }$ is independent of θ₁ and θ₂. This completes the derivation of θ₁ and θ₂, as required by the release control problem.

Proposition 3.9. (Sensitivity of the infection prevalence ratio) Let θ₁ and θ₂ the fractions of infected and uninfected mosquitoes, such that :

$\begin{array}{l}{\theta }_1^{*}(t)=\frac{({\gamma }_{12}+{\gamma }_{21})\beta -{\lambda }_1}{{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta }\\ -{\displaystyle {\int }_{ℝ}}\frac{y\left[1-{(1+{\theta }_{1}y)}^{\gamma -1}\right]}{{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta }\vartheta (dy)& (66)\end{array}$

Therefore, as the growth rate β increases the infection prevalence ratio θ₁ increases. That is, for all β ∈ (0, 1)

$\begin{array}{ll}\frac{\partial {\theta }_1}{\partial \beta }>0.& (67)\end{array}$

Proof. Assume that γ₁₁ and γ₂₂ are different from $\frac{1}{2}({\gamma }_{12}+{\gamma }_{21})$.

(i) Considering the case where the Lévy measure ϑ = 0, by setting θ₁ = f(β) such that for $\beta \in (0,1)-\left\{\frac{{\sigma }^2(1-\gamma )}{2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21})}\right\}$ the function is well defined,

Thus, if ${\sigma }^2>\frac{{\lambda }_1(({\gamma }_{11}-{\gamma }_{12})+({\gamma }_{11}-{\gamma }_{21}))}{(1-\gamma )({\gamma }_{12}+{\gamma }_{21})}$

$\begin{array}{l}{f}^{\prime }(\beta )=\frac{|\begin{array}{ll}{\gamma }_{12}+{\gamma }_{21}& -{\lambda }_1\\ -(({\gamma }_{12}+{\gamma }_{21})-2{\gamma }_{11})& {\sigma }^2(1-\gamma )\end{array}|}{{\left[{\sigma }^2(1-\gamma )-(2{\gamma }_{11}-({\gamma }_{12}+{\gamma }_{21}))\beta \right]}^2}>0.\end{array}$

(ii) In the same way for ϑ ≠ 0, f′(β) > 0 whenever ${\gamma }_{22}\ne \frac{1}{2}({\gamma }_{12}+{\gamma }_{21})$.

3.2 Numerical results

3.2.1 Numerical simulation of mosquito population dynamics

This section provides a comprehensive description of the numerical methods used to simulate the dynamics of mosquito populations influenced by stochastic factors and vector control strategies, incorporating Lévy-driven stochastic differential equations (SDEs). The simulations are performed using Python, leveraging its robust libraries for numerical computations and stochastic modeling.

Algorithm 1

Algorithm 1. Numerical simulation of mosquito population dynamics using Lévy-driven SDEs - Part 1

Algorithm 2

Algorithm 2. Numerical simulation of mosquito population dynamics using Lévy-driven SDEs - Part 2

3.3 Numerical simulations of the optimal control problem

In our study, we employ numerical simulations to investigate the intricacies of stochastic optimal control in mosquito release strategies, focusing on the effects of Lévy measures and seasonality. While we do not present the cases with zero Lévy noise and non-zero Lévy noise without seasonality, they serve as baseline references to understand the infection dynamics under more controlled conditions. The first simulation, depicted in Figure 3, incorporates seasonal dynamics under a zero Lévy measure. This simulation demonstrates how infection prevalence and the ratio of uninfected individuals fluctuate with the seasons, while also highlighting the necessary adjustments in stochastic growth rates to optimize the release process. It emphasizes the critical role of seasonality in the success of mosquito control efforts, where timing and environmental conditions are pivotal to strategic planning. In Figure 4, we extend the analysis by introducing non-zero Lévy measures, in combination with seasonality, to explore their joint impact on mosquito release strategies. This simulation provides a multifaceted view of how Lévy-driven fluctuations and seasonal variations influence infection dynamics and control strategies. The results underscore the complexity of designing effective mosquito release strategies that must account for both ecological variability and stochastic uncertainties.

Figure 3

Figure 3. Stochastic optimal control and seasonal dynamics in mosquito release: a zero Lévy measure simulation.

Figure 4

Figure 4. Seasonality and non-zero Lévy measure in stochastic control of mosquito release.

3.4 Numerical illustration of the infection prevalence ratio θ₁

This section presents both 2D and 3D analyses to deepen our understanding of how symbiont prevalence in mosquito populations correlates with key environmental and biological factors. Specifically, the 2D analyses (Figures 5A, B) explore the relationships between symbiont prevalence, temperature, and the stochastic growth rate β. These correlations provide insights into the individual effects of temperature and growth rates on symbiont spread within mosquito populations. Extending beyond this, the 3D visualizations (Figures 5C, D) synthesize these variables, offering a comprehensive view that highlights the intricate interplay between environmental conditions and the efficacy of mosquito release strategies. Through this multifaceted analytical approach, we gain a clearer understanding of the factors that influence symbiont prevalence, aiding in the optimization of mosquito release methodologies.

Figure 5

Figure 5. Interplay of environmental conditions and symbiont prevalence in mosquito releases: a 2D and 3D analytical perspective. (A) θ₁_values with Levy measure ϑ = 0 in 2D. (B) θ₁_values with Levy measure ϑ ≠ 0 in 2D. (C) θ₁_values with ϑ = 0 in 3D. (D) θ₁_values with ϑ ≠ 0 in 3D.

3.4.1 Spatial representations based on temperature

This part of the paper presents numerical simulations and spatial illustrations that complement our analytical findings, while also considering the biological and ecological implications of our study. To ensure the model's compatibility with biological settings, we have adjusted the infection prevalence ratio using a logistic transformation. Specifically, we normalize the symbiont prevalence rate during release with the logistic function $\frac{1}{1+exp(-y)}$. This guarantees that the ratio remains non-negative and is bounded between 0 and 1, which is essential for maintaining the practical applicability and interpretability of our model in real-world biological scenarios. The proportion θ₁ represents the infection prevalence ratio, influenced by temperature and varying growth rates over time. Conversely, θ₂ denotes the uninfected ratio, reflecting the susceptibility of wild mosquitoes to the Plasmodium parasite. This becomes critical in estimating the maximum malaria prevalence during periods when the symbiont prevalence drops or becomes inefficient, possibly due to environmental changes. The simulations were conducted with both zero and non-zero Lévy measures to demonstrate the varying effects under different conditions. To define the infection prevalence ratio quantitatively, we derive the formula for θ₁ as follows:

$\begin{array}{ll}{\theta }_1^{*}\left(\beta \left(t\right),T\right)=\{\begin{array}{ll}\frac{1}{1+exp\left(-\left(\frac{\left({\gamma }_{12}+{\gamma }_{21}\right)\beta \left(t\right)-\frac{1}{J{T}^2+KT+L}}{{\sigma }^2\left(1-\gamma \right)-\left(2{\gamma }_{11}-\left({\gamma }_{12}+{\gamma }_{21}\right)\right)\beta \left(t\right)}\right)\right)}& if\vartheta =0\\ \frac{1}{1+exp\left(-\left(\frac{2{\gamma }_{22}\beta \left(t\right)-\frac{1}{J{T}^2+KT+L}}{\left(2{\gamma }_{22}-\left({\gamma }_{12}+{\gamma }_{21}\right)\right)\beta \left(t\right)}\right)\right)}& if\vartheta \ne 0.\end{array}& (68)\end{array}$

To map the malaria transmission potential risk area, we consider the uninfected ratio ${\theta }_2^{*}$, which can be represented by:

$\begin{array}{ll}{\theta }_2^{*}(\beta (t),T)=1-{\theta }_1^{*}(\beta (t),T)& (69)\end{array}$

3.5 Spatial illustrations

In this subsection, we present spatial illustrations depicting the areas of malaria transmission potential risk based on the variability in symbiont prevalence over a year of release. This risk is closely linked to the effectiveness of the release strategy.

The maps show variations in malaria transmission risks across different study areas, illustrating the effects of mosquito releases in Benin, Niger, Ghana, and Kenya on reducing malaria transmission risks. At day 28 of the releases, the lowest transmission risks observed are 82% for Benin (Figure 6A), 76% for Niger (Figure 7A), 77% for Ghana (Figure 8A), and 76% for Kenya (Figure 9A). By day 252, there is a notable decrease in malaria transmission risks: 44% for Benin (Figure 6B), 37% for Niger (Figure 7B), 37% for Ghana (Figure 8B), and 37% for Kenya (Figure 9B). This demonstrates the impact of the release strategy on malaria transmission dynamics, primarily influenced by vector population densities. The green areas indicate low transmission risk (Figures 6–9), while the red areas indicate high transmission risk, highlighting regions where the release strategy has less impact.

Figure 6

Figure 6. Mapping malaria transmission risks in Benin: evaluating mosquito release strategy efficiency. (A) 28 days. (B) 252 days. This map illustrates malaria transmission risk areas in Benin, influenced by temperature variations, geographical location, and the timing and quantity of mosquito releases.

Figure 7

Figure 7. Mapping malaria transmission risks in Niger: evaluating mosquito release strategy efficiency. (A) 28 days. (B) 252 days. This map shows malaria transmission risk areas, impacted by temperature differences, geographical location, and the timing and quantity of mosquito releases.

Figure 8

Figure 8. Mapping malaria transmission risks in Ghana: evaluating mosquito release strategy efficiency. (A) 28 days. (B) 252 days. This map depicts malaria transmission risk areas, affected by temperature variations, geographical location, and the timing and quantity of mosquito releases.

Figure 9

Figure 9. Mapping malaria transmission risks in Kenya: evaluating mosquito release strategy efficiency. (A) 28 days. (B) 252 days. This map highlights malaria transmission risk areas, shaped by temperature differences, geographical location, and the timing and quantity of mosquito releases.

3.6 Validation

The integrity of our model was rigorously tested using data on the prevalence of the microbe Microsporidia MB collected in Kenya. Figure 10 illustrates a direct correlation between symbiont prevalence in mosquito populations and malaria transmission risk across Kenya. The map highlights regions based on confirmed prevalence and overlays these with malaria transmission risks, providing a compelling visual validation of our model. Our analysis mapped Microsporidia MB prevalence data against areas of significant malaria transmission potential, identifying specific zones of concern. The validation revealed that regions with higher symbiont prevalence exhibit lower malaria transmission risks, underscoring the importance of symbiont activity in mitigating malaria transmission. In conclusion, our model's validation highlights the necessity for tailored strategies that address malaria transmission challenges, including the innovative use of symbionts for vector control.

Figure 10

Figure 10. Correlating Microsporidia MB prevalence with low malaria transmission risk area: a validation of the model. This figure shows how the prevalence of Microsporidia MB in mosquito populations reduces malaria transmission risk in various regions of Kenya, validating our model's effectiveness.

3.6.1 Model parameter values

The sensitivity of the stochastic growth rate β and the population of mosquitoes X(t) under the parameters in Tables 2, 3 are in Table 4 in Appendix.

Table 2

Table 2. Parameters for mosquito population simulation.

Table 3

Table 3. Parameters for CIR process with jumps.

4 Discussion

This study investigated the use of stochastic optimal control to manage the release of Anopheles mosquitoes infected with Microsporidia MB, aiming to reduce malaria transmission. By employing Hamilton–Jacobi–Bellman (HJB) Equation 17, we identified the optimal daily release and symbiont prevalence necessary for effective malaria control. Our results demonstrate that an increase in the growth rate of Anopheles mosquitoes enhances the infection prevalence ratio (Proposition 3.9) and reduces the proportion of uninfected mosquitoes, thereby amplifying the potential for disease control.

A key improvement in this study is the refined calculation of the infection prevalence ratio using a sigmoid function, providing a more accurate representation of biological and ecological realities. This ratio, influenced by stochastic factors such as time and temperature, highlights the variability in infection rates within mosquito populations. The relationship between the optimal infection prevalence (θ₁) and the uninfected mosquito ratio (θ₂) illustrates malaria transmission risk dynamics, as shown in Figures 6–9. Our spatial analysis across these regions reveals a general reduction in malaria transmission following the release of infected mosquitoes. However, this reduction is not uniform; it varies with time, temperature, and geographical conditions, with fluctuations reflecting the adaptation period of released mosquitoes into wild populations.

Recent advances in malaria control emphasize the ongoing need for innovative strategies. For example, the CDC recently reported the effectiveness of vector control measures such as pyrethroid-chlorfenapyr-impregnated nets, despite emerging challenges like insecticide resistance [46]. Concurrently, mathematical modeling has been pivotal in evaluating various interventions. Notably, one recent study integrated machine learning with transmission dynamics to optimize malaria intervention strategies [47]. Similarly, Ozodiegwu et al. developed a model incorporating drug resistance and preventive measures to assess malaria dynamics in Nigeria [48]. In this context, our study introduces a novel approach, leveraging stochastic optimal control and mosquitoes infected with Microsporidia MB to disrupt Plasmodium transmission. This strategy complements existing interventions, reinforcing flexibility and adaptation to local environmental conditions. Moreover, by integrating stochastic elements such as white and Lévy noise, we have enhanced the model's capacity to account for diverse environmental variabilities, offering a more comprehensive depiction of mosquito population dynamics.

The implications of our findings are significant for public health, particularly in malaria-endemic regions of West and East Africa. The adaptability of the infected mosquitoes and the variability in symbiont prevalence, influenced by environmental and ecological factors, underscore the complexity of implementing effective vector control strategies. Additionally, challenges related to potential ecological impacts, logistical concerns, and public acceptance must be carefully considered to ensure successful implementation of this approach. Ongoing monitoring will be essential to assess the ecological impacts and efficacy of the releases over time. While this study does not explicitly address the economic costs of mosquito releases, future research should incorporate economic modeling to evaluate the cost-effectiveness of various release strategies. This would aid in optimizing the timing, frequency, and scale of releases, thereby maximizing their impact while ensuring financial viability for large-scale implementation.

Future studies should also explore additional factors that influence the success of Microsporidia MB-infected mosquito releases, such as the interaction between symbionts and insecticide-resistant mosquito populations, potential ecological side effects, and long-term sustainability. Incorporating more extensive field data, especially from regions with varied climates, could further refine the model and increase its applicability across different malaria-endemic zones.

In conclusion, this research advances malaria vector control strategies by employing stochastic optimal control under a Lévy process to optimize the release of Microsporidia MB-infected mosquitoes. Our model demonstrates the significant influence of symbiont activity on mosquito population dynamics, affecting both infection prevalence and malaria transmission risks. By accounting for environmental factors such as seasonality and incorporating Lévy noise, the model captures a more realistic and dynamic picture of mosquito behavior under varying conditions. Spatial analysis across Ghana, Niger, Benin, and Kenya reveals the gradual impact of infected mosquitoes on local populations, supporting the potential for integrated approaches in malaria control. These findings underscore the importance of adaptable, region-specific strategies, with further research needed to explore economic feasibility and the interplay between biological and ecological factors for long-term success.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Ethics statement

The manuscript presents research on animals that do not require ethical approval for their study.

Author contributions

SA: Conceptualization, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing. HT: Conceptualization, Supervision, Validation, Visualization, Writing – review & editing. PN: Conceptualization, Methodology, Supervision, Writing – review & editing. BK: Conceptualization, Investigation, Methodology, Supervision, Visualization, Writing – review & editing. SA: Formal analysis, Investigation, Supervision, Visualization, Writing – review & editing. JH: Conceptualization, Funding acquisition, Resources, Supervision, Validation, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study, conducted by the International Centre of Insect Physiology and Ecology, is part of the Symbiont Modeling and Deployment Strategies project, generously supported by the Bill and Melinda Gates Foundation (INV-022584). The opinions, findings, conclusions, or recommendations expressed in this publication are those of the authors and do not necessarily reflect the view of the donor.

Acknowledgments

The authors express their gratitude to the Microsporidia team for their support in this work.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2024.1465153/full#supplementary-material

References

1. CDC. Comparison of the Plasmodium Species Which Cause Human Malaria. (2013). Available at: http://www.cdc.gov/dpdx/malaria/index.html (accessed November 17, 2023).

Google Scholar

2. Marshall JM, Taylor CE. Malaria control with transgenic mosquitoes. PLoS Med. (2009) 6:e1000020. doi: 10.1371/journal.pmed.1000020

PubMed Abstract | Crossref Full Text | Google Scholar

3. World Health Organization. Fact Sheet Malaria. WHO Media Centre (2017). Available at: https://apps.who.int/iris/bitstream/handle/10665/204183/Fact_Sheet_WHD_2014_EN_1632.pdf (accessed November 19, 2023).

Google Scholar

4. World Health Organization. World Malaria Report 2021. World Health Organization (2021). Licence: CC BY-NC-SA 3.0 IGO. Available at: https://www.who.int/teams/global-malaria-programme/reports/world-malaria-report-2021 (accessed November 19, 2023).

Google Scholar

5. Mwamtobe PM, Abelman S, Tchuenche JM, Kasambara A. Optimal (control of) intervention strategies for malaria epidemic in Karonga District, Malawi. In: Abstract and Applied Analysis. Hindawi Publishing Corporation (2014), p. 1–20. doi: 10.1155/2014/594256

Crossref Full Text | Google Scholar

6. World Health Organization. World Malaria Report: 20 years of global progress and challenges. World Health Organization (2020). Available at: https://www.who.int/publications/i/item/9789240015791 (accessed November 19, 2023).

Google Scholar

7. Gomes FM, Barillas-Mury C. Infection of anopheline mosquitoes with Wolbachia: implications for malaria control. PLoS Pathog. (2018) 14:e1007333. doi: 10.1371/journal.ppat.1007333

PubMed Abstract | Crossref Full Text | Google Scholar

8. Hughes GL, Koga R, Xue P, Fukatsu T, Rasgon JL. Wolbachia infections are virulent and inhibit the human malaria parasite Plasmodium falciparum in Anopheles gambiae. PLoS Pathog. (2011) 7:e1002043. doi: 10.1371/journal.ppat.1002043

PubMed Abstract | Crossref Full Text | Google Scholar

9. Herren JK, Mbaisi L, Mararo E, Makhulu EE, Mobegi VA, Butungi H, et al. A microsporidian impairs Plasmodium falciparum transmission in Anopheles arabiensis mosquitoes. Nat Commun. (2020) 11:2187. doi: 10.1038/s41467-020-16121-y

PubMed Abstract | Crossref Full Text | Google Scholar

10. Akorli J, Akorli EA, Tetteh SNA, Amlalo GK, Opoku M, Pwalia R, et al. Microsporidia MB is found predominantly associated with Anopheles gambiae s.s. and Anopheles coluzzii in Ghana. Sci Rep. (2021) 11:18658. doi: 10.1038/s41598-021-98268-2

PubMed Abstract | Crossref Full Text | Google Scholar

11. Moustapha LM, Sadou IM, Arzika II, Maman LI, Gomgnimbou MK, Konkobo M, et al. First identification of microsporidia MB in Anopheles coluzzii from Zinder City, Niger. Parasit Vectors. (2024) 17:39. doi: 10.1186/s13071-023-06059-7

PubMed Abstract | Crossref Full Text | Google Scholar

12. Ahouandjinou MJ, Sovi A, Sidick A, Sewadé W, Koukpo CZ, Chitou S, et al. First report of natural infection of Anopheles gambiae s.s. and Anopheles coluzzii by Wolbachia and Microsporidia in Benin: a cross-sectional study. Malaria J. (2024) 23:72. doi: 10.1186/s12936-024-04906-1

PubMed Abstract | Crossref Full Text | Google Scholar

13. Fox RM, Weiser J. A microsporidian parasite of Anopheles gambiae in Liberia. J Parasitol. (1959) 45:21–30. doi: 10.2307/3274782

PubMed Abstract | Crossref Full Text | Google Scholar

14. Reynolds D. Laboratory studies of the microsporidian Plistophora culicis (Weiser) infecting Culex pipiens fatigans Wied. Bull Entomol Res. (1970) 60:339–49. doi: 10.1017/S0007485300040852

PubMed Abstract | Crossref Full Text | Google Scholar

15. Schenker W, Maier W, Seitz H. The effects of Nosema algerae on the development of Plasmodium yoelii nigeriensis in Anopheles stephensi. Parasitol Res. (1992) 78:56–9. doi: 10.1007/BF00936182

PubMed Abstract | Crossref Full Text | Google Scholar

16. Bano L. Partial inhibitory effect of Plistophora culicis on the sporogonic cycle of Plasmodium cynomolgi in Anopheles stephensi. Nature. (1958) 181:430–430. doi: 10.1038/181430a0

PubMed Abstract | Crossref Full Text | Google Scholar

17. Huang M, Tang M, Yu J. Wolbachia infection dynamics by reaction-diffusion equations. Sci China Math. (2015) 58:77–96. doi: 10.1007/s11425-014-4934-8

PubMed Abstract | Crossref Full Text | Google Scholar

18. Matsufuji T, Seirin-Lee S. The optimal strategy of incompatible insect technique (IIT) using Wolbachia and the application to malaria control. J Theor Biol. (2023) 569:111519. doi: 10.1016/j.jtbi.2023.111519

PubMed Abstract | Crossref Full Text | Google Scholar

19. Zhang X, Tang S, Cheke RA. Models to assess how best to replace dengue virus vectors with Wolbachia-infected mosquito populations. Math Biosci. (2015) 269:164–77. doi: 10.1016/j.mbs.2015.09.004

PubMed Abstract | Crossref Full Text | Google Scholar

20. Zheng B, Tang M, Yu J. Modeling Wolbachia spread in mosquitoes through delay differential equations. SIAM J Appl Math. (2014) 74:743–70. doi: 10.1137/13093354X

Crossref Full Text | Google Scholar

21. Cai L, Bao L, Rose L, Summers J, Ding W. Malaria modeling and optimal control using sterile insect technique and insecticide-treated net. Appl Anal. (2022) 101:1715–34. doi: 10.1080/00036811.2021.1999419

Crossref Full Text | Google Scholar

22. Fister KR, McCarthy ML, Oppenheimer SF. Diffusing wild type and sterile mosquitoes in an optimal control setting. Math Biosci. (2018) 302:100–15. doi: 10.1016/j.mbs.2018.05.015

PubMed Abstract | Crossref Full Text | Google Scholar

23. Multerer L, Smith T, Chitnis N. Modeling the impact of sterile males on an Aedes aegypti population with optimal control. Math Biosci. (2019) 311:91–102. doi: 10.1016/j.mbs.2019.03.003

PubMed Abstract | Crossref Full Text | Google Scholar

24. Fister KR, McCarthy ML, Oppenheimer SF, Collins C. Optimal control of insects through sterile insect release and habitat modification. Math Biosci. (2013) 244:201–12. doi: 10.1016/j.mbs.2013.05.008

PubMed Abstract | Crossref Full Text | Google Scholar

25. Stone CM. Transient population dynamics of mosquitoes during sterile male releases: modelling mating behaviour and perturbations of life history parameters. PLoS ONE. (2013) 8:e76228. doi: 10.1371/journal.pone.0076228

PubMed Abstract | Crossref Full Text | Google Scholar

26. Thomé RC, Yang HM, Esteva L. Optimal control of Aedes aegypti mosquitoes by the sterile insect technique and insecticide. Math Biosci. (2010) 223:12–23. doi: 10.1016/j.mbs.2009.08.009

PubMed Abstract | Crossref Full Text | Google Scholar

27. Diaz H, Ramirez A, Olarte A, Clavijo C. A model for the control of malaria using genetically modified vectors. J Theor Biol. (2011) 276:57–66. doi: 10.1016/j.jtbi.2011.01.053

PubMed Abstract | Crossref Full Text | Google Scholar

28. Rafikov M, Bevilacqua L, Wyse A. Optimal control strategy of malaria vector using genetically modified mosquitoes. J Theor Biol. (2009) 258:418–25. doi: 10.1016/j.jtbi.2008.08.006

PubMed Abstract | Crossref Full Text | Google Scholar

29. Meetei MZ, Zafar S, Zaagan AA, Mahnashi AM, Idrees M. Dengue transmission dynamics: a fractional-order approach with compartmental modeling. Fractal Fract. (2024) 8:207. doi: 10.3390/fractalfract8040207

Crossref Full Text | Google Scholar

30. Alyobi S, Jan R. Qualitative and quantitative analysis of fractional dynamics of infectious diseases with control measures. Fractal Fract. (2023) 7:400. doi: 10.3390/fractalfract7050400

Crossref Full Text | Google Scholar

31. Tang TQ, Jan R, Khurshaid A, Shah Z, Vrinceanu N, Racheriu M. Analysis of the dynamics of a vector-borne infection with the effect of imperfect vaccination from a fractional perspective. Sci Rep. (2023) 13:14398. doi: 10.1038/s41598-023-41440-7

PubMed Abstract | Crossref Full Text | Google Scholar

32. Jan R, Boulaaras S, Alyobi S, Rajagopal K, Jawad M. Fractional dynamics of the transmission phenomena of dengue infection with vaccination. Discrete Contin Dyn Syst Ser S. (2022) 16:2096–117. doi: 10.3934/dcdss.2022154

PubMed Abstract | Crossref Full Text | Google Scholar

33. Jan A, Jan R, Khan H, Zobaer M, Shah R. Fractional-order dynamics of Rift Valley fever in ruminant host with vaccination. Commun Math Biol Neurosci. (2020) 2020. doi: 10.28919/cmbn/5017

Crossref Full Text | Google Scholar

34. Allen E. Environmental variability and mean-reverting processes. Discrete Contin Dyn Syst Ser B. (2016) 21:2073–89. doi: 10.3934/dcdsb.2016037

Crossref Full Text | Google Scholar

35. Allen EJ, Allen LJ, Schurz H. A comparison of persistence-time estimation for discrete and continuous stochastic population models that include demographic and environmental variability. Math Biosci. (2005) 196:14–38. doi: 10.1016/j.mbs.2005.03.010

PubMed Abstract | Crossref Full Text | Google Scholar

36. Luquin-Covarrubias MA. Morales-Bojórquez E. Effects of stochastic growth on population dynamics and management quantities estimated from an integrated catch-at-length assessment model: Panopea globosa as case study. Ecol Modell. (2021) 440:109384. doi: 10.1016/j.ecolmodel.2020.109384

Crossref Full Text | Google Scholar

37. Ayoubi T, Bao H. Persistence and extinction in stochastic delay Logistic equation by incorporating Ornstein-Uhlenbeck process. Appl Math Comput. (2020) 386:125465. doi: 10.1016/j.amc.2020.125465

Crossref Full Text | Google Scholar

38. Chidzalo P, Ngare PO, Mung'atu JK. Pricing weather derivatives under a tri-variate stochastic model. Sci Afr. (2023) 21:e01768. doi: 10.1016/j.sciaf.2023.e01768

Crossref Full Text | Google Scholar

39. Øksendal B, Sulem A. Stochastic control of jump diffusions stochastic control. In: Applied Stochastic Control of Jump Diffusions. Cham: Springer (2019), p. 93–155. doi: 10.1007/978-3-030-02781-0_5

Crossref Full Text | Google Scholar

40. Capocelli R, Ricciardi L. A diffusion model for population growth in random environment. Theor Popul Biol. (1974) 5:28–41. doi: 10.1016/0040-5809(74)90050-1

PubMed Abstract | Crossref Full Text | Google Scholar

41. Turelli M. Random environments and stochastic calculus. Theor Popul Biol. (1977) 12:140–78. doi: 10.1016/0040-5809(77)90040-5

PubMed Abstract | Crossref Full Text | Google Scholar

42. Nipa KF, Jang SRJ, Allen LJ. The effect of demographic and environmental variability on disease outbreak for a dengue model with a seasonally varying vector population. Math Biosci. (2021) 331:108516. doi: 10.1016/j.mbs.2020.108516

PubMed Abstract | Crossref Full Text | Google Scholar

43. Wyse APP, Bevilacqua L, Rafikov M. Simulating malaria model for different treatment intensities in a variable environment. Ecol Modell. (2007) 206:322–30. doi: 10.1016/j.ecolmodel.2007.03.038

Crossref Full Text | Google Scholar

44. Nunno GD, Øksendal B, Proske F. Malliavin Calculus for Lévy Processes with Applications to Finance. Cham: Springer (2008). doi: 10.1007/978-3-540-78572-9

Crossref Full Text | Google Scholar

45. Parham PE, Michael E. Modelling climate change and malaria transmission. In: Michael E, Spear RC, editors. Modelling Parasite Transmission and Control. Cham: Springer (2010), p. 184–99. doi: 10.1007/978-1-4419-6064-1_13

PubMed Abstract | Crossref Full Text | Google Scholar

46. Centers for Disease Control and Prevention. Research contributions. (2024). Available at: https://www.cdc.gov/malaria/php/public-health-strategy/research-contributions.html (accessed October 04, 2024).

Google Scholar

47. Golumbeanu M, Yang GJ, Camponovo F, Stuckey EM, Hamon N, Mondy M, et al. Leveraging mathematical models of disease dynamics and machine learning to improve development of novel malaria interventions. Infect Dis Poverty. (2022) 11:37–53. doi: 10.1186/s40249-022-00981-1

PubMed Abstract | Crossref Full Text | Google Scholar

48. Ozodiegwu ID, Ambrose M, Galatas B, Runge M, Nandi A, Okuneye K, et al. Application of mathematical modelling to inform national malaria intervention planning in Nigeria. Malar J. (2023) 22:137. doi: 10.1186/s12936-023-04563-w

PubMed Abstract | Crossref Full Text | Google Scholar