Stochastic invasion thresholds for transmission-blocking symbionts and malaria control under environmental noise

1 Introduction

Malaria remains one of the most consequential vector-borne diseases globally, with the burden concentrated overwhelmingly in sub-Saharan Africa. Recent World Health Organization reports indicate that progress toward the 2025 Global Technical Strategy milestones remains substantially off track, despite major reductions in morbidity and mortality since 2000 and despite an increasingly broad intervention portfolio that includes long-lasting insecticidal nets, indoor residual spraying, improved parasitological diagnosis, artemisinin-based combination therapies, seasonal malaria chemoprevention and malaria vaccines [1, 2]. In many high-burden settings, however, the effectiveness of these tools is constrained by unequal access, heterogeneous coverage, residual and outdoor transmission, behavioral variation in mosquito–human contact, insecticide resistance, emerging antimalarial drug resistance, chronic funding gaps and fragile health systems exposed to conflict, displacement and climate-related shocks [3–8]. These limitations have made complementary biological and ecological interventions increasingly important, particularly approaches that reduce vector competence, alter mosquito population structure or suppress transmission under ecological and operational uncertainty [5, 7].

The mathematical study of malaria transmission continues to build on the Ross–Macdonald framework, which formalizes transmission in terms of mosquito density relative to humans, biting rate, vector competence, mosquito survival through the extrinsic incubation period and the duration of human infectiousness [9–12]. Modern extensions have incorporated seasonality, climate suitability, spatial heterogeneity, acquired immunity, heterogeneous exposure and intervention coverage, reflecting the fact that malaria transmission is governed by coupled epidemiological, ecological and environmental processes [13–16]. This ecological extension is essential because mosquito abundance, survival, biting behavior, larval development and parasite maturation depend nonlinearly on temperature, rainfall, humidity, and aquatic habitat availability [17–22]. From a mathematical perspective, such environmental variability can shift persistence thresholds, alter extinction probabilities and change the stationary distribution of vector populations. Consequently, stochastic population dynamics provide a natural framework for evaluating biological interventions that must establish from low frequency in noisy mosquito populations [23–26].

Symbiont-mediated vector control provides a particularly compelling setting in which stochastic threshold dynamics are biologically and operationally relevant. The most prominent example is Wolbachia-based population replacement in Aedes aegypti, where maternally inherited infections can spread through cytoplasmic incompatibility and reduce arboviral transmission under field conditions [27–29]. These successes have established microbial replacement as a viable public-health strategy, but translation to malaria vectors has been more challenging. Stable and scalable Wolbachia-based replacement has not yet been achieved in major Anopheles systems, motivating interest in native mosquito symbionts, microbiome-mediated interference, and paratransgenic approaches [30–34]. Within this broader landscape, Microsporidia MB has emerged as a promising malaria-relevant symbiont because it blocks Plasmodium falciparum transmission in Anopheles, is vertically transmitted, can spread horizontally through mating, localizes to reproductive tissues and appears to impose limited fitness costs under laboratory and semi-field conditions [35–39].

These biological features lead naturally to a two-threshold problem. First, the symbiont must invade and persist within the mosquito population through vertical inheritance, horizontal transmission or both. Second, its stationary prevalence must be sufficiently high to reduce the effective malaria reproduction number below one. These two thresholds are related but not equivalent: a symbiont may establish without reaching the prevalence required for epidemiological suppression, while a highly effective transmission blocker may fail operationally if inheritance, horizontal spread or environmental suitability is too weak. Existing mathematical work on microbial replacement, especially for Wolbachia, has clarified invasion thresholds, bistability, spatial spread, seasonal forcing and optimal release strategies [40–45]. Stochastic and control-based models have also begun to quantify how environmental noise and release design shape establishment success [46]. Nevertheless, a general stochastic threshold framework linking symbiont invasion, environmental noise, inter-class noise transfer and Ross–Macdonald malaria suppression remains underdeveloped for malaria-relevant symbionts such as Microsporidia MB.

In this paper, we develop a scale-free stochastic differential equation model for uninfected and Microsporidia MB-infected adult female Anopheles mosquitoes. The model incorporates logistic density dependence, vertical transmission, mating-mediated horizontal conversion, differential mortality and multiplicative environmental noise. From the rare-infection limit, we derive a noise-adjusted invasion index, RQstoch, which decomposes into vertical and horizontal transmission components and is penalized by infected mortality and environmental variance. We establish positivity, persistence and extinction threshold results using stochastic stability arguments, and we apply Malliavin calculus with Hörmander bracket analysis to characterize hypoelliptic regularity and noise transfer between uninfected and infected mosquito classes. Finally, we embed the stationary symbiont prevalence into a Ross–Macdonald transmission framework to define an effective malaria reproduction number and a biologically interpretable win zone in which the symbiont both establishes and makes malaria transmission subcritical. The resulting framework connects stochastic vector ecology, microbial population dynamics and malaria transmission theory, providing a mathematically rigorous basis for spatial screening, experimental prioritization and the design of future symbiont-based malaria-control strategies.

2 Materials and methods

2.1 Mathematical setting and modeling objective

We formulate a stochastic threshold model for a transmission-blocking symbiont circulating in adult female Anopheles mosquitoes. The biological motivation is Microsporidia MB, a malaria-relevant symbiont reported to reduce Plasmodium falciparum transmission and to spread through vertical and horizontal transmission routes in members of the Anopheles gambiae complex [35–38]. The mathematical objective is to connect two distinct thresholds. The first is an ecological invasion threshold: when introduced at low frequency, the symbiont must have positive long-run stochastic growth and persist in the mosquito population. The second is an epidemiological control threshold: the stationary symbiont prevalence must be sufficiently high to reduce the effective Ross–Macdonald reproduction number below one. Thus, establishment and malaria suppression are treated as related but non-equivalent events.

All stochastic processes are defined on a complete filtered probability space

(Ω,F,(Ft)t≥0,ℙ),

where the filtration (Ft)t≥0 satisfies the usual conditions, namely completeness and right-continuity. Thus, F0 contains all ℙ-null sets of F, and

Ft=⋂s>tFs, t≥0.

These conditions are standard in the theory of Itô processes and stochastic differential equations [47, 48]. The space supports two independent one-dimensional standard Brownian motions, denoted by W_*P* and W_*Q*, which represent baseline class-specific environmental fluctuations acting on the uninfected and symbiont-infected mosquito classes, respectively. The state process is a two-dimensional Itô diffusion

Xt=(Pt,Qt)⊤∈ℝ+2,

where P_*t* and Q_*t* denote uninfected and symbiont-infected adult female mosquito abundances. In compact form, the model is written as

dXt=b(Xt)dt+G(Xt)dWt, Wt=(WP,t,WQ,t)⊤,

where b:ℝ+2→ℝ2 is the drift vector and G:ℝ+2→ℝ2×2 is the diffusion matrix. For f∈C2(ℝ+2), the infinitesimal generator associated with X_*t* is (Equation 1).

Lf(x)=∇f(x)·b(x)+12Tr(G(x)G(x)⊤∇2f(x)),(1)

where ∇f is the gradient and ∇²f is the Hessian matrix.

This setting allows the invasion problem to be studied through the rare-symbiont limit near the boundary {Q = 0}, where the sign of the average stochastic growth rate of Q_*t* determines whether the symbiont can invade. It also allows the persistence problem to be formulated in terms of invariant measures on the positive state space, and the epidemiological control problem to be expressed by mapping the stationary symbiont prevalence into an effective Ross–Macdonald reproduction number. The subsequent analysis therefore proceeds from population-level stochastic dynamics to epidemiological threshold reduction.

2.2 State variables and stochastic symbiont–mosquito model

Let P(t) denote the uninfected adult female mosquito population and Q(t) denote the symbiont-infected adult female mosquito population. The total adult female abundance is

N(t)=P(t)+Q(t).

Recruitment is density-dependent with carrying capacity K_γ > 0, baseline recruitment rate r > 0, and truncated logistic factor

ϕ(u)=u+=max{u,0}.

The vertical transmission probability is τ ∈ [0, 1]. Thus, a fraction τ of offspring from infected females inherit the symbiont and enter the infected class, whereas the remaining fraction 1−τ enter the uninfected class. Horizontal transmission is represented by the effective mass-action term βPQ. This term should be interpreted as a reduced-form female-level conversion rate rather than as a literal female–female mating process. Biologically, mating-mediated dissemination involves males, but male compartments are not modeled explicitly because the epidemiological output of interest is the prevalence of the symbiont among adult females. The coefficient β therefore summarizes the average effect of infected males, mating contacts, and any other effective horizontal route after male dynamics have been averaged out. This reduction is most appropriate when male abundance and infection prevalence equilibrate on a faster time scale than female demographic change, or when male infection prevalence is sufficiently correlated with infected female abundance that the net conversion pressure can be approximated by a function of P and Q. A fully mechanistic extension would include uninfected and infected male compartments and an explicit mating function. Adult mortality rates are λ_*P*, λ_*Q* > 0. Environmental fluctuations are represented by multiplicative Brownian perturbations with intensities σ_*P*, σ_*Q*≥0.

The truncation ϕ(u) = u⁺ prevents artificial negative recruitment when total abundance exceeds the carrying-capacity scale. It is locally Lipschitz, which is sufficient for existence, uniqueness and positivity. However, it is not differentiable at u = 0. Consequently, the Malliavin–Hörmander calculations below are interpreted locally on compact subsets of the positive quadrant that do not intersect the switching hypersurface p+q = 1, or equivalently for a smooth approximation of u⁺. On the unsaturated region p+q < 1, the model reduces to the smooth algebraic logistic form ϕ(1−p−q) = 1−p−q; on the saturated region p+q > 1, ϕ(1−p−q) = 0, and the drift is again smooth. The dimensional stochastic model is

dPt=[r{Pt+(1-τ)Qt}ϕ(1-Pt+QtKγ)-λPPt-βPtQt] dt+σPPtdWP,t,dQt=[rτQtϕ(1-Pt+QtKγ)-λQQt+βPtQt]dt+σQQt dWQ,t.(2)

The conversion terms −βP_*t*Q_*t* and +βP_*t*Q_*t* preserve the total abundance N(t) and represent transfer from the uninfected class to the symbiont-infected class. The stochastic perturbations are written in Itô differential form; we do not use dW_*t*/dt, since Brownian paths are almost surely nowhere differentiable. The biological structure of system (Equation 2) is summarized in Figure 1. For compactness in the schematic, we write

Φ=ϕ(1-NKγ).

2.3 Scale-free normalization

Because field carrying capacity is difficult to estimate and varies across larval habitats, we work with carrying-capacity-normalized variables

pt=PtKγ, qt=QtKγ, nt=pt+qt.

The identifiable scaled horizontal transmission parameter is

βeff=βKγ.

The carrying capacity K_γ is used here as a local reference scale for adult female abundance rather than as an explicitly resolved larval-habitat process. After the normalization p = P/K_γ and q = Q/K_γ, the threshold analysis depends on the identifiable combination β_eff = βK_γ, and the absolute value of K_γ no longer appears explicitly in the scale-free system. In the spatial analysis, temperature, and rainfall are therefore allowed to modify recruitment and mortality, which enter directly into the rare-symbiont growth exponent, while K_γ is treated as a fixed local scaling quantity. This is a simplification: in natural systems, rainfall, land use, aquatic habitat persistence and larval-site availability may also change carrying capacity. A fully habitat-explicit extension could replace K_γ by K_γ(T, R, x), but this would require independent spatial information on larval habitat availability and is beyond the threshold-screening objective of the present study.

Biologically, β_eff is the carrying-capacity-scaled effective horizontal conversion coefficient. It combines the dimensional contact or mating-mediated conversion rate β with the local abundance scale K_γ. Because males are not modeled explicitly, β_eff should be interpreted as a compound female-level parameter summarizing the net contribution of mating-mediated and other effective horizontal routes to the conversion of uninfected adult females into symbiont-infected adult females. After normalization, system (Equation 2) becomes

dpt=[r{pt+(1-τ)qt}ϕ(1-nt)-λPpt-βeffptqt]dt +σPptdWP,t,dqt=[rτqtϕ(1-nt)-λQqt+βeffptqt]dt+σQqtdWQ,t.(3)

This normalized system is the main object of analysis. We define

A=λP+12σP2, D=λQ+12σQ2,

which are the Itô-corrected effective mortality rates of the uninfected and infected classes, respectively.

2.4 Assumptions

We use the following standing assumptions.

• (A0) Parameters satisfy
• (A1) The initial condition satisfies
• (A2) W_*P* and W_*Q* are independent standard Brownian motions adapted to (Ft).
• (A3) The resident uninfected population is viable:

The independence assumption is used as a baseline decomposition of class-specific environmental variability. It allows the direct contribution of uninfected-class and infected-class noise to be interpreted separately. In reality, shared climatic drivers may induce correlated fluctuations in both classes. A correlated-noise extension can be obtained by replacing the independence condition with d〈_{W_*P*, W_*Q*〉t} = ρ_*W*dt, where ρ_*W* ∈ [−1, 1]. We retain the independent case to keep the invasion threshold and Malliavin covariance decomposition transparent.

Assumption (Equation 4) states that the uninfected mosquito population can persist after accounting for the Itô correction induced by environmental noise.

2.5 Positivity and global existence

The first mathematical requirement is that the model remains biologically meaningful. In particular, solutions must be non-negative and must not explode in finite time.

Theorem 2.1 (Positivity and global well-posedness). *Under assumptions (A0)–(A2), for every initial condition p₀, q₀≥0, system (Equation 3) admits a unique global strong solution. Moreover,*

pt≥0, qt≥0

almost surely for all t≥0.

Proof. We write the drift of Equation 3 as b=(bp,bq)⊤, where

bp(p,q)=r{p+(1-τ)q}ϕ(1-p-q)-λPp-βeffpq,

and

bq(p,q)=rτqϕ(1-p-q)-λQq+βeffpq.

The diffusion matrix is diagonal:

G(p,q)=(σPp00σQq).

Since ϕ(u) = u⁺ = max{u, 0} is locally Lipschitz, and since the remaining terms are polynomial in (p, q), the coefficients of Equation 3 are locally Lipschitz on ℝ². Hence, for every initial condition (p0,q0)∈ℝ2, there exists a unique local strong solution up to the explosion time

ζ=limR→∞τR, τR=inf{t≥0:||(pt,qt)||≥R}.

This follows from the standard local existence and pathwise uniqueness theorem for Itô SDEs with locally Lipschitz coefficients [47, 48].

We next prove positivity. Let

x-=max{-x,0}.

For R > 0, define the localization time

τR=inf{t≥0:pt+qt≥R}.

On [0, τ_*R*], all coefficients are bounded. Applying Tanaka's formula to pt- and qt- gives

pt∧τR-=p0--∫0t∧τR1{ps≤0}dps+12Lt∧τR0(p),

and

qt∧τR-=q0--∫0t∧τR1{qs≤0}dqs+12Lt∧τR0(q),

where Lt0(p) and Lt0(q) are the local times at zero.

The local-time terms vanish. Indeed, the quadratic variations satisfy

d〈p〉t=σP2pt2dt, d〈q〉t=σQ2qt2dt.

By the occupation-density formula,

Lt0(p)=limε↓012ε∫0t1{|ps|≤ε}σP2ps2ds=0,

and similarly Lt0(q)=0, because ps2≤ε2 on {|p_*s*| ≤ ε}, and the same argument applies to q_*s*. Therefore the local time at the boundary is zero.

Substituting the SDE into Tanaka's formula and using p₀, q₀≥0, we obtain

pt∧τR-+qt∧τR-=-∫0t∧τR1{ps≤0}bp(ps,qs)ds -∫0t∧τR1{qs≤0}bq(ps,qs)ds+Mt∧τR,

where

Mt=-∫0t1{ps≤0}σPpsdWP,s-∫0t1{qs≤0}σQqsdWQ,s.

On [0, τ_*R*], M is a true martingale because its integrands are bounded.

We now bound the drift terms. Since ϕ≥0, 0 ≤ τ ≤ 1, and p_*s*, q_*s* are bounded on [0, τ_*R*], there exists a constant C_*R* > 0 such that

-1{ps≤0}bp(ps,qs)-1{qs≤0}bq(ps,qs)≤CR(ps-+qs-).

To see this explicitly, note that every term in the drift that can be positive on the negative half-space is at most linear in ps- and qs- after localization, because |p_*s*| ≤ R, |q_*s*| ≤ R, and ϕ(1−p_*s*−q_*s*) is bounded on the stopped interval. Hence the above inequality follows from local boundedness of the coefficients.

Taking expectations and using the martingale property of M_{t∧_τ*R*}, we obtain

E[pt∧τR-+qt∧τR-]≤CR∫0tE[ps∧τR-+qs∧τR-]ds.

Since p0-=q0-=0, Gronwall's inequality gives

E[pt∧τR-+qt∧τR-]=0

for every t≥0. Therefore

pt∧τR≥0, qt∧τR≥0

almost surely. Letting R → ∞, we conclude that

pt≥0, qt≥0

almost surely for all t < ζ.

It remains to rule out finite-time explosion. Let

nt=pt+qt.

On the non-negative quadrant, adding the two equations in Equation 3 cancels the horizontal conversion terms and gives

dnt=[rntϕ(1-nt)-λPpt-λQqt]dt+σPptdWP,t+σQqtdWQ,t.

Let

λmin=min{λP,λQ}.

Since p_*t*, q_*t*≥0 up to ζ,

λPpt+λQqt≥λmin(pt+qt)=λminnt.

Moreover, for n≥0,

0≤nϕ(1-n)≤1.

Indeed, if 0 ≤ n ≤ 1, then nϕ(1−n) = n(1−n) ≤ 1/4 ≤ 1, while if n > 1, then ϕ(1−n) = 0. Hence

rntϕ(1-nt)-λPpt-λQqt≤r-λminnt.

Consider the Lyapunov function

V(p,q)=1+p+q=1+n.

Its generator satisfies

LV(p,q)≤r-λminn=r+λmin-λminV(p,q).

Therefore,

LV≤c1-c2V, c1=r+λmin, c2=λmin>0.

Equivalently,

LV≤c1+c2V

with finite constants, which is sufficient for non-explosion; the stronger negative drift also gives a uniform first-moment bound.

Applying Dynkin's formula to V(p_{t∧_τ*R*}, q_{t∧_τ*R*}), we obtain

E[V(pt∧τR,qt∧τR)]=V(p0,q0)+E∫0t∧τRLV(ps,qs)ds.

Using the drift inequality

E[V(pt∧τR,qt∧τR)]≤V(p0,q0)+∫0t[c1-c2E[V(ps∧τR,qs∧τR)]]ds.

Solving this differential inequality gives

E[V(pt∧τR,qt∧τR)]≤e-c2tV(p0,q0)+c1c2(1-e-c2t).

In particular, for every fixed T > 0,

supR≥1sup0≤t≤TE[V(pt∧τR,qt∧τR)]<∞.

Suppose, for contradiction, that explosion occurs before time T with positive probability. Then, on the event {ζ ≤ T}, the stopped process must satisfy p_{τ_*R*}+q_{τ_*R*}≥R for all sufficiently large R. Hence

V(pτR,qτR)≥1+R on{τR≤T}.

By Markov's inequality,

(1+R)ℙ(τR≤T)≤E[V(pτR∧T,qτR∧T)]≤CT,

where C_*T* < ∞ is independent of R. Therefore

ℙ(τR≤T)≤CT1+R→0 asR→∞.

Since {ζ ≤ T}⊆liminf_{*R* → ∞}{τ_*R* ≤ T}, we obtain

ℙ(ζ≤T)=0.

As T > 0 was arbitrary, ζ = ∞ almost surely.

Thus the solution is global, unique, and remains in the non-negative quadrant for all time.

2.6 Resident uninfected mosquito dynamics

To derive the invasion threshold, we first study the resident system obtained by setting q_*t* = 0. The resident equation is

dpt=[(r-λP)pt-rpt2]dt+σPptdWP,t.(5)

This is a stochastic logistic diffusion on (0, ∞).

Theorem 2.2 (Resident long-run dynamics). Assume r > A, where

A=λP+12σP2.

If σ_P > 0, then the resident diffusion (Equation 5) is positive Harris recurrent on (0, ∞), admits a unique invariant probability measure π, and satisfies

Eπ[p]=1-Ar.(6)

If σ_P* = 0, then the resident dynamics are deterministic and every solution with p₀ > 0 converges to*

p⋆=1-λPr.

Proof. Suppose first that σ_*P* > 0. For p > 0, the diffusion coefficient σ_*P*p is non-zero, and the drift is locally Lipschitz. The process is therefore irreducible on (0, ∞). To establish recurrence, consider a Lyapunov function

V(p)=p+clog(1+p), c>0.

The generator associated with Equation 5 is

Lf(p)=[(r-λP)p-rp2]f′(p)+12σP2p2f″(p).

Since

V′(p)=1+c1+p, V″(p)=-c(1+p)2,

we obtain

LV(p)=[(r-λP)p-rp2](1+c1+p)-cσP2p22(1+p)2.

For large p, the dominant term is −rp², so there exist constants a, b > 0 and a compact set C⊂(0, ∞) such that

LV(p)≤-aV(p)+b1C(p).

The Foster–Lyapunov criterion implies positive Harris recurrence and uniqueness of the invariant law π [24, 49].

To compute the stationary mean, apply Itô's formula to ln p_*t*:

dln pt=[r-λP-rpt-12σP2]dt+σPdWP,t.

Starting the process in stationarity and taking expectations gives

0=r-λP-rEπ[p]-12σP2.

Solving for E_π[p] gives Equation 6. When σ_*P* = 0, Equation 5 reduces to the deterministic logistic equation

ṗt=(r-λP)pt-rpt2,

whose positive equilibrium is globally asymptotically stable whenever r > λ_*P*.

2.7 Rare-symbiont invasion exponent

When the symbiont is rare, q_*t* is close to zero and p_*t* follows the resident stationary dynamics. Applying Itô's formula to ln q_*t* in Equation 3 gives

dln qt=[rτϕ(1-pt-qt)-λQ+βeffpt-12σQ2]dt+σQdWQ,t.(7)

In the rare-symbiont limit q_*t*≈0, the long-run growth exponent is

ΛQ=rτEπ[ϕ(1-p)]-λQ+βeffEπ[p]-12σQ2.(8)

For the algebraic logistic case ϕ(1−n) = 1−n, Equation 6 gives

Eπ[1-p]=Ar.

Therefore

ΛQ=τA+βeff(1-Ar)-D.(9)

Definition 2.3 (Noise-adjusted symbiont invasion index). The stochastic invasion index of the symbiont is (Equation 8).

RQstoch=τA+βeff(1-Ar)D.(10)

The quantity RQstoch is a stochastic invasion index rather than a reproduction number in the classical next-generation-operator sense. It is obtained by rescaling the long-run per-capita logarithmic growth exponent of a rare symbiont-infected lineage in the resident uninfected environment. Thus, RQstoch>1 means that the rare infected class has positive stochastic growth when introduced at low frequency, whereas RQstoch<1 means that the rare infected class has negative stochastic growth. This ecological threshold is distinct from the Ross–Macdonald reproduction number used later to evaluate malaria transmission.

Theorem 2.4 (Stochastic invasion threshold). Under assumptions (A0)–(A3), for the algebraic logistic case ϕ(1−n) = 1−n,

ΛQ=D(RQstoch-1).

Consequently,

RQstoch>1 ⇔ ΛQ>0,

and

RQstoch<1 ⇔ ΛQ<0.

Proof. Starting from Equation 9,

ΛQ=τA+βeff(1-Ar)-D.

Factoring D > 0, we obtain

ΛQ=D[τA+βeff(1-Ar)D-1]=D(RQstoch-1).

Since D > 0, the sign of Λ_*Q* is exactly the sign of RQstoch-1.

Proposition 2.5 (Generator identity and threshold balance). *Assume the resident viability condition r > A, where*

A=λP+12σP2, D=λQ+12σQ2.

Let π denote the invariant probability measure of the resident uninfected mosquito diffusion on the boundary {q* = 0}. In the algebraic logistic case ϕ(1−n) = 1−n, the boundary generator applied to ln q satisfies*

∫L∂ln qdπ=D(RQstoch-1),

*where L∂ denotes the rare-infection generator obtained by linearizing the q-equation around q = 0. Equivalently,*

RQstoch+1λQ+12σQ2∫L∂ln(1q)dπ=1.(11)

Thus, the noise-adjusted reproductive gain and the normalized mean boundary drift of ln(1/q) exactly complement each other to one. Invasion corresponds to a negative average drift of ln(1/q), whereas extinction corresponds to a positive average drift.

Proof. Applying Itô's formula to ln q_*t* in the rare-symbiont limit gives

L∂ln q=rτ(1-p)-λQ+βeffp-12σQ2.

Averaging this expression with respect to the resident invariant distribution π yields

∫L∂ln qdπ=rτ∫(1-p)dπ-λQ+βeff∫pdπ-12σQ2.

From the resident ergodicity result,

∫pdπ=1-Ar, ∫(1-p)dπ=Ar.

Therefore,

∫L∂ln qdπ=rτAr-λQ+βeff(1-Ar)-12σQ2.

Hence

∫L∂ln qdπ=τA+βeff(1-Ar)-(λQ+12σQ2).

Since

D=λQ+12σQ2

and

RQstoch=τA+βeff(1-Ar)D,

we obtain

∫L∂ln qdπ=D(RQstoch-1).

Finally, since

ln(1q)=-ln q,

we have

∫L∂ln(1q)dπ=-D(RQstoch-1).

Dividing by D > 0 gives

1D∫L∂ln(1q)dπ=1-RQstoch,

which is equivalent to

RQstoch+1λQ+12σQ2∫L∂ln(1q)dπ=1.

This proves the identity.

Proposition 2.5 gives a useful Lyapunov interpretation of the stochastic threshold. The function ln(1/q) measures distance from symbiont establishment: it becomes large when the infected class is close to extinction and decreases when the infected class grows. Therefore, a negative mean drift of ln(1/q) along the resident boundary indicates invasion, while a positive mean drift indicates extinction. The identity (Equation 11) shows that this boundary drift is exactly balanced by the stochastic invasion index.

2.8 Biological decomposition and comparative statics

Equation 10 decomposes as

RQstoch=τAD︸vertical contribution+βeff(1-Ar)D︸horizontal contribution.(12)

The vertical component increases with inheritance fidelity τ, whereas the horizontal component increases with mating- or contact-mediated conversion β_eff and with the resident mean abundance 1−A/r. The denominator D penalizes infected mortality and infected environmental noise (Equations 12–15).

If r > A, then

∂RQstoch∂τ=AD>0, ∂RQstoch∂βeff=1-A/rD>0,

and

∂RQstoch∂σQ2=-RQstoch2D<0.

Thus infected-class noise always reduces invasion potential. The effect of resident-class noise is more subtle because it changes the resident background. Since A=λP+σP2/2,

∂RQstoch∂σP2=12D(τ-βeffr).

Therefore, resident noise increases the invasion index under vertical dominance τ > β_eff/r, but decreases it under horizontal dominance τ < β_eff/r. The balance

βeff=rτ

will also reappear in the Hörmander bracket structure below.

2.9 Extinction and persistence

The invasion exponent gives the boundary behavior of the infected class.

Theorem 2.6 (Extinction and persistence dichotomy). *Assume Equation 4. If RQstoch<1, then q_*t* → 0 exponentially in probability. If RQstoch>1, then the infected class grows when rare. Under the Foster–Lyapunov recurrence condition induced by logistic crowding and interior non-degeneracy, the process admits an invariant probability measure supported on the interior {p > 0, q > 0}.*

Proof. Integrating Equation 7 gives

qt=q0exp{∫0t[rτϕ(1-ps-qs)-λQ+βeffps-12σQ2]ds}. +σQWQ,t

When q_*t* is close to zero, ϕ(1−p_*t*−q_*t*) = ϕ(1−p_*t*)+o(1). By the ergodic theorem for the resident process and W_{*Q, t*}/t → 0 almost surely,

limt→∞1tlnqtq0=ΛQ.

If RQstoch<1, then Λ_*Q* < 0 by Theorem 2.4; hence q_*t* decays exponentially in probability. If RQstoch>1, then Λ_*Q* > 0, so the boundary q = 0 is repelling for rare infected mosquitoes. The Lyapunov bound used in Theorem 2.1 gives tightness at large total abundance. Together with irreducibility on the interior, this yields the existence of an interior invariant measure by standard Markov-process recurrence theory [24, 49].

2.10 Malliavin derivative, covariance, and local Hörmander analysis

For the scaled diffusion

dXt=b(Xt)dt+G1(Xt)dWP,t+G2(Xt)dWQ,t, Xt=(pt,qt)⊤,

with

G1(p,q)=(σPp0), G2(p,q)=(0σQq),

the Malliavin–Hörmander analysis is local. Because ϕ(u) = u⁺ is not differentiable at u = 0, the smooth vector-field calculations are made either on compact subsets of ℝ++2 that do not intersect the hypersurface p+q = 1, or for a smooth approximation of the positive-part function. In particular, on the biologically relevant unsaturated region p+q < 1, one has ϕ(1−p−q) = 1−p−q, and the drift and diffusion vector fields are smooth. We therefore work on the interior state space ℝ++2=(0,∞)2, with b,G1,G2∈C∞ on the local smooth region considered, σ_*P*, σ_*Q* > 0, and trajectories kept in ℝ++2 by the positivity result. The Lie-bracket calculations should be read as local interior calculations for the smooth logistic regime used in the threshold derivation, not as a global smoothness claim across the switching surface of the truncated model. In this setting, the Malliavin derivative is interpreted as the first-order sensitivity of the future mosquito state X_*t* to an infinitesimal environmental shock applied at time s in one of the two noise channels.

Let J_{*t, s*} = ∂X_*t*/∂X_*s* denote the Jacobian flow. It solves

Jt,s=I2+∫st∇b(Xu)Ju,sdu+∑k=12∫st∇Gk(Xu)Ju,sdWu(k),Js,s=I2.(13)

Since

∇G1=(σP000), ∇G2=(000σQ),

the Jacobian admits the factorization

Jt,s=Γt,sTexp(∫stΓu,s-1∇b(Xu)Γu,sdu),

where Texp denotes the time-ordered exponential and

Γt,s=(exp(σPΔWP;s,t-12σP2(t-s))00exp(σQΔWQ;s,t-12σQ2(t-s))),

with ΔW_{*P*; *s, t*} = W_{*P, t*}−W_{*P, s*} and ΔW_{*Q*; *s, t*} = W_{*Q, t*}−W_{*Q, s*}. This representation separates direct geometric amplification of class-specific environmental noise from subsequent ecological coupling through ∇b.

For 0 ≤ s ≤ t, the Malliavin derivatives in the two noise directions are

Ds(1)Xt=Jt,sG1(Xs), Ds(2)Xt=Jt,sG2(Xs),(14)

and Ds(i)Xt=0 for s > t. Hence

Ds(1)qt=σPpsJt,s21, Ds(2)pt=σQqsJt,s12,

measure cross-noise transfer: the first quantity is the effect on the infected class of a shock that initially perturbs the uninfected class, whereas the second measures the reverse pathway.

The Malliavin covariance matrix is

Mt=∑i=12∫0tDs(i)Xt(Ds(i)Xt)⊤ds=∫0tJt,sΣ(Xs)Jt,s⊤ds,(15)

where

Σ(Xs)=(σP2ps200σQ2qs2).

Thus Mt aggregates all past infinitesimal environmental perturbations after they have been transported through recruitment, competition, vertical transmission, and horizontal conversion. It is also useful to introduce the reduced covariance

Ct=∫0tJs,0-1Σ(Xs)(Js,0-1)⊤ds, Mt=Jt,0CtJt,0⊤,

which separates direct noise injection from subsequent stretching of the flow.

To interpret cross-class noise propagation, we next compute the first Lie bracket between the uninfected-noise vector field G₁ and the drift vector field b. Writing

b=(bp,bq)⊤, bq(p,q)=rτq(1-p-q)-λQq+βeffpq,

the Lie bracket is defined by

[G1,b]=∇bG1-∇G1b.

Since

G1(p,q)=(σPp,0)⊤,

its infected-class component is

([G1,b])q=σPp∂pbq=σPpq(βeff-rτ).

([G1,b])q=σPpq(βeff-rτ).(16)

The sign of H: = β_eff−rτ determines the leading bracket direction by which environmental variability in the uninfected class is transferred into infected-class dynamics.

Theorem 2.7 (Local Hörmander regularity on the positive quadrant). *Assume that Xt∈ℝ++2 almost surely for all t > 0. Then, for every compact set K⋐ℝ++2, the diffusion matrix*

a(x)=G1(x)G1(x)⊤+G2(x)G2(x)⊤=(σP2p200σQ2q2)

*is uniformly positive definite on K. Hence the diffusion is locally elliptic on K, the parabolic Hörmander condition is automatically satisfied there, and for every t > 0 the transition kernel admits a C^∞ density on K. Consequently, the Markov semigroup is strong Feller on ℝ++2. The bracket formula (Equation 16) is therefore not needed to prove interior smoothing when σ_*P*, σ_*Q* > 0, but it remains the correct first-order diagnostic of cross-noise transfer and becomes especially informative near the axes or in reduced-noise limits*.

Proof.[Proof sketch] Fix K⋐ℝ++2. Since p≥p_*K* > 0 and q≥q_*K* > 0 on K,

ξ⊤a(x)ξ=σP2p2ξ12+σQ2q2ξ22≥min{σP2pK2,σQ2qK2}|ξ|2, x∈K.

Thus the stopped diffusion is uniformly elliptic on K. Classical elliptic/hypoelliptic regularity therefore yields a smooth transition density; under positivity, the associated semigroup is smoothing and hence strong Feller on the interior. The bracket identity follows directly from

[G1,b]=∇bG1-∇G1b

together with

∂pbq=q(βeff-rτ).

For climate-driven screening we therefore map the bracket driver

H(T,R)=βeff-r(T,R)τ

and the effective coefficient

Ξ(T,R)=σPp⋆(T,R)q⋆(T,R)H(T,R),

because these are interpretable first-order summaries of cross-noise transfer. By contrast, the full path-dependent matrices Mt or C_*t* should be estimated only when the goal is finite-time uncertainty quantification, likelihood-based calibration, filtering, or local identifiability analysis from longitudinal abundance and prevalence data.

2.10.1 Malliavin–Hörmander interpretation

The Malliavin analysis should be interpreted as a regularity and sensitivity layer rather than as an additional epidemiological threshold. Because the model contains two direct multiplicative noise channels, local smoothing on the positive quadrant follows primarily from interior ellipticity when σ_*P*, σ_*Q* > 0 and p, q > 0. The Lie-bracket driver (Equations 17–21).

([G1,b])q=σPpq(βeff-rτ)

therefore serves mainly as a biologically interpretable first-order measure of cross-noise transfer: it quantifies how environmental perturbations to the uninfected mosquito class generate an infected-class direction through the ecological drift. For climate screening, the most interpretable quantities are consequently H(T, R) = β_eff−r(T, R)τ and Ξ(T, R) = σ_*P*p_⋆q_⋆H(T, R), whereas full estimation of the Malliavin covariance is better reserved for data-rich calibration, filtering, or local identifiability analyses [50–53].

2.11 Ross–Macdonald mapping

The stochastic invasion threshold determines whether the symbiont persists in the mosquito population. To evaluate public-health relevance, we map the stationary symbiont prevalence to an effective Ross–Macdonald reproduction number. We denote the symbiont prevalence among adult female mosquitoes by

ψt=qtpt+qt, 0≤ψt≤1.(17)

The notation ψ_*t* is used to avoid confusion with the crowding function ϕ, the invariant distribution π, and generic functions f.

Let R_{0, base} denote the baseline Ross–Macdonald reproduction number in the absence of the symbiont:

R0,base=ma2bce-μmgμhμm,(18)

where m is the mosquito-to-human ratio, a is the biting rate, b is the vector-to-human transmission probability, c is the human-to-vector transmission probability, μ_*m* is the mosquito mortality rate, g is the extrinsic incubation period, and μ_*h* is the human recovery or removal rate [10–12].

If symbiont-infected mosquitoes have relative transmission potential ρ ∈ [0, 1], then the effective malaria reproduction number in a mixed mosquito population is

R0eff(ψ)=R0,base{(1-ψ)+ρψ}.(19)

Equivalently,

R0eff(ψ)=R0,base{1-(1-ρ)ψ}.(20)

The expression above is a competence-weighted approximation. It assumes that the main epidemiological effect of the symbiont is to reduce the probability that infected mosquitoes transmit Plasmodium, while the biting rate, adult survival through the extrinsic incubation period, and mosquito-to-human ratio are represented by their baseline values. This assumption is useful for separating ecological establishment from transmission blocking, but it is not the only possible mapping. If the symbiont also changes adult survival, fecundity, larval recruitment, biting behavior, or dispersal, then those traits would enter the Ross–Macdonald expression through m, a, μ_*m*, and potentially the extrinsic-incubation survival term. In that case, the effective reproduction number should be constructed by weighting the infected and uninfected mosquito sub-populations by their trait-specific biting, survival and competence parameters rather than by competence alone.

In the competence-only case, ρ = κ, where κ ∈ [0, 1] is the relative vector competence of symbiont-infected mosquitoes. Then

R0eff(ψ)=R0,base{1-(1-κ)ψ}.(21)

Solving R0eff(ψ)<1 gives the critical symbiont prevalence required to make malaria transmission subcritical:

ψcrit=1-1R0,base1-ρ, ρ<1.(22)

Thus, malaria transmission is subcritical whenever

ψ>ψcrit.(23)

If ρ≥1, the symbiont does not reduce the mosquito transmission potential, and no prevalence level can make R0eff smaller than R_{0, base} through this mechanism.

2.12 Definition of the win zone

We define a successful symbiont-control regime, or win zone, as the simultaneous satisfaction of two thresholds:

RQstoch>1 and R0eff(ψ⋆)<1,(24)

where

ψ⋆=Eπ⋆[qp+q](25)

is the stationary symbiont prevalence under the interior invariant distribution π_⋆. The first inequality in Equation 24 is ecological: the symbiont must invade and persist. The second inequality is epidemiological: the resulting stationary symbiont prevalence must be sufficient to reduce malaria transmission below the Ross–Macdonald threshold (Equation 23).

Using Equation 22, the win zone can be written equivalently as

RQstoch>1 and ψ⋆>ψcrit.(26)

This formulation is useful because it separates the ecological problem of establishment from the epidemiological problem of transmission reduction.

Near the invasion boundary, Λ_*Q* is small and positive. Replacing p_*t* by its resident mean

p⋆=1-Ar,

the infected dynamics can be approximated by the stochastic logistic normal form

dqt≈qt(ΛQ-rτqt)dt+σQqtdWQ,t.(27)

Thus, to first order,

E[q]≈ΛQrτ.(28)

Consequently, when q is small relative to p_⋆, the stationary prevalence satisfies the approximation

ψ⋆≈ΛQrτp⋆=ΛQrτ(1-Ar).(29)

This near-threshold expression gives a transparent interpretation: stationary symbiont prevalence grows approximately linearly with the distance above the stochastic invasion threshold.

2.13 Climate-driven Monte Carlo scenario analysis

To visualize the ecological and epidemiological thresholds, we performed a climate-driven Monte Carlo scenario analysis over sub-Saharan Africa. The purpose of this analysis was not to produce a calibrated map of observed Microsporidia MB prevalence, but to illustrate how temperature, rainfall and biological uncertainty jointly affect the stochastic invasion threshold, the Malliavin–Hörmander bracket structure and the probability of symbiont-mediated malaria control (Equations 25–42).

Monthly temperature and precipitation rasters were obtained from WorldClim and clipped to sub-Saharan Africa using Natural Earth country boundaries [54, 55]. We used six representative months,

m∈{January,March,May,July,September,November},

to capture seasonal contrasts. For each raster cell x and month m, let T = T(x, m) denote mean temperature and R = R(x, m) denote monthly rainfall.

The climate-dependent adult mortality rates were represented by

λX(T)=1LX(T),LX(T)=max{LX,min,LX,max-qX(T-TX,opt)2}, X∈{P,Q}.(30)

The lower bound L_{*X*, min} prevents non-biological adult lifespans at extreme temperatures. Recruitment was modeled as a product of fecundity, temperature suitability, rainfall suitability, and development response:

r(T,R)=βEpE(T)pL(T)pP(T)SR(R)DT(T),(31)

where

pj(T)=pj,maxexp{-12(T-Tj,optwj)2}, j∈{E,L,P},

SR(R)=[max{4R(RL-R)RL2,0}]3,

and

DT(T)=(c1T+c2)exp[-(c1T+c2)]2(c1T+c2)+1.

The rainfall term was truncated at zero outside its biologically meaningful support.

For each Monte Carlo draw, uncertain biological and ecological parameters were sampled independently from the ranges shown in Table 1. The Itô-corrected mortalities were then computed as

A(T)=λP(T)+12σP2, D(T)=λQ(T)+12σQ2.(32)

The climate-dependent invasion exponent was

ΛQ(T,R)=τA(T)+βeff(1-A(T)r(T,R))-D(T),(33)

and the corresponding stochastic invasion index was

RQstoch(T,R)=τA(T)+βeff(1-A(T)r(T,R))D(T).(34)

Raster cells for which the resident viability condition r(T, R) > A(T) failed were treated as non-viable for symbiont invasion in that Monte Carlo draw.

Table 1. Uncertain parameters sampled in the spatial Monte Carlo analysis.

Ranges are illustrative and literature-informed; they are intended for stochastic threshold exploration rather than field-calibrated prediction.

To connect the spatial analysis to the Malliavin–Hörmander component, we evaluated the climate-driven bracket driver

H(T,R)=βeff-r(T,R)τ.(35)

This quantity controls the leading cross-transfer coefficient

([G1,b])q=σPpq{βeff-r(T,R)τ}=σPpqH(T,R).(36)

Thus, H(T, R) > 0 indicates dominance of horizontal conversion in the leading cross-noise transfer term, H(T, R) < 0 indicates dominance of the climate-driven recruitment and vertical-crowding contribution, and H(T, R)≈0 identifies the knife-edge where first-order transfer of environmental noise from the uninfected class to the infected class weakens.

Because the spatial maps were intended to reflect the analytical threshold rather than long transient dynamics, we used the near-threshold approximation derived from the stochastic normal form. The resident level was approximated by

p⋆(T,R)=1-A(T)r(T,R), r(T,R)>A(T),(37)

and the infected level was approximated by

q⋆(T,R)≈ΛQ(T,R)r(T,R)τ, ΛQ(T,R)>0.(38)

A small rare-introduction floor q₀ was included when RQstoch>1. The stationary symbiont prevalence proxy was then approximated by

ψ⋆(T,R)=q⋆(T,R)p⋆(T,R)+q⋆(T,R).(39)

The public-health component used the effective Ross–Macdonald expression

R0eff(ψ⋆)=R0,base{(1-ψ⋆)+ρψ⋆},(40)

where ρ ∈ [0, 1] is the relative transmission potential of symbiont-infected mosquitoes. The win event was defined as

win={RQstoch>1}∩{R0eff(ψ⋆)<1}.(41)

Equivalently, when ρ < 1, this can be written as

win={RQstoch>1}∩{ψ⋆>ψcrit}, ψcrit=1-1/R0,base1-ρ.

At each raster cell and month, we computed Monte Carlo estimates of

E[ΛQ], E[RQstoch], ℙ(RQstoch>1),

E[H], ℙ(H>ε), ℙ(|H|≤ε), E[Ξ],

E[ψ⋆], E[R0eff], ℙ(win).

The final spatial maps used N_MC = 10, 000 Monte Carlo samples per raster cell and month. A simple temperature–rainfall suitability mask was applied to remove very low-suitability cells from the main threshold maps, while the suitability score itself was also reported.

2.14 Monte Carlo stability and sensitivity analysis

To assess numerical stability, we repeated the Monte Carlo estimation on representative raster cells using

NMC∈{2,500,5,000,7,500,10,000}.

Representative cells were selected from different quantiles of the temperature–rainfall suitability distribution so that low, moderate and high environmental suitability regimes were included. For each N_MC, we computed the sample mean, variance, standard deviation, standard error and 95% Monte Carlo confidence interval for the main stochastic outputs:

ΛQ, RQstoch, 1{RQstoch>1}, H, Ξ, ψ⋆, R0eff, 1{win}.

For a generic Monte Carlo output Y₁, …, Y_*N*, we used

μ^N=1N∑i=1NYi, SE^N=Var^(Y)N,

and

μ^N±1.96SE^N.(42)

For probability outputs, Y_*i* was the corresponding binary indicator. Stability was assessed by the narrowing of confidence intervals and by the relative change in the Monte Carlo mean between successive sample sizes.

As an additional screening step, we computed rank-based sensitivity measures using Spearman correlations between sampled parameters and model outputs. For each representative raster cell, each sampled parameter θ_*j* was correlated with each output Y, giving

ρS(θj,Y).

Parameters were ranked by |ρ_*S*(θ_*j*, Y)|, and rankings were compared across representative cells to identify parameters that consistently contributed to uncertainty in the stochastic invasion exponent, symbiont prevalence, bracket coefficient, and win probability. This analysis was used only as a qualitative screening diagnostic. It does not imply causality and does not replace formal variance-based sensitivity analysis.

3 Results

3.1 Stochastic invasion threshold and climate dependence

Before interpreting the spatial outputs, we emphasize that all maps reported below are climate-driven stochastic threshold-screening maps. They are not calibrated maps of observed Microsporidia MB prevalence, mosquito abundance, or malaria burden. The biological and epidemiological parameters are sampled from plausible uncertainty ranges rather than fitted jointly to field observations. Consequently, the maps identify where model conditions are more or less favorable for invasion, prevalence build-up, and epidemiological control under the assumed parameter ranges; they should not be interpreted as predictions of current symbiont distribution or as direct forecasts of intervention impact.

The analytical threshold shows that the long-run growth rate of a rare symbiont-infected lineage is governed by

ΛQ=D(RQstoch-1),

where

RQstoch=τA+βeff(1-Ar)D.

This expression separates the invasion mechanism into a vertical-transmission component,

τAD,

and a horizontal-conversion component,

βeff(1-Ar)D.

The numerator describes the combined reproductive gain of infected mosquitoes through vertical inheritance and horizontal conversion, while the denominator D=λQ+σQ2/2 penalizes infected mortality and infected-class environmental noise.

The spatial analysis extends this threshold by allowing r, λ_*P*, and λ_*Q* to vary with temperature and rainfall. Hence, even under the same biological uncertainty ranges, the sign and magnitude of Λ_*Q*(T, R) can shift across space and season. Positive values of Λ_*Q* identify locations and months in which a rare symbiont lineage has positive stochastic growth, while negative values indicate expected loss when rare. Figure 2 shows the monthly spatial distribution of the mean invasion exponent λ_*Q*, while Figure 3 shows the corresponding stochastic invasion probability P(RQstoch>1).

3.2 Spatial probability of stochastic symbiont invasion

The probability of stochastic invasion was estimated as

ℙ(RQstoch>1)=E[1{RQstoch>1}],

where the expectation is taken over biological parameter uncertainty. This probability summarizes the robustness of invasion to uncertainty in vertical transmission, horizontal conversion, climate-driven recruitment, mortality, and environmental noise.

The resulting maps identify where the stochastic threshold is frequently satisfied under the sampled parameter ranges. Regions with high invasion probability correspond to climate–parameter combinations where recruitment, resident viability and symbiont transmission are sufficient to overcome infected mortality and noise. Regions with low invasion probability correspond to regimes where RQstoch≤1 for most Monte Carlo draws.

3.3 Spatial illustration of Malliavin–Hörmander noise transfer

The Malliavin–Hörmander analysis identified the leading cross-transfer coefficient

([G1,b])q=σPpq{βeff-r(T,R)τ}.

Therefore, the climate-dependent bracket driver

H(T,R)=βeff-r(T,R)τ

determines the sign of first-order noise transfer from the uninfected mosquito class into the infected mosquito class. When H > 0, the horizontal-conversion contribution dominates the leading cross-transfer term. When H < 0, the recruitment–vertical-transmission contribution dominates. When |H| ≤ ε, the system is close to the knife-edge where first-order cross-noise transfer weakens. Because τ and β_eff are treated as sampled biological parameters rather than climate-fitted functions, the sign of H should not be interpreted as direct evidence that climate changes maternal or mating-mediated transmission.

For the main spatial result, we mapped the effective bracket coefficient

Ξ=σPp⋆q⋆{βeff-r(T,R)τ}.

This quantity combines the bracket driver with the local stochastic state approximation. Hence, Ξ measures not only the sign of the cross-transfer direction, but also whether sufficient uninfected and infected mosquito mass is present for the transfer to be dynamically relevant. The corresponding bracket-sign maps are reported in Supplementary Figure S1. Figure 4 shows the monthly spatial distribution of the effective Lie-bracket coefficient Xi.

3.4 Win-zone probability

The stationary symbiont prevalence proxy was computed as

ψ⋆=q⋆p⋆+q⋆,

where p_⋆ and q_⋆ were obtained from the near-threshold approximation. This prevalence was then mapped to the effective Ross–Macdonald reproduction number

R0eff(ψ⋆)=R0,base{(1-ψ⋆)+ρψ⋆}.

The win event combines the ecological and epidemiological thresholds:

win={RQstoch>1}∩{R0eff(ψ⋆)<1}.

The first condition requires that the symbiont invade and persist in the mosquito population. The second requires that the resulting symbiont prevalence reduce malaria transmission below the Ross–Macdonald threshold. The win probability was therefore estimated as

ℙ(win)=E[1{RQstoch>1}1{R0eff(ψ⋆)<1}].

The maps show a clear distinction between ecological establishment and epidemiological success. Areas with high ℙ(RQstoch>1) do not necessarily have high ℙ(win), because the epidemiological condition R0eff(ψ⋆)<1 is stricter than the ecological invasion condition. This indicates that successful establishment of the symbiont is necessary but not sufficient for malaria transmission suppression; sufficiently high stationary symbiont prevalence and strong transmission blocking are also required. Figure 5 shows the monthly spatial distribution of the win probability P(win).

3.5 Monte Carlo stability

Monte Carlo stability was assessed on four representative raster cells selected across the climate-suitability gradient. The selected sites covered low-to-moderate suitability, moderate suitability, high suitability and very high suitability environments. Table 2 reports the temperature, rainfall, and suitability score for these representative cells. Before reporting the numerical stability results, Table 3 summarizes the biological interpretation of each stochastic threshold metric used in the Monte Carlo analysis.

Table 2. Representative raster cells used for Monte Carlo stability and sensitivity checks.

Suitability denotes the temperature–rainfall suitability score used for selecting representative environments.

Table 3. Biological interpretation of the stochastic threshold metrics reported in the Monte Carlo stability analysis.

For each representative cell, estimates were recomputed using

NMC∈{2,500,5,000,7,500,10,000}.

For each sample size, we recorded the mean, variance, standard error, 95% confidence interval and relative change from the previous sample size. Table 4 gives the final N_MC = 10, 000 summaries, while the full convergence table for all sample sizes is provided in Supplementary Table S1.

Table 4. Monte Carlo stability summary at N_MC = 10, 000 for representative raster cells.

Values are reported as mean with 95% Monte Carlo confidence interval in parentheses.

The representative-site analysis shows a clear separation between ecological invasion and epidemiological impact. As shown in Table 4, the invasion exponent Λ_*Q* was positive in all four representative environments, with mean values ranging from 0.3058 to 0.3946 at N_MC = 10, 000. The stochastic invasion probability was also high, ranging from 0.9461 to 1.0000. Thus, under the sampled parameter ranges, symbiont establishment was frequently supported by the ecological threshold.

However, the win probability was much lower, ranging from 0.0421 to 0.3329. This reduction occurred because the epidemiological requirement R0eff(ψ⋆)<1 is stricter than the ecological requirement RQstoch>1. In particular, Sites 2–4 had high invasion probabilities but lower win probabilities because the mean effective reproduction number remained above one. This suggests that symbiont establishment alone is not sufficient for malaria suppression unless the stationary symbiont prevalence and transmission-blocking effect are large enough.

The bracket quantities also varied across the suitability gradient. Site 1 had a positive mean bracket driver, H = 0.1649, indicating dominance of horizontal conversion in the leading cross-transfer term. Sites 2–4 had negative mean bracket drivers, indicating dominance of the climate-driven recruitment and vertical-crowding contribution. The effective bracket coefficient Ξ followed the same sign pattern. This confirms that the Malliavin–Hörmander noise-transfer structure is not spatially uniform, but changes with the local temperature–rainfall environment.

To summarize numerical precision, Table 5 compares selected confidence interval widths at N_MC = 2, 500 and N_MC = 10, 000. Increasing the sample size reduced the confidence interval widths by approximately one half for the main continuous quantities, as expected from the NMC-1/2 scaling of Monte Carlo standard errors. For example, at Site 1 the confidence interval width for Λ_*Q* decreased from 0.0146 to 0.0072, and the confidence interval width for ℙ(win) decreased from 0.0370 to 0.0185. Similar reductions were observed across Sites 2–4, supporting the numerical stability of the reported N_MC = 10, 000 maps.

Table 5. Reduction in Monte Carlo confidence interval width between N_MC = 2, 500 and N_MC = 10, 000 for selected outputs.

The reduction is consistent with the expected NMC-1/2 decrease in Monte Carlo standard error.

3.6 Numerical simulation of dual-population dynamics

To complement the analytical threshold results and the climate-driven spatial maps, we simulated the full stochastic dual-population system for uninfected and symbiont-infected mosquitoes. This numerical experiment was designed to show the transient stochastic dynamics underlying the threshold quantities and to illustrate how vertical transmission and horizontal conversion jointly shape stationary symbiont prevalence and malaria-control probability.

For the dynamic simulation, we used a representative suitable-environment parameter set and simulated N_paths = 4, 000 stochastic trajectories of the normalized system

dpt=[r{pt+(1-τ)qt}(1-pt-qt)-λPpt-βeffptqt]dt +σPptdWP,t,

dqt=[rτqt(1-pt-qt)-λQqt+βeffptqt]dt+σQqtdWQ,t.

The symbiont prevalence was computed as

ψt=qtpt+qt.

The initial conditions were p₀ = 0.65 and q₀ = 0.02, representing a rare symbiont introduction into a resident mosquito population. The baseline scenario used

τ=0.78, βeff=0.35, r=0.42,

λP=0.095, λQ=0.105, σP=0.10, σQ=0.12.

The epidemiological parameters were R_{0, base} = 1.60 and ρ = 0.25. These values were chosen as an illustrative suitable-environment scenario and were not intended to represent a calibrated field site.

Figure 6A shows the Monte Carlo mean trajectories of p_*t*, q_*t*, and ψ_*t*, with shaded 95% simulation bands. The infected class increased from rare introduction and stabilized at a positive level, while the prevalence trajectory approached an apparent stationary regime. This provides a dynamic numerical illustration of the analytical condition RQstoch>1: when the stochastic invasion threshold is satisfied, the symbiont-infected class can persist in the full dual-population stochastic system.

We next evaluated a two-dimensional parameter-space scenario over vertical transmission probability τ and effective horizontal conversion rate β_eff. For each grid point, all remaining uncertain parameters were sampled from the same ranges used in the Monte Carlo threshold analysis. To maintain consistency with the spatial maps, the heat map in Figure 6B was computed using N_MC = 10, 000 Monte Carlo samples per (τ, β_eff) grid point. The plotted quantity is the win probability

ℙ(win)=ℙ({RQstoch>1}∩{R0eff(ψ⋆)<1}).

The heat map shows that the probability of successful symbiont-mediated control increases when both vertical transmission and horizontal conversion are sufficiently high. Regions with low vertical transmission or weak horizontal conversion have lower win probability because the symbiont either fails to invade or does not reach a prevalence high enough to reduce R0eff below one. Thus, the heat map confirms that the epidemiological threshold is stricter than the ecological invasion threshold.

Finally, we generated N_MC = 20, 000 Monte Carlo samples of the stationary prevalence proxy ψ_⋆ under the baseline (τ, β_eff) scenario. Figure 6C shows the resulting distribution, together with the mean critical prevalence

ψcrit=1-1/R0,base1-ρ.

The position of the prevalence distribution relative to ψ_crit explains the difference between ecological invasion and epidemiological success. Even when RQstoch>1, malaria suppression requires that the realized stationary prevalence exceeds the critical prevalence needed for R0eff<1. Therefore, symbiont establishment alone is not sufficient: the stationary prevalence and transmission-blocking effect must also be large enough.

The numerical simulation supports the interpretation of the spatial results. The stochastic invasion probability measures whether the symbiont can establish in the mosquito population, whereas the win probability measures whether establishment is sufficient to reduce malaria transmission below the Ross–Macdonald threshold. The dual-population simulation therefore provides a dynamic bridge between the analytical threshold, the parameter-space scenario analysis and the climate-driven sub-Saharan Africa maps.

4 Discussion

The principal problem addressed in this study is that ecological establishment of a transmission-blocking symbiont does not automatically imply epidemiological control of malaria. A symbiont may invade and persist in a mosquito population but remain at a prevalence too low to reduce malaria transmission, whereas a highly effective transmission blocker may fail to establish if inheritance, horizontal spread, survival, or environmental suitability are insufficient. This distinction is especially important for Microsporidia MB, whose potential as a malaria-control agent depends jointly on vertical transmission, mating-mediated dissemination, mosquito fitness effects, environmental conditions, and the magnitude of parasite blocking. The framework developed here separates these two requirements through the ecological invasion condition

RQstoch>1

and the epidemiological control condition

R0eff(ψ⋆)<1.

The first condition is derived from the long-run stochastic growth rate of a rare symbiont lineage in the resident mosquito environment, whereas the second evaluates whether the stationary symbiont prevalence is sufficiently high to make malaria transmission subcritical. This two-threshold structure is more informative than a single deterministic establishment criterion because stochastic persistence is governed by the long-run per-capita logarithmic growth rate rather than by arithmetic mean growth alone [56–58].

The use of stochastic differential equations is justified by the ecological setting of the problem. Mosquito recruitment, mortality, breeding-site availability, and adult abundance vary through time under the influence of temperature, rainfall, habitat conditions, and demographic variability. These fluctuations alter not only trajectory variance but also the invasion threshold itself through the Itô correction. In the present model, the infected-class variance enters through

D=λQ+12σQ2,

so stronger environmental variability acting on infected mosquitoes reduces their geometric mean growth rate and therefore lowers invasion potential. Resident-class noise acts more indirectly by changing the stationary environment in which rare infected mosquitoes reproduce and convert uninfected individuals. This feature distinguishes the present framework from deterministic replacement models in which climate is often represented only through seasonal forcing or fixed scenarios. The stochastic formulation is therefore essential for quantifying establishment risk, uncertainty, and the difference between asymptotic invasion criteria and finite-time establishment success [56, 58].

The derived invasion index is also experimentally interpretable because it separates the principal biological mechanisms that promote or inhibit establishment. Vertical transmission contributes through τA, horizontal conversion through

βeff(1-Ar),

and losses through infected mortality and environmental variance. This decomposition shows that successful establishment can be improved by stronger maternal transmission, more efficient mating-mediated conversion, lower infected-mosquito mortality, or reduced environmental sensitivity. However, the relative importance of these mechanisms depends on the local recruitment environment. This conclusion is consistent with empirical evidence showing that Microsporidia MB prevalence and dissemination vary across mosquito taxa, ecological settings, and environmental gradients, and that long-term maintenance requires favorable transmission together with limited fitness penalties [35, 37–39]. The model therefore supports treating τ, β_eff, λ_*Q*, and the relative transmission parameter ρ as site-, season-, and population-specific quantities rather than as universal constants transferable across the whole region without adjustment.

The Malliavin–Hörmander analysis provides an additional theoretical contribution by identifying how stochastic forcing is transferred between mosquito classes. The leading cross-noise term is

([G1,b])q=σPpq(βeff-rτ),

which motivates the bracket driver

H=βeff-rτ.

When H > 0, the horizontal-conversion contribution dominates the vertical-recruitment coupling in the leading bracket term; when H < 0, the vertical-recruitment contribution dominates; and H≈0 identifies a first-order cancellation or balance regime. This extends the role of Hörmander theory beyond regularity of transition densities by giving the bracket a biologically meaningful interpretation as a measure of inter-class noise transfer. Nevertheless, the sign of H should not be interpreted as direct evidence that temperature or rainfall increases mating or maternal transmission. In the present formulation, climate affects H mainly through the recruitment function r(T, R), while τ and β_eff are sampled independently of climate. Direct climatic interpretation of vertical or horizontal transmission would require experimentally supported functions τ(T, R) and β_eff(T, R).

The epidemiological mapping further demonstrates why biological establishment alone is an incomplete measure of intervention success. The critical prevalence

ψcrit=1-1/R0,base1-ρ

increases with baseline malaria transmission and decreases as the symbiont's blocking effect becomes stronger. Consequently, the same symbiont may have substantial public-health value in one transmission setting but limited value in another, even when ecological invasion probabilities are similar. This result is consistent with the Ross–Macdonald literature, in which transmission thresholds depend jointly on vector abundance, survival, competence, and local transmission intensity [11]. It also complements previous microbial-replacement studies, particularly those developed for Wolbachia, which have emphasized bistability, release thresholds, traveling waves, spatial spread, and optimized dissemination strategies [40–42, 45]. The present framework differs by linking stochastic symbiont invasion to a stationary-prevalence requirement for malaria suppression and by incorporating noise-transfer diagnostics through Malliavin–Hörmander analysis.

The sub-Saharan African maps illustrate how the framework can be used as a regional threshold-screening tool. The mapped quantities are not observed Microsporidia MB prevalence surfaces and should not be interpreted as direct forecasts of intervention impact. Rather, they indicate where ecological invasion, bracket-regime dominance, stationary symbiont prevalence, and malaria-control success are more or less plausible under the assumed biological parameter ranges and climate-dependent recruitment and mortality functions. This distinction is important because contemporary malaria cartography estimates realized burden using surveillance, intervention, demographic, and environmental data, whereas the maps developed here quantify mechanistic threshold plausibility [59–62]. Their most defensible use is therefore to identify broad ecological zones where additional entomological and symbiont data would be especially informative. Areas with high invasion probability but low win probability should be prioritized for measurements of transmission-blocking efficacy and long-run symbiont prevalence, whereas areas close to H = 0 require more precise estimation of vertical transmission, horizontal conversion, and their seasonal variation.

Several biological and modeling limitations must be considered when applying the results. First, the model tracks adult female mosquitoes and represents horizontal spread through an effective conversion term. This mass-action term is biologically motivated by mating-mediated or contact-mediated dissemination, but it is also chosen because it gives a tractable first-order incidence form for threshold analysis. It does not explicitly include uninfected and infected males, mating competition, mate choice, male survival, or sex-specific symbiont transmission. Consequently, β_eff should be interpreted as a compound parameter summarizing the net female-level effect of mating-mediated or contact-mediated spread. Second, the carrying capacity used in the scale-free normalization is treated as a fixed local reference scale. This is appropriate for deriving an identifiable threshold in normalized variables, but it does not resolve rainfall-driven changes in aquatic habitat availability. Future spatial versions should couple the adult symbiont model to larval-habitat data or use a climate- and land-use-dependent carrying capacity K_γ(T, R, x). Third, the near-threshold approximation used for q_⋆ and ψ_⋆ is analytically useful but may be inaccurate far from the invasion boundary. Long-run stochastic simulation, numerical invariant-measure estimation, or likelihood-based inference should therefore replace this approximation when calibrated data become available. Fourth, the Ross–Macdonald coupling assumes that the dominant epidemiological action of the symbiont is reduced vector competence. If the symbiont also affects biting behavior, adult survival, fecundity, dispersal, or larval recruitment, then R0eff may not be a simple competence-weighted linear function of symbiont prevalence. Fifth, the climate functions represent recruitment and survival suitability rather than observed mosquito abundance or confirmed malaria transmission; climatically masked regions should not be described as malaria-free without independent epidemiological evidence. Sixth, τ, β_eff, and ρ are treated as uncertain but not dynamically climate-dependent, which may underestimate interactions between temperature, mating activity, maternal transmission, and parasite blocking. Seventh, the model omits larval-stage dynamics, mosquito dispersal, release pulses, species composition, human infection compartments, and feedback between intervention coverage and transmission. Finally, the Spearman sensitivity analysis is appropriate for preliminary screening but does not provide a full variance decomposition or identify higher-order parameter interactions. These limitations mean that the outputs should be interpreted as first-stage threshold-screening results rather than as calibrated maps of symbiont prevalence or malaria impact.

The results lead to several practical recommendations. Field and semi-field studies across contrasting ecological zones of sub-Saharan Africa should jointly estimate maternal transmission, mating-mediated conversion, infected-mosquito survival, symbiont prevalence, and transmission-blocking efficacy across seasons rather than measuring each trait independently. Site-specific observations could then replace broad parameter ranges with locally calibrated distributions and allow direct estimation of climate-conditioned functions such as τ(T, R), β_eff(T, R), and λ_*Q*(T, R). Spatial prioritization should use invasion probability and win probability together: the former identifies environments favorable for ecological establishment, whereas the latter identifies environments where establishment is also likely to reduce malaria transmission below threshold. Regions with high invasion probability but low win probability may require stronger transmission blocking, supplementary releases, or integration with existing vector-control tools. Future theoretical work should replace the near-threshold prevalence approximation with numerically estimated invariant distributions, incorporate explicit seasonal forcing, correlated environmental noise, mosquito dispersal, species heterogeneity, and release dynamics, and examine spatially optimized dissemination strategies. Future statistical work should combine the stochastic model with hierarchical Bayesian or state-space methods using climate, mosquito abundance, symbiont prevalence, and malaria surveillance data, while variance-based global sensitivity analysis should quantify nonlinear interactions among transmission, mortality, recruitment, and noise parameters. These developments would move the present threshold-screening framework toward a calibrated regional decision-support system for experimental design, ecological surveillance, and symbiont-based malaria-control planning.

5 Conclusion

This study develops a stochastic threshold framework for transmission-blocking symbionts in malaria vectors by separating the ecological problem of symbiont establishment from the epidemiological problem of malaria suppression. The analysis shows that persistence alone is not enough for control: vertical transmission, mating-mediated horizontal conversion, infected-mosquito mortality, and environmental variability determine whether the symbiont invades, while long-run symbiont prevalence determines whether the effective Ross–Macdonald reproduction number falls below one. Malliavin–Hörmander analysis further clarifies how environmental noise is redistributed across coupled mosquito classes, giving a mechanistic interpretation of stochastic coupling in symbiont ecology. The climate-driven maps and numerical experiments show that invasion potential and malaria-control potential need not coincide. The framework is therefore best interpreted as a threshold-screening tool for identifying favorable environments, prioritizing key biological traits, and guiding dissemination strategies under uncertainty. Future work should calibrate the model with field data on mosquito abundance, symbiont prevalence, species composition, and local climate, and extend it to dispersal, release pulses, male dynamics, and explicit human infection dynamics. This provides a concrete foundation for field-calibrated evaluation of Microsporidia MB and related transmission-blocking symbionts.

Generative AI statement

The author(s) declared that Generative AI was used in the creation of this manuscript. Generative AI was used to assist with language editing, LaTeX formatting, code refinement, and figure-caption preparation. All scientific content and conclusions were reviewed and approved by the author(s), who take full responsibility for the manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

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.2026.1924647/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Nomenclature

No new zoological names, nomenclatural acts, Life Science Identifiers (LSIDs), or Resource Identification Initiative identifiers (RRIDs) are introduced in this study. All biological terms are used in their standard epidemiological and entomological sense. The mathematical symbols used throughout the manuscript are defined when first introduced; the principal notation is summarized in Table 6.

Table 6. Principal notation used in the stochastic symbiont-mediated mosquito model.

Footnotes

Edited by: Haitao Song, Shanxi University, China

Reviewed by: Jason M. Graham, University of Scranton, United States

Abdoulaye Kabore, Université de Ouahigouya, Burkina Faso