# Optimizing insecticide deployment strategies to delay quantitative resistance in mosquito populations

Sylvère Kezeta-Bondja, Charles S Wondji, Ramsès Djidjou-Demasse

► To cite this version:

Sylvère Kezeta-Bondja, Charles S Wondji, Ramsès Djidjou-Demasse. Optimizing insecticide deployment strategies to delay quantitative resistance in mosquito populations. Journal of Mathematical Biology, 2026, 92 (37), 10.1007/s00285-026-02343-z. hal-05523015

HAL Id: hal-05523015

https://hal.science/hal-05523015v1

Submitted on 23 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Optimizing insecticide deployment strategies to delay quantitative resistance in mosquito populations

Sylvère Kezeta-Bondja<fnref n="1" />,<fnref n="2" /> · Charles S. Wondji<fnref n="1" />,<fnref n="3" /> ·

Ramsès Djidjou-Demasse<fnref n="4" />,<fnref n="5" /> ID

Received: 24 August 2025 / Revised: 11 December 2025 / Accepted: 3 January 2026

© The Author(s) 2026

## Abstract

The large-scale use of insecticides remains a cornerstone of malaria vector control, but its long-term effectiveness is undermined by the evolution of quantitative insecticide resistance (qIR) in mosquito populations. We develop and analyze a mathematical model to identify optimal deployment strategies for two insecticides that differ only in their relative efficacy against target mosquito populations. Resistance is represented as a continuous phenotypic trait influencing mosquito fecundity and mortality, and the model accounts for successive deployment periods. Our results show that when mutational variance is high, the optimal strategy is to deploy the most effective insecticide at full coverage, regardless of its relative efficacy or pre-deployment exposure history. By contrast, when mutational variance is low, optimal deployment requires a transient reduction in coverage during early periods, with a threshold effect driven by both relative efficacy and initial exposure rates. Crucially, we find that, under the hypothesis that the first insecticide is ineffective against mosquitoes, simultaneous use of both insecticides is rarely optimal. Instead, sequential deployment—using one insecticide until resistance reaches a critical threshold, followed by optimal use of the second—delays resistance evolution and improves long-term control. These findings provide a theoretical foundation for adaptive qIR management strategies aimed at prolonging the effectiveness of insecticides in malaria vector control.

<footnote>Ramsès Djidjou-Demasse
ramses.djidjoudemasse@ird.fr</footnote>

<footnote n="1">¹ Centre for Research in Infectious Diseases, Yaoundé, Yaoundé, Cameroon</footnote>

<footnote n="2">² Department of Mathematics, Faculty of Sciences, University of Yaoundé I, Yaoundé, Yaoundé, Cameroon</footnote>

<footnote n="3">³ Department of Vector Biology, Liverpool School of Tropical Medicine, Liverpool, UK</footnote>

<footnote n="4">⁴ MIVEGEC, Univ. Montpellier, CNRS, IRD, Montpellier, France</footnote>

<footnote n="5">⁵ École Polytechnique de Thiès, Thiès, Sénégal, France</footnote>
**Keywords** Nonlinear dynamical systems · Structured mathematical model · Quantitative insecticide resistance · Optimal control strategies · Adaptive insecticide management

**Mathematics Subject Classification** 37N25 · 35BXX · 92D15 · 92D45

# 1 Introduction

Vector-borne diseases (VBDs) remain a major global health threat, causing over 700,000 deaths annually and affecting nearly 80% of the world's population (World Health Organization 2017). Insecticide-based vector control remains one of the cornerstones of prevention, as emphasized in the Global Vector Control Response 2017–2030 (World Health Organization 2017). Among VBDs, malaria is particularly burdensome, especially in sub-Saharan Africa, with an estimated 263 million cases and 597,000 deaths in 2023 (WHO 2024). Control efforts have relied heavily on pyrethroid-based tools, notably insecticide-treated nets (ITNs) and indoor residual spraying (IRS), which accounted for roughly 78% of the malaria decline between 2000 and 2015 (Bhatt et al. 2015).

However, sustained pyrethroid use has driven widespread insecticide resistance (IR) in mosquito populations (Alout et al. 2017; World Health Organization 2018; Hemingway et al. 2016), undermining the effectiveness of current tools and compromising progress toward malaria elimination (World Health Organization 2012; Ranson and Lissenden 2016; Churcher et al. 2016; Killeen and Ranson 2018; Strode et al. 2014). Unlike agriculture, public health vector control has access to only a limited arsenal of insecticides (Ranson and Lissenden 2016), highlighting the urgent need for insecticide resistance management (IRM) strategies that can prolong insecticide efficacy (World Health Organization 2012). Overall, two main IRM approaches can be considered: rotation, in which insecticides with different modes of action are alternated over time, and mixture, in which different insecticides are applied simultaneously in the same location and at the same time.

The rapid rise of the IR has prompted numerous theoretical studies examining its evolutionary dynamics and the performance of IRM strategies. Many of these models, grounded in population genetics, assume a monogenic basis for IR, i.e., mutations in one or two genes confer resistance (Birget and Koella 2015; Levick et al. 2017; Gabrielkuniyoshi and Piodossantos 2017; Sudo et al. 2018; South and Hastings 2018; Mohammed-Awel and Gumel 2019, 2023). While such models have yielded important insights, they often overlook the transient evolutionary processes that precede fixation. Moreover, growing empirical evidence suggests that IR is frequently polygenic, with multiple genes contributing additively (Tatchou-Nebangwa et al. 2024; Weedall et al. 2020; Opondo et al. 2025). In this setting, IR can be viewed as a quantitative trait, represented by a continuous phenotypic value $x \in \Omega$, that influences mosquitoes life-history traits such as egg-laying (or fecundity) and mortality (Kezeta-Bondja et al. 2025).

Several studies compare the effectiveness of various IRM strategies, eg., (South and Hastings 2018; Helps et al. 2020; Madgwick and Kanitz 2022, 2024; Hobbs et al. 2023; Hobbs and Hastings 2025b, a). Collectively, these studies treated IR as a qualitative trait representing varying levels of reduced susceptibility and explore the evolutionary

<figure><img src="image_2.png" /></figure>
and operational dimensions of IRM, and provide complementary insights into how
resistance genetics, spatial heterogeneity, and insecticide properties influence strategy
effectiveness. In South and Hastings (2018); Helps et al. (2020), the authors account
for selection pressure, genetic mixing during reproduction, and mutation, using these
factors to compare the long-term performance of IRM. In Madgwick and Kanitz
(2022, 2024), previous work is extended by incorporating spatial heterogeneity and
mosquito dispersal between treated and untreated patches, with simulations assessing
how different spatial deployment patterns affect both vector suppression and the spread
of IR. Hobbs et al. (2023) developed a hybrid discrete-continuous model with seasonal
forcing to investigate how environmental variability modifies IRM outcomes. Finally,
(Hobbs and Hastings 2025a, b) incorporated stochasticity and more realistic genetic
architectures (including partial dominance and non-additive effects) into continuous-
trait models, highlighting the role of evolutionary uncertainty and trait variance in
determining insecticide durability.

In line with the quantitative nature of the IR, here, we consider a scenario in which an initial insecticide ($I_1$) has lost effectiveness due to resistance, prompting the introduction of a second insecticide ($I_2$) with a specified relative efficacy compared to $I_1$. We develop an age-structured integro-differential model to describe the transient evolution of resistance in both exposed and unexposed adult mosquitoes. Three IRM strategies are evaluated for their ability to suppress the adult mosquito population: (i) simple rotation (denoted as Rot), in which insecticides $I_1$ and $I_2$ are alternated in successive update periods; (ii) simple mosaic (denoted as Mos), in which each cycle's insecticide landscape is evenly divided, with 50% $I_1$ and 50% $I_2$; and (iii) optimal mosaic (denoted as Mos*), in which the proportions of $I_1$ and $I_2$ are optimized to minimize the mosquito population size over the deployment period. Building on earlier work on the quantitative nature of IR (Kezeta-Bondja et al. 2025), our study incorporates a continuous resistance trait within a structured modelling framework, thereby offering a novel tool for assessing both epidemiological and evolutionary outcomes.

The remainder of the paper is organized as follows. Section 2 describes the model formulation. Section 3 presents preliminary analytical results, including the existence of a unique maximal bounded semiflow and the conditions for steady states in the baseline model. Section 4 details parameterization and quantitative variables. Section 5 reports simulation outcomes, sensitivity analyses, and the influence of key parameters. Finally, Section 6 discusses the implications of our findings.

## 2 The model description

The model captures the dynamics of a mosquito population under insecticide pres-
sure. This population is captured by two developmental stages, with *E* representing
the number of eggs and *A* the number of adult mosquitoes. Two types of insecticides,
referred to as *I*₁ and *I*₂, are taken into account. At any given time *t*, the mosquito
population is partitioned into three distinct groups: those unexposed to insecticides
(represented by *E*₀(*t*, *x*) for eggs and *A*₀(*t*, *a*, *x*) for adults), those exposed to insecti-
cide *I*₁ (denoted *E*₁(*t*, *x*) and *A*₁(*t*, *a*, *x*)), and those exposed to insecticide *I*₂ (denoted
*E*₂(*t*, *x*) and *A*₂(*t*, *a*, *x*)). Here, *a* ∈ R⁺ indicates the age of adult mosquitoes, while
$x \in \Omega \subset \mathbb{R}$ corresponds to their resistance level. A comprehensive list of the model's variables and parameters is provided in Table 1.

## 2.1 Baseline model with a spatially uniform insecticide landscape

We define the following aggregate quantity:

$$E(t) = \sum_{j=0}^{2} \int_{\Omega} E_j(t, x)dx,$$

representing the total eggs population at time $t$. The reproductive capacity of the population is assumed to be regulated by a density-dependent function $H(E(t))$, given by:

$$H(E(t)) = (1 + E(t))^{-\kappa}, \quad (2.1)$$

where the constant $\kappa > 0$ modulates the strength of density dependence.

Eggs laid by adult female mosquitoes (AFMs) with resistance level $x$ experience mortality at rate $\mu_j(x)$ and hatch at rate $\gamma_j(x)$, where $j \in \{0, 1, 2\}$ denotes unexposed mosquitoes ($j = 0$), mosquitoes exposed to $I_1$ ($j = 1$), and mosquitoes exposed to $I_2$ ($j = 2$). Of the hatched eggs, only a fraction $\tau(x)$ emerges as adult mosquitoes.

We assume that the insecticide present in the environment, although affecting the life-history traits of the different developmental stages of the mosquito, does not induce mortality in eggs or larvae. Instead, it directly increases the mortality of adult mosquitoes–i.e., vectors that have emerged from the aquatic developmental stage. Consequently, insecticides $I_1$ and $I_2$ are present in the environment such that newly emerged adult vectors may be exposed to $I_1$ with probability $p_1$ or to $I_2$ with probability $p_2$. Under the assumption that an adult mosquito cannot be simultaneously exposed to both insecticides, ($p_1 + p_2$) represents the probability that an emerging vector is exposed to either $I_1$ or $I_2$, and thus ($1 - p_1 - p_2$) represents the probability that an emerging vector escapes exposure to both $I_1$ and $I_2$.

The boundary conditions for the adult compartments at age zero are thus:

$$\begin{cases} A_0(t, a = 0, x) = (1 - p_1 - p_2) \tau(x) [\gamma_0(x)E_0(t, x) + \gamma_1(x)E_1(t, x) + \gamma_2(x)E_2(t, x)], \\ A_1(t, a = 0, x) = p_1 \tau(x) [\gamma_0(x)E_0(t, x) + \gamma_1(x)E_1(t, x) + \gamma_2(x)E_2(t, x)], \\ A_2(t, a = 0, x) = p_2 \tau(x) [\gamma_0(x)E_0(t, x) + \gamma_1(x)E_1(t, x) + \gamma_2(x)E_2(t, x)]. \end{cases} \quad (2.2)$$

With these definitions, the model dynamics with a spatially uniform insecticide landscape is governed by the following age-structured integro-differential equations, for $j \in \{0, 1, 2\}$:

$$\left\{ \begin{aligned} \partial_t E_j(t, x) &= H(E(t)) \int_{\Omega} \int_0^{\infty} m_j(x, y)r(a, y)A_j(t, a, y)dady - (\mu_j(x) + \gamma_j(x))E_j(t, x), \\ (\partial_t + \partial_a) A_j(t, a, x) &= -d_j(a, x)A_j(t, a, x), \end{aligned} \right. \quad (2.3)$$
<table><caption>**Table 1** Model notations, state variables and parameters</caption><thead><tr><th colspan="2">Model notations</th><th rowspan="2">Biological meaning (Unit)</th></tr><tr><th>Notation</th><th></th></tr></thead><tbody><tr><td>t</td><td></td><td>Chronological time</td></tr><tr><td>a</td><td></td><td>Age of adult mosquito</td></tr><tr><td>x</td><td></td><td>Level of mosquito insecticide resistance</td></tr><tr><td>T&#773;</td><td></td><td>Number of years of insecticide I<sub>2</sub> deployment</td></tr><tr><th colspan="2">State variables</th><th rowspan="2">Description</th></tr><tr><th>Variable</th><th></th></tr><tr><td>E<sub>j</sub>(t, x)</td><td></td><td>Number of eggs laid by AFM in insecticide state j at time t and resistance level x</td></tr><tr><td>A<sub>j</sub>(t, a, x)</td><td></td><td>Number of adult mosquitoes in insecticide state j, aged a at time t with resistance level x</td></tr><tr><th colspan="2">Fixed parameters</th><th rowspan="2">Value [Ref]</th></tr><tr><th>Notation</th><th>Biological meaning (Unit)</th></tr><tr><td>γ<sub>j</sub>(x)</td><td>Proportion of eggs laid by AFM in insecticide state j that hatch (percentage/day)</td><td>0.75 [Assumed]</td></tr><tr><td>τ(x)</td><td>Proportion of hatched eggs laid by AFM that reach adulthood (percentage/day)</td><td>0.67 Bayoh and Lindsay (2003)</td></tr><tr><td>μ<sub>j</sub>(x)</td><td>Natural death rate of eggs laid by AFM in insecticide state j (percentage/day)</td><td>0.2 [Assumed]</td></tr><tr><td>A</td><td>Adult mosquitoes age limit (in days)</td><td>30 Carnevale and Robert (2009)</td></tr></tbody></table>
<table><caption>**Table 1 continued**</caption><thead><tr><th>Fixed parameters</th><th>Biological meaning (Unit)</th><th>Value [Ref]</th></tr></thead><tbody><tr><td>Notation</td><td></td><td></td></tr><tr><td>a<sub>L</sub></td><td>Average age at which AFMs start laying eggs (in days)</td><td>3 Carnevale and Robert (2009)</td></tr><tr><td>a<sub>L,S</sub></td><td>Average lifespan of Adult mosquitoes (in days)</td><td>21 Carnevale and Robert (2009)</td></tr><tr><td>r<sub>m</sub></td><td>Maximum number of eggs laid by the reference sensitive mosquito</td><td>200 CDC (2024)</td></tr><tr><td>r<sub>0</sub></td><td>Average number of eggs laid by the reference sensitive mosquito</td><td>0.875 × r<sub>m</sub> [Assumed]</td></tr><tr><td>r<sub>1</sub></td><td>Average number of eggs laid by the reference resistant mosquito</td><td>0.750 × r<sub>m</sub> [Assumed]</td></tr><tr><th>Variable parameters</th><th>Biological meaning (Unit)</th><th>Range</th></tr><tr><td>Notation</td><td></td><td></td></tr><tr><td>p<sub>ℓ</sub></td><td>Exposure rate to insecticide I<sub>ℓ</sub> (dimensionless)</td><td>(0,1)</td></tr><tr><td>r<sub>eff</sub></td><td>The relative effectiveness of insecticide I<sub>2</sub> compared to I<sub>1</sub> (dimensionless)</td><td>(0,1)</td></tr><tr><td>m<sub>j</sub>(x, y)</td><td>Probability of mutation from IR level y to x for AFM in insecticide state j (dimensionless)</td><td>defined by (4.1)</td></tr><tr><td>r<sub>(a, x)</sub></td><td>Egg-laying rate for an AFM (number/day)</td><td>defined by (4.2)</td></tr><tr><td>d<sub>j</sub>(a, x)</td><td>Natural death rate of adult mosquitoes in insecticide state j (percentage/day)</td><td>defined by (4.3)-(4.4)</td></tr><tr><td>p<sub>1</sub><sup>S</sup></td><td>Probability of surviving I<sub>1</sub>-insecticide exposure during one day</td><td>(0,1)</td></tr><tr><td>p<sub>2</sub><sup>S</sup></td><td>Probability of surviving I<sub>2</sub>-insecticide exposure during one day</td><td>defined by (4.6)</td></tr><tr><td>var<sub>0</sub></td><td>Mutational variance within unexposed AFM population</td><td>{10<sup>-2</sup>; 2 × 10<sup>-3</sup>; 10<sup>-3</sup>}</td></tr><tr><td>var</td><td>Mutational variance ratio</td><td>{2; 3}</td></tr><tr><td>var<sub>ℓ</sub></td><td>Mutational variance of I<sub>ℓ</sub> exposed AFM population</td><td>var<sub>ℓ</sub> = var × var<sub>0</sub></td></tr><tr><td>d̄<sub>0,ℓ</sub></td><td>Death rate of reference sensitive mosquito due to I<sub>ℓ</sub> exposure</td><td>− log(p<sub>ℓ</sub><sup>S</sup>)</td></tr><tr><td>d̄<sub>ℓ</sub></td><td>Death rate of reference I<sub>ℓ</sub> resistant mosquito due to I<sub>ℓ</sub> exposure</td><td>− log(1 − p<sub>ℓ</sub><sup>S</sup>)</td></tr><tr><td>κ</td><td>Scaling positive constant</td><td>(0,1)</td></tr></tbody></table>

AFM=Adult Female Mosquito, $j \in \{0, 1, 2\}$, $ℓ \in \{1, 2\}$

<figure><img src="image_3.png" /></figure>
where the death rates of adult mosquitoes are denoted by $d_j(a, x)$, for $j \in \{0, 1, 2\}$. The term

$$H(E(t)) \int_{\Omega} \int_{0}^{\infty} m_j(x, y)r(a, y)A_j(t, a, y)dady$$

describes the rate at which eggs with resistance level $x$ are produced by AFMs of resistance level $y$ and age $a$ in insecticide state $j$. The function $m_j(x, y)$ encodes the probability of mutation from resistance level (or phenotype) $y$ to resistance level (or phenotype) $x$. Thus, mutations randomly displace strains in the phenotype space at each eggs-laying generation according to the mutation kernel $m_j$. The function $r(a, y)$ denotes the age-specific oviposition rate.

The system is complemented by the following initial conditions:

$$E_j(0, x) = E_{j,0}(x), \quad A_j(0, a, x) = A_{j,0}(a, x).$$

## 2.2 The model with a variable insecticide landscape

Model (2.2)–(2.3) is based on the assumption of a spatially and temporally uniform insecticide landscape. However, to better capture field conditions, we now consider a more general framework allowing for landscape variability. We assume that during the $T$ years of insecticide deployment period, the composition of the insecticide landscape (comprising $I_1$ and $I_2$) is revised $N$ times at regular intervals of length $D_{ur}$, at times $T_0 = 0, T_1, \dots, T_{N-1}$, such that $T_n = n \times D_{ur}$, for $n = 0, \dots, (N-1)$, with $N = \bar{T}/D_{ur}$. Therefore, during each period $[T_{n-1}, T_n]$, with $n = 1, \dots, N$, the insecticide landscape composition is defined by the proportions $p_1^n$ and $p_2^n$, representing the proportion of insecticides $I_1$ and $I_2$, respectively, with $p_1^n + p_2^n \le 1$. Consistent with previous notations, for $n \in \{1, \dots, N\}$, $E_j^n(t, \cdot) = A_j(t, \cdot)$ and $A_j^n(t, \cdot, \cdot) = A_j(t, \cdot, \cdot)$, for $t \in [T_{n-1}, T_n)$. We then have, for $n \in \{1, \dots, N\}$ and $t \in [T_{n-1}, T_n]$:

$$\left\{ \begin{array}{l} A_0^n(t, a = 0, x) = (1 - p_1^n - p_2^n) \tau(x) [\gamma_0(x)E_0^n(t, x) + \gamma_1(x)E_1^n(t, x) + \gamma_2(x)E_2^n(t, x)], \\ A_1^n(t, a = 0, x) = p_1^n \tau(x) [\gamma_0(x)E_0^n(t, x) + \gamma_1(x)E_1^n(t, x) + \gamma_2(x)E_2^n(t, x)], \\ A_2^n(t, a = 0, x) = p_2^n \tau(x) [\gamma_0(x)E_0^n(t, x) + \gamma_1(x)E_1^n(t, x) + \gamma_2(x)E_2^n(t, x)], \end{array} \right. \quad (2.4)$$

and for $j \in \{0, 1, 2\}$:

$$\left\{ \begin{array}{l} \partial_t E_j^n(t, x) = H(E^n(t)) \int_{\Omega} \int_{0}^{\infty} m_j(x, y)r(a, y)A_j^n(t, a, y)dady - (\mu_j(x) + \gamma_j(x))E_j^n(t, x), \\ (\partial_t + \partial_a) A_j^n(t, a, x) = -d_j(a, x)A_j^n(t, a, x), \end{array} \right. \quad (2.5)$$

with the initial conditions $E_j^1(0, x) = E_{j,0}(x)$, $A_j^1(0, a, x) = A_{j,0}(a, x)$ and $E_j^n(T_{n-1}, x) = E_j^{n-1}(T_{n-1}, x)$, $A_j^n(T_{n-1}, a, x) = A_j^{n-1}(T_{n-1}, a, x)$ when $n \in \{2, 3, \dots, N\}$.

Systems (2.2)–(2.3) and (2.4)–(2.5) will be considered under the following assumptions:

**Assumption 2.1** For $j \in \{0, 1, 2\}$,

1. $\tau, \mu_j$ and $\gamma_j$ are positive, continuous and bounded functions on $\Omega$.
2. The functions $r$ and $d_j$ are positive, continuous and bounded on $(0, \infty) \times \Omega$.

3. The mutation kernel $m_j$, is strictly positive almost everywhere, Lipschitz continuous, integrable, bounded on $\Omega \times \Omega$, and has a unit mass, ie., $\int_{\Omega} m_j(x, y)dx = 1$ for $y \in \Omega$.

**Assumption 2.2** The mutation kernel $m_j$, $j \in \{0, 1, 2\}$,

1. is symmetric on $\Omega \times \Omega$, ie., $m_j(x, y) = m_j(y, x)$.

2. decays rather rapidly towards infinity in the sense that $\lim_{|x| \to \infty} |x|^{\ell} m_j(x, \cdot) = 0$, for all $\ell \in \mathbb{N}$.

## 3 Preliminary results under constant insecticide exposure rate

This section is devoted to the main results of Model (2.2)–(2.3). Our analysis reveals the existence of a unique maximal bounded semiflow, a significant result that underscores the underlying dynamics of System (2.2)–(2.3) under constant exposure. Furthermore, we provide the necessary and sufficient conditions for the existence of equilibria for System (2.2)–(2.3) under constant exposure. The autonomous nature of System (2.2)–(2.3) facilitates the derivation of these analytical results just as rewriting System (2.2)–(2.3) into a compact form.

By setting, $u(t, x) = (E_0(t, x), E_1(t, x), E_2(t, x))$ and $v(t, a, x) = (A_0(t, a, x), A_1(t, a, x), A_2(t, a, x))$, Model (2.2)–(2.3) rewrites as:

$$ egin{cases} 
\partial_t u(t, x) = h(u)(t) \int_{\Omega} \int_0^{\infty} m(x, y)r(a, y)v(t, a, y)dady - \mathcal{N}(x)u(t, x), \\ 
v(t, a = 0, x) = \mathcal{P}C(x)u(t, x), \\ 
(\partial_t + \partial_a)v(t, a, x) = -\mathcal{D}(a, x)v(t, a, x). 
\end{cases} 
\quad (3.1) $$

with,

$$ m(x, y) = \operatorname{diag}(m_j(x, y))_j, \quad \mathcal{N}(x) = \operatorname{diag}(\mu_j(x) + \gamma_j(x))_j, $$

$$ \mathcal{D}(a, x) = \operatorname{diag}(d_j(a, x))_j, $$

$$ \mathcal{P} = \operatorname{diag}(1 - p_1 - p_2, p_1, p_2), \quad C(x) = \tau(x) \begin{pmatrix} \gamma_0(x) & \gamma_1(x) & \gamma_2(x) \\ \gamma_0(x) & \gamma_1(x) & \gamma_2(x) \\ \gamma_0(x) & \gamma_1(x) & \gamma_2(x) \end{pmatrix}, $$

and, where

$$ h(u)(t) = H(E(t)) = (1 + \bar{h}(u)(t))^{-\kappa}, \quad \text{with} \quad \bar{h}(u)(t) = E(t). \quad (3.2) $$

### 3.1 Semiflow generated by System (2.2)–(2.3)

To establish the global well-posedness, positivity, and dissipativity of the solutions of System (2.2)–(2.3), or equivalently of (3.1), we formulate it as an abstract Cauchy
problem. Let us introduce the Banach space $X = L^1(\Omega, \mathbb{R}^3) \times L^1(\Omega, \mathbb{R}^3) \times L^1((0, \infty) \times \Omega, \mathbb{R}^3)$, endowed with the usual product norm $\|\cdot\|_X$, as well as its positive cone $X_+ = L_+^1(\Omega, \mathbb{R}^3) \times L_+^1(\Omega, \mathbb{R}^3) \times L_+^1((0, \infty) \times \Omega, \mathbb{R}^3)$. Consider the linear operator $\mathcal{A}: D(\mathcal{A}) \subset X \to X$ defined by $D(\mathcal{A}) = W^{1,1}(\Omega, \mathbb{R}^3) \times \{0_{L^1(\Omega, \mathbb{R}^3)}\} \times W^{1,1}((0, \infty) \times \Omega, \mathbb{R}^3)$ and

$$ \mathcal{A} \begin{pmatrix} u \\ 0_{L^1} \\ v \end{pmatrix} = \begin{pmatrix} -\mathcal{N}u \\ -v(0, \cdot) \\ -\partial_a v - \mathcal{D}v \end{pmatrix}. \quad (3.3) $$

Note that $\mathcal{A}$ is not densely defined in $X$ as $\overline{D(\mathcal{A})} = X_0 \subset X$. We note $X_{0+} = X_0 \cap X_+$ the positive cone of $X_0$. Let $F: X_0 \to X$ be the non-linear map defined by:

$$ F \begin{pmatrix} u \\ 0_{L^1} \\ v \end{pmatrix} = \begin{pmatrix} h(u) \int_{\Omega} \int_{0}^{\infty} m(\cdot, y) r(a, y) v(\cdot, a, y) \, \mathrm{d}a \, \mathrm{d}y \\ \mathcal{P}Cu \\ 0_{L^1((0, \infty) \times \Omega, \mathbb{R}^3)} \end{pmatrix}. \quad (3.4) $$

By setting $z(t) = (u(t, \cdot), 0_{L^1}, v(t, \cdot, \cdot))^T$ and $z_0 = (u_0, 0_{L^1}, v_0)^T$ the associated initial condition, System (2.2)–(2.3) rewrites as the following abstract Cauchy problem:

$$ \left\{ \begin{aligned} \frac{dz(t)}{dt} &= Az(t) + F(z(t)), & t > 0 \\ z(0) &= z_0. \end{aligned} \right. \quad (3.5) $$

We have the following result.

**Theorem 3.1** Let Assumption 2.1 be satisfied. There exists a unique strongly continuous semiflow $\{T(t, \cdot) : X_0 \to X_0\}_{t \ge 0}$ such that, for each $z_0 \in X_{0+}$, the map $z \in C((0, \infty), X_{0+})$ defined by $z = T(\cdot, z_0)$ is a mild solution of (3.5), ie., $\int_0^t z(s)ds \in X_0$ and $z(t) = z_0 + \mathcal{A} \int_0^t z(\sigma)d\sigma + \int_0^t F(z(\sigma))d\sigma$ for all $t \ge 0$. Moreover the semiflow $\{T(t, \cdot)\}_{t \ge 0}$ satisfies the following properties:

1. Let $T(t, z_0) = (u(t, \cdot), 0_{L^1}, v(t, \cdot, \cdot))^T$, then the following Volterra formulation holds true

$$ v(t, a, x) = \begin{cases} \Pi(a - t, a, x) v_0(a - t, x) & \text{if } t \le a, \\ \Pi(0, a, x) \mathcal{P}C(x) u(t - a, x) & \text{if } t > a, \end{cases} $$

coupled with the u-equation of (3.1), and where, for $0 \le \alpha \le \beta$,

$$ \Pi(\alpha, \beta, x) = \mathrm{diag}(\pi_j(\alpha, \beta, x))_j, \quad \text{with } \pi_j(\alpha, \beta, x) := e^{-\int_{\alpha}^{\beta} d_j(\sigma, x) d\sigma}. \quad (3.6) $$

2. For all $z_0 \in X_0$, we have

$$ E(t) \le H^{-1}\left(\frac{\bar{c}\beta}{v}\right), \quad \forall t > 0 $$
and

$$ \sum_{j=0}^{2} \int_{\Omega} \frac{1}{\tau(x)} \int_{0}^{\infty} A_{j}(t, a, x) \mathrm{d}ax \leq w_{0}e^{-t\beta} + \frac{\alpha}{\beta} H^{-1}\left(\frac{\bar{c}\beta}{\nu}\right) (1 - e^{-t\beta}), \quad \forall t > 0 $$

where $w_0 = \sum_{j=0}^{2} \int_{\Omega} \frac{1}{\tau(x)} \int_{0}^{\infty} A_{j,0}(a, x) \mathrm{d}ax$. The constants $\nu$, $\alpha$, and $\beta$ are defined by: $\nu = \|r\|_{\infty} \max_{j \in \{0,1,2\}} \left\{ \int_{\Omega} \sup_{y} m_{j}(x, y) \mathrm{d}x \right\}$, $\alpha = \max_{j \in \{0,1,2\}} \{\|\gamma_{j}\|_{\infty}\}$, $\beta = \frac{\min_{j \in \mathcal{J}} \{\inf_{(a,x) \in \mathbb{R}_{+} \times \Omega} d_{j}\}}{\max_{j \in \mathcal{J}} \{\|\tau\|_{\infty}\}}$. The constant $\bar{c}$ satisfies: $0 < \bar{c} < \min \left( \frac{\min_{j \in \{0,1,2\}} \left( \frac{\mu_j + \gamma_j}{\alpha} \right)}{\beta}, \frac{\nu}{\beta} \right)$. Here, $\underline{f} := \inf_{x \in \Omega} f$.

3. The semiflow $\{T(t, \cdot) : X_0 \to X_0\}_{t \ge 0}$ is bounded dissipative, ie., there exists a bounded set $B \subset X_0$ that attracts any bounded set $\mathcal{R} \subset X_0$, basically there exists $\zeta = \zeta(\mathcal{R}, B) \ge 0$ such that $T(t, \mathcal{R}) \subset B$ for $t \ge \zeta$.

4. The semiflow $\{T(t, \cdot)\}_t$ is asymptotically smooth in $X_+$, ie., for any nonempty, closed, bounded and positively invariant set $B \subset X_+$, there exists a compact set $\mathcal{R} \subset X_+$ such that $\lim_{t \to \infty} d_h(T(t, B), \mathcal{R}) = 0$, where $d_h$ is the Hausdorff semi-distance defined as in Hale (2014).

**Proof of Theorem 3.1** It is straightforward to verify that the operator $\mathcal{A}$ defined in (3.3) is a Hille-Yosida operator, and that the nonlinear map $F$ defined in (3.4) is positive, continuous, and locally Lipschitz on $X_0$. Standard semigroup theory then ensures the existence and uniqueness of a mild solution to (3.5) (see (Magal and Ruan 2018; Thieme 2011, 1990)). The equivalent Volterra formulation is also well-established; see, e.g., (Iannelli 1995; Webb 1985). The proof of Theorem 3.1 closely follows that of Theorem 3.1 in Kezeta-Bondja et al. (2025).

For the boundedness, let us set

$$ E(t) = \sum_{j=0}^{2} \int_{\Omega} E_{j}(t, x) \mathrm{d}x, \quad A(t) = \sum_{j=0}^{2} \int_{\Omega} \frac{1}{\tau(x)} \int_{\mathbb{R}_{+}} A_{j}(t, a, x) \mathrm{d}ax. $$

Then, by (2.2)-(2.3), it comes

$$ \begin{aligned} \dot{E}(t) = &H(E(t)) \sum_{j=0}^{2} \int_{\Omega} \int_{\Omega} \int_{\mathbb{R}_{+}} m_{j}(x, y)r(a, y)A_{j}(t, a, y) \mathrm{d}adydx \\ & - \sum_{j=0}^{2} \int_{\Omega} (\mu_{j}(x) + \gamma_{j}(x))E_{j}(t, x) \mathrm{d}x, \end{aligned} $$

$$ \dot{A}(t) = \sum_{j=0}^{2} \int_{\Omega} \gamma_{j}(x) E_{j}(t, x) \mathrm{d}x - \int_{\Omega} \int_{\mathbb{R}_{+}} \frac{d_{j}(a, x)}{\tau(x)} A_{j}(t, a, x) \mathrm{d}dx. $$
By Assumption 2.1, we find that

$$ 
\begin{aligned}
\dot{E}(t) &\leq \nu H(E(t)) A(t) - \zeta E(t), \\
\dot{A}(t) &\leq \alpha E(t) - \beta A(t),
\end{aligned}
\qquad (3.7) $$

where $\zeta = \min_{j \in \{0, 1, 2\}} (\underline{\mu}_j + \underline{\gamma}_j)$, $\nu = \|r\|_{\infty} \max_{j \in \{0, 1, 2\}} \left\{ \int_{\Omega} \sup_y m_j(x, y) dx \right\}$, $\alpha = \max_{j \in \{0, 1, 2\}} \{\|\gamma_j\|_{\infty}\}$, and $\beta = \frac{\min_{j \in \mathcal{J}} \{\inf_{(a,x) \in \mathbb{R}_+ \times \Omega} d_j\}}{\max_{j \in \mathcal{J}} \{\|\tau\|_{\infty}\}}$. Here, $\underline{f} := \inf_{x \in \Omega} f$. Let $\bar{c} > 0$, and set $W = E + \bar{c}A$. Estimates (3.7) give

$$ \dot{W}(t) \leq (\bar{c}\alpha - \zeta) E(t) + (\nu H(E(t)) - \bar{c}\beta) A(t). \qquad (3.8) $$

Since the function $H$, introduced in (2.1), is decreasing and takes values in $(0, 1)$, we have

$$ \nu H(E(t)) - \bar{c}\beta < 0 \quad \text{iff} \quad E(t) > H^{-1}\left(\frac{\bar{c}\beta}{\nu}\right), $$

where, of course, the above estimate holds on the necessary condition $\frac{\bar{c}\beta}{\nu} \in (0, 1)$. Thus, by choosing $\bar{c}$, such that

$$ 0 < \bar{c} < \min\left(\frac{\zeta}{\alpha}, \frac{\nu}{\beta}\right), \qquad (3.9) $$

estimate (3.8) leads to $\dot{W}(t) < 0$ as soon as $E(t) > H^{-1}\left(\frac{\bar{c}\beta}{\nu}\right)$. From where we find that $E$ is ultimately bounded and, for all $t$,

$$ E(t) \leq H^{-1}\left(\frac{\bar{c}\beta}{\nu}\right), $$

with $\bar{c}$ satisfying (3.9).

Finally, by (3.7), we find for all $t \geq 0$:

$$ A(t) \leq A(0)e^{-t\beta} + \frac{\alpha}{\beta}H^{-1}\left(\frac{\bar{c}\beta}{\nu}\right)(1 - e^{-t\beta}). $$

□

Still under a constant insecticide exposure rate in System (2.2)–(2.3), other technical results–such as the asymptotic behavior and uniform persistence–can also be investigated. These analyses, however, lie beyond the primary scope of this study. They can nevertheless be addressed by closely following the methodology presented in (Kezeta-Bondja et al., 2025, Section 4).
## 3.2 Steady states of System (2.2)–(2.3)

We now turn to the existence of steady states for System (2.2)–(2.3), or equivalently, for the abstract formulation (3.5). The system always admits a mosquito-free steady state, denoted by $\mathcal{E}^0 = 0$. Additional non-trivial steady states may also exist under suitable conditions.

Let us set, for all $x \in \Omega$ and $i, j \in \{0, 1, 2\}$:

$$ \Gamma_{i,j}(x) = \frac{\tau(x)\gamma_j(x)}{\mu_j(x) + \gamma_j(x)} \int_0^\infty r(a, x)\pi_i(0, a, x)da. \quad (3.10) $$

Note that $\Gamma_{i,j}(x)$–hereafter referred to as the fitness function–represents the reproductive number of AFMs with insecticide resistance level $x$ in insecticide state $j$, given that they were previously in insecticide state $i$. In the above expression, the survival probability $\pi_i(0, a, x)$, introduced in (3.6), is the likelihood that an AFM individual with IR level $x$ survives up to age $a$. When multiplied by $r(a, x)$, it gives the individual’s contribution to egg production at age $a$. Integrating this product over all ages yields the total number of eggs produced during the individual’s lifetime. Note that, under Assumption 2.1, the functions $\Gamma_{i,j}: \Omega \to \mathbb{R}_+$ are positive and continuous.

When $i = j$, we simplify the notation by writing $\Gamma_j = \Gamma_{j,j}$, which represents the fitness of an AFM individual within the $j$-insecticide state.

The existence of steady states of System (2.2)–(2.3) are closely linked to the spectral properties of linear operators $\mathcal{L}_j: L_+^1(\Omega, \mathbb{R}^2) \to L_+^1(\Omega, \mathbb{R}^2)$ and $\mathcal{L}: L_+^1(\Omega, \mathbb{R}^3) \to L_+^1(\Omega, \mathbb{R}^3)$ defined by:

$$ \mathcal{L}_j[\varphi](\cdot) = \int_{\Omega} \operatorname{diag}((1 - p_j)m_0(\cdot, y), p_j m_j(\cdot, y)) (\Gamma_{i,\ell}(y))_{i,\ell \in \{0,j\}} \varphi(y)dy, \quad j \in \{1, 2\}, $$

and

$$ \mathcal{L}[\varphi](\cdot) = \mathcal{P} \int_{\Omega} m(\cdot, y) (\Gamma_{i,j}(y))_{i,j \in \{0,1,2\}} \phi(y)dy. \quad (3.11) $$

Under Assumptions 2.1 and 2.2, note that the linear operators $\mathcal{L}_j$, in addition to being positive, are also compact and irreducible (see (Djidjou-Demasse et al., 2017, Theorem 4.1)). Therefore, by the Krein–Rutman theorem, the spectral radius $r(\mathcal{L}_j)$ of $\mathcal{L}_j$ is positive, and the associated eigenfunction is strictly positive. The following result is then obtained.

**Theorem 3.2** Let Assumptions 2.1 and 2.2 be satisfied. Then,

1. When $p_1 = 0$ and $r(\mathcal{L}_2) > 1$, System (2.2)–(2.3) admits a steady state $\mathcal{E}_1 = (E_j^*, A_j^*)_{j=0,1,2}$ that is free of mosquitoes exposed to insecticide $I_1$, such that:
$E_1 \equiv 0, A_1 \equiv 0$, and

$$E_j^* = ((r(\mathcal{L}_2))^{1/\kappa} - 1) C_{\text{cte2}}^{-1} (\mu_j + \gamma_j)^{-1} \psi_2^j, \quad j = 0, 2,$$

$$A_0^*(a, x) = (1 - p_2) ((r(\mathcal{L}_2))^{1/\kappa} - 1) C_{\text{cte2}}^{-1} \tau(x) \pi_0(0, a, x) \sum_{j=0,2} \frac{\gamma_j}{\mu_j + \gamma_j} \psi_2^j(x),$$

$$A_2^*(a, x) = p_2 ((r(\mathcal{L}_2))^{1/\kappa} - 1) C_{\text{cte2}}^{-1} \tau(x) \pi_2(0, a, x) \sum_{j=0,2} \frac{\gamma_j}{\mu_j + \gamma_j} \psi_2^j(x),$$

where $\psi_2 = (\psi_2^j)_{j=0,2} > 0$ is the eigen vector associated to $r(\mathcal{L}_2) > 1$, and $C_{\text{cte2}} = C_{\text{cte2}}(\psi_2)$ is a positive constant define later by (3.20).

2. When $p_2 = 0$ and $r(\mathcal{L}_1) > 1$, System (2.2)-(2.3) has a steady state $\mathcal{E}_2 = (E_j^*, A_j^*)_{j=0,1,2}$ that is free of mosquitoes exposed to insecticide $I_2$, such that: $E_2 \equiv 0, A_2 \equiv 0$, and

$$E_j^* = ((r(\mathcal{L}_1))^{1/\kappa} - 1) C_{\text{cte1}}^{-1} (\mu_j + \gamma_j)^{-1} \psi_1^j, \quad j = 0, 1,$$

$$A_0^*(a, x) = (1 - p_1) ((r(\mathcal{L}_1))^{1/\kappa} - 1) C_{\text{cte1}}^{-1} \tau(x) \pi_0(0, a, x) \sum_{j=0}^{1} \frac{\gamma_j}{\mu_j + \gamma_j} \psi_1^j(x),$$

$$A_1^*(a, x) = p_1 ((r(\mathcal{L}_1))^{1/\kappa} - 1) C_{\text{cte1}}^{-1} \tau(x) \pi_1(0, a, x) \sum_{j=0}^{1} \frac{\gamma_j}{\mu_j + \gamma_j} \psi_1^j(x),$$

with $\psi_1 = (\psi_1^j)_{j=0,1} > 0$ is the eigenfunction associated to $r(\mathcal{L}_1) > 1$, and $C_{\text{cte1}} = C_{\text{cte1}}(\psi_1)$ is a positive constant define later by (3.21).

3. When $p_1, p_2 \in (0, 1)$ and $r(\mathcal{L}) > 1$, System (2.2)-(2.3) possess a steady state $\mathcal{E}^* = (E_j^*, A_j^*)_{j=0,1,2}$, with mosquitoes exposed to both insecticides, such that

$$E_j^* = ((r(\mathcal{L}))^{1/\kappa} - 1) C_{\text{cte}}^{-1} (\mu_j + \gamma_j)^{-1} \psi^j, \quad j = 0, 1, 2,$$

$$A_j^*(a, x) = \mathcal{P}_{j,j} ((r(\mathcal{L}))^{1/\kappa} - 1) C_{\text{cte}}^{-1} \tau(x) \pi_j(0, a, x) \sum_{j=0}^{2} \frac{\gamma_j}{\mu_j + \gamma_j} \psi^j(x), \quad j = 0, 1, 2,$$

where $\psi = (\psi^j)_{j=0,1,2} > 0$ is the eigenfunction associated to $r(\mathcal{L}) > 1$, and $C_{\text{cte}} = C_{\text{cte}}(\psi)$ is a positive constant define later by (3.22).

The next result proceed further towards a more explicit expression of the spectral radii $r(\mathcal{L}_j)$ and $r(\mathcal{L})$ of the operators $\mathcal{L}_j$ and $\mathcal{L}$.

**Theorem 3.3** Let Assumptions 2.1 and 2.2 be satisfied. Additionally, assume that $\gamma_\ell = \gamma$ and $\mu_\ell = \mu$, for all $\ell \in \{0, 1, 2\}$. Then, the spectral radii $r(\mathcal{L}_j)$ and $r(\mathcal{L})$ of the operators $\mathcal{L}_j$, with $j = 1, 2$, and $\mathcal{L}$ can be more explicitly expressed as:

$$r(\mathcal{L}_j) = r(\mathcal{M}_j), \text{ for } j = 1, 2, \text{ and } r(\mathcal{L}) = r(\mathcal{M}),$$
with linear operators $\mathcal{M}_j : L^1(\Omega, \mathbb{R}) \to L^1(\Omega, \mathbb{R})$ and $\mathcal{M} : L^1(\Omega, \mathbb{R}) \to L^1(\Omega, \mathbb{R})$ given by

$$ 
\begin{aligned}
\mathcal{M}_j \varphi &= \int_{\Omega} [m_0(\cdot, y)(1 - p_j)\Gamma_0(y) + m_j(\cdot, y)p_j\Gamma_j(y)] \varphi(y)dy, \\
\mathcal{M} \varphi &= \int_{\Omega} [m_0(\cdot, y)(1 - p_1 - p_2)\Gamma_0(y) + m_1(\cdot, y)p_1\Gamma_1(y) + m_2(\cdot, y)p_2\Gamma_2(y)] \varphi(y)dy.
\end{aligned}
$$

**Remark 3.4** The existence of steady states of System (2.2)-(2.3) are linked to the spectral radii $r(\mathcal{L}_j)$ and $r(\mathcal{L})$ of the operators $\mathcal{L}_j$ and $\mathcal{L}$, for mutation kernels $m_j$'s satisfying Assumptions 2.1 and 2.2. In particular, when considering a family of mutation kernels $m_j^\zeta$ that depend on a small positive parameter $\zeta$ (with $\zeta \ll 1$), $r(\mathcal{L}_j)$'s and $r(\mathcal{L})$ can be directly approximated from the fitness functions $\Gamma_j$'s. Indeed, assume that

$$ m_j^\zeta(x) = \zeta^{-1} m_j(\zeta^{-1} x), $$

where $\zeta \ll 1$ captures the variance of mutations across the phenotypic space. Let us set

$$ \Theta_{[p_1, p_2]} = (1 - p_1 - p_2)\Gamma_0 + p_1\Gamma_1 + p_2\Gamma_2. \quad (3.12) $$

The function $\Theta_{[p_1, p_2]}$ denotes the fitness of a mosquito in a mixed landscape where insecticides $I_1$ and $I_2$ are present in proportions $p_1$ and $p_2$, respectively. We then introduce the sets

$$ \Theta_{[p_1, p_2]}^{\max} = \{x \in \Omega : \Theta_{[p_1, p_2]}(x) = \|\Theta_{[p_1, p_2]}\|_{\infty}\}, \quad (3.13) $$

and let $\mathcal{M}_j^\zeta$'s as well as $\mathcal{M}^\zeta$ denote the operator $\mathcal{M}_j$'s and $\mathcal{M}$ with $m_j$ replaced by $m_j^\zeta$. According to (Djidjou-Demasse et al., 2017, Theorem 2.2), for sufficiently small $\zeta$, the spectral radii satisfied:

$$ r(\mathcal{M}_1^\zeta) = r(\mathcal{L}_1) = \Theta_{[p_1, 0]}^2(\bar{x}) + \mathcal{O}(\zeta), \quad \text{for any } \bar{x} \in \Theta_{[p_1, 0]}^{\max}, $$

$$ r(\mathcal{M}_2^\zeta) = r(\mathcal{L}_2) = \Theta_{[0, p_2]}^2(\bar{x}) + \mathcal{O}(\zeta), \quad \text{for any } \bar{x} \in \Theta_{[0, p_2]}^{\max}, $$

$$ r(\mathcal{M}^\zeta) = r(\mathcal{L}) = \Theta_{[p_1, p_2]}^2(\bar{x}) + \mathcal{O}(\zeta), \quad \text{for any } \bar{x} \in \Theta_{[p_1, p_2]}^{\max}. $$

### 3.2.1 Proof of Theorem 3.2

Finding steady states of Model (2.2)-(2.3) comes down to solving the following System:

$$ 
\left\{
\begin{array}{l}
h(u^*) \displaystyle\int_0^\infty \displaystyle\int_\Omega m(x, y)r(a, y)v^*(a, y)dady = \mathcal{N}(x)u^*(x), \\
v^*(0, x) = \mathcal{PC}(x)u^*(x), \\
\partial_a v^*(a, x) = -\mathcal{D}(a, x)v^*(a, x),
\end{array}
\right.
\quad (3.14) $$
with $u^*(\cdot) = (E_j^*(\cdot))_{j=0,1,2}$ and $v^*(\cdot, \cdot) = (A_j^*(\cdot, \cdot))_{j=0,1,2}$.

Solving (3.14) for $v^*$ yields,

$$v^*(a, x) = \Pi(0, a, x)PC(x)u^*(x), \quad (3.15)$$

where $\Pi(0, a, x)$ is defined by (3.6).

The first equation of (3.14) combined with (3.15) leads to:

$$\int_{\Omega} m(x, y) \left( \int_{0}^{\infty} r(a, y) \Pi(0, a, y) da \right) PC(y) N^{-1}(y) u^{**}(y) dy = \frac{1}{h(u^*)} u^{**}(x), \quad (3.16)$$

where $u^{**} = Nu^*$, ie.,

$$u^{**} = \mathrm{diag}(\mu_j + \gamma_j)_{j=0,1,2} (E_j^*)_{j=0,1,2}.$$

Equality (3.16) rewrites as:

$$\mathcal{L}[u^{**}](x) = \frac{1}{h(u^*)} u^{**}(x), \quad (3.17)$$

and then, the existence of $u^{**}$ is strongly linked to the spectral properties of the linear operator $\mathcal{L}$ defined by (3.11).

Let us set $u^{**} = (u_0^{**}, u_1^{**}, u_2^{**})$. We then consider three cases.

**Case 1:** $p_1 = 0$ and $p_2 \in (0, 1)$. By (3.14) and (3.10), we have $u_1^{**} \equiv 0$ (equivalently $u_1^* = E_1^* \equiv 0$) and

$$\mathcal{L}_2[(u_0^{**}, u_2^{**})] = \zeta^*(u_0^{**}, u_2^{**}), \quad (3.18)$$

with

$$\zeta^* = h(u^*)^{-1}, \quad (3.19)$$

with the function $h$ introduced by (3.2), and $\mathcal{L}_2$ the linear operator defined as

$$\mathcal{L}_2 : L_+^1(\Omega, \mathbb{R}^2) \ni \varphi \mapsto \mathcal{L}_2[\varphi](\cdot) = \int_{\Omega} m(\cdot, y) \Gamma_2(y) \varphi(y) dy \in L_+^1(\Omega, \mathbb{R}^2),$$

and

$$\Gamma_2(y) = \mathrm{diag}(1 - p_2, p_2) (\Gamma_{i,j}(y))_{i,j \in \{0,2\}}.$$

As the linear operator $\mathcal{L}_2$ is positive, compact, and irreducible (see (Djidjou-Demasse et al., 2017, Theorem 4.1)), thanks to Krein-Rutman's theorem, the spectral radius $r(\mathcal{L}_2)$ of $\mathcal{L}_2$ is positive, and there exists a function $\Psi_2 = (\psi_2^0, \psi_2^2) \in L^1(\Omega, \mathbb{R}^2)$, with $\Psi_2 > 0$ a.e., such that $\mathcal{L}_2[\Psi_2] = r(\mathcal{L}_2)\Psi_2$. Moreover, by (3.18), we have (see (Djidjou-Demasse et al., 2017, Theorem 4.1))

$$r(\mathcal{L}_2) = \zeta^*, \quad \text{and} \quad \mathrm{diag}(\mu_j + \gamma_j)_{j=0,2} (E_0^*, E_2^*) = (u_0^{**}, u_2^{**}) = c\Psi_2,$$
with c a positive constant. From where, (3.15) gives

$$
\begin{align*}
A_1^* &\equiv 0, \quad A_0^*(a, x) = c(1 - p_2)\tau(x)\pi_0(0, a, x)\omega_2(x), \\
A_2^*(a, x) &= cp_2\tau(x)\pi_2(0, a, x)\omega_2(x),
\end{align*}
$$

with

$$
\omega_2 = \frac{\gamma_0}{\mu_0 + \gamma_0} \psi_2^0 + \frac{\gamma_2}{\mu_2 + \gamma_2} \psi_2^2.
$$

Therefore, by (3.19), it comes

$$
c = \left( (r(\mathcal{L}_2))^{1/\kappa} - 1 \right) C_{\text{cte2}}^{-1},
$$

with

$$
\begin{equation}
\begin{aligned}
C_{\text{cte2}} = & \sum_{j=0,2} \int_{\Omega} (\mu_j + \gamma_j)^{-1} \psi_2^j dx + (1-p_2) \int_{\Omega} \int_0^{\infty} \pi_0(0, \cdot, \cdot) \omega_2 da dx \\
& + p_2 \int_{\Omega} \int_0^{\infty} \pi_2(0, \cdot, \cdot) \omega_2 da dx.
\end{aligned}
\tag{3.20}
\end{equation}
$$

Finally, we have

$$
\begin{align*}
E_1^* &\equiv 0, \quad A_1^* \equiv 0, \\
E_j^* &= \left( (r(\mathcal{L}_2))^{1/\kappa} - 1 \right) C_{\text{cte2}}^{-1} (\mu_j + \gamma_j)^{-1} \psi_2^j, \quad j = 0, 2, \\
A_0^*(a, x) &= (1 - p_2) \left( (r(\mathcal{L}_2))^{1/\kappa} - 1 \right) C_{\text{cte2}}^{-1} \tau(x) \pi_0(0, a, x) \omega_2(x), \\
A_2^*(a, x) &= p_2 \left( (r(\mathcal{L}_2))^{1/\kappa} - 1 \right) C_{\text{cte2}}^{-1} \tau(x) \pi_2(0, a, x) \omega_2(x).
\end{align*}
$$

Case 2: $p_1 \in (0, 1)$ and $p_2 = 0$. By (3.14) and (3.10), we have $u_2^{**} \equiv 0$ (equiva-
lently $u_2^* = E_2^* \equiv 0$) and

$$
\mathcal{L}_1[(u_0^{**}, u_1^{**})] = \xi^*(u_0^{**}, u_1^{**}),
$$

with

$\xi^* = h(u^*)^{-1},$

$$
\mathcal{L}_1 : L_+^1(\Omega, \mathbb{R}^2) \ni \varphi \mapsto \mathcal{L}_1[\varphi](\cdot) = \int_{\Omega} m(\cdot, y) \Gamma_1(y) \varphi(y) dy \in L_+^1(\Omega, \mathbb{R}^2),
$$

and

$$
\Gamma_1(y) = \operatorname{diag} (1 - p_1, p_1) (\Gamma_{i,j}(y))_{i,j \in \{0,1\}}.
$$

<figure><img src="image_4.png" /></figure>
Applying the same reasoning as for case 1, we deduce that

$$E_2^* \equiv 0, \quad A_2^* \equiv 0,$$

$$E_j^* = ((r(\mathcal{L}_1))^{1/\kappa} - 1) C_{\text{cte1}}^{-1} (\mu_j + \gamma_j)^{-1} \psi_1^j, \quad j = 0, 1,$$

$$A_0^*(a, x) = (1 - p_1) ((r(\mathcal{L}_1))^{1/\kappa} - 1) C_{\text{cte1}}^{-1} \tau(x) \pi_0(0, a, x) \omega_1(x),$$

$$A_1^*(a, x) = p_1 ((r(\mathcal{L}_1))^{1/\kappa} - 1) C_{\text{cte1}}^{-1} \tau(x) \pi_1(0, a, x) \omega_1(x),$$

with

$$\omega_1 = \frac{\gamma_0}{\mu_0 + \gamma_0} \psi_1^0 + \frac{\gamma_1}{\mu_1 + \gamma_1} \psi_1^1,$$

and

$$C_{\text{cte1}} = \sum_{j=0}^{1} \int_{\Omega} (\mu_j + \gamma_j)^{-1} \psi_1^j dx + (1 - p_1) \int_{\Omega} \int_{0}^{\infty} \pi_0(0, \cdot, \cdot) \omega_1 da dx \\ + p_1 \int_{\Omega} \int_{0}^{\infty} \pi_1(0, \cdot, \cdot) \omega_1 da dx. \qquad (3.21)$$

**Case 3:** $p_1, p_2 \in (0, 1)$. Thanks to (Djidjou-Demasse et al., 2017, Theorem 4.1) and Krein-Rutman's theorem, the spectral radius $r(\mathcal{L})$ of $\mathcal{L}$ is positive, and there exists a function $\Psi = (\psi^0, \psi^1, \psi^2) \in L^1(\Omega, \mathbb{R}^3)$, with $\Psi > 0$ a.e., such that $\mathcal{L}[\Psi] = r(\mathcal{L})\Psi$. Moreover, by (3.17)

$$r(\mathcal{L}) = \frac{1}{h(u^*)}, \quad \text{and} \quad u^* = c\mathcal{N}^{-1}\Psi,$$

with $c$ a positive constant defined by,

$$c = ((r(\mathcal{L}))^{1/\kappa} - 1) C_{\text{cte}}^{-1},$$

where,

$$C_{\text{cte}} = \bar{h}(\mathcal{N}^{-1}\psi) \qquad (3.22)$$

Finally,

$$E_j^* = ((r(\mathcal{L}))^{1/\kappa} - 1) C_{\text{cte}}^{-1} (\mu_j + \gamma_j)^{-1} \psi^j, \quad j = 0, 1, 2,$$

$$A_0^*(a, x) = (1 - p_1 - p_2) ((r(\mathcal{L}))^{1/\kappa} - 1) C_{\text{cte}}^{-1} \tau(x) \pi_0(0, a, x) \omega(x),$$

$$A_1^*(a, x) = p_1 ((r(\mathcal{L}))^{1/\kappa} - 1) C_{\text{cte}}^{-1} \tau(x) \pi_1(0, a, x) \omega(x),$$

$$A_2^*(a, x) = p_2 ((r(\mathcal{L}))^{1/\kappa} - 1) C_{\text{cte}}^{-1} \tau(x) \pi_2(0, a, x) \omega(x),$$

with

$$\omega = \frac{\gamma_0}{\mu_0 + \gamma_0} \psi^0 + \frac{\gamma_1}{\mu_1 + \gamma_1} \psi^1 + \frac{\gamma_2}{\mu_2 + \gamma_2} \psi^2.$$

This concludes the proof of Theorem 3.2.
### 3.2.2 Proof of Theorem 3.3

Under assumption $\gamma_\ell = \gamma$ and $\mu_\ell = \mu$, observe that the operator $\mathcal{L}_2$ is then of the form: $\mathcal{L}_2 = \begin{pmatrix} L_2 & L_2 \\ \bar{L}_2 & \bar{L}_2 \end{pmatrix}$, with

$$L_2\varphi = (1 - p_2) \int_{\Omega} m_0(\cdot, y) \Gamma_{0,0}(y) \varphi(y) dy, \quad \bar{L}_2\varphi = p_2 \int_{\Omega} m_2(\cdot, y) \Gamma_{2,2}(y) \varphi(y) dy.$$

To proceed further, we claim that (the proof will be provided later):

**Claim 3.5** Let $X$ be a Banach space, $M: X \to X$ and $N: X \to X$ two linear, positive, compact and irreducible operators. Set $W = \begin{pmatrix} N & N \\ M & M \end{pmatrix}$. Denote by $r(W)$ the spectral radius of $W$ as well as $\psi$ the corresponding eigenfunction. Then,

$$r(W) = r(M + N), \quad \text{and} \quad \psi = \left\| \begin{pmatrix} N\phi \\ M\phi \end{pmatrix} \right\|^{-1} \begin{pmatrix} N\phi \\ M\phi \end{pmatrix},$$

where $\phi > 0$ is the eigenfunction associated to spectral radius $r(M + N)$ of $(M + N)$.

Consequently, by Claim 3.5, we have $r(\mathcal{L}_2) = r(L_2 + \bar{L}_2)$. Similar arguments apply for the cases of $r(\mathcal{L}_1)$ and $r(\mathcal{L})$.

We conclude the proof of Theorem 3.3 by providing the proof of Claim 3.5.

**Proof of Claim 3.5** Since $W$ is a linear positive, compact and irreducible operator, by the Krein-Rutman's theorem, we can find $\psi = (u, v) > 0$, with $\|\psi\| = 1$, such that $W\psi = r(W)\psi$, with $r(W) > 0$ the spectral radius of $W$, ie.,

$$\begin{pmatrix} N(u + v) \\ M(u + v) \end{pmatrix} = r(W) \begin{pmatrix} u \\ v \end{pmatrix}. \qquad (3.23)$$

Therefore, $(N + M)(u + v) = r(W)(u + v)$.

Since $N + M$ is a positive, compact and irreducible operator, by the Krein-Rutman's theorem, we have $r(W) = r(M + N)$ and $u + v = c_0\phi$, where $c_0$ is a positive constant and $\phi$ is the eigenvector associated the spectral radius $r(M + N)$ of $N + M$.

By (3.23), we have

$$u = \frac{c_0}{r(M + N)} N\phi, \quad v = \frac{c_0}{r(M + N)} M\phi.$$

As $\|\psi\| = \|(u, v)\| = 1$, it comes

$$u = \left\| \begin{pmatrix} N\phi \\ M\phi \end{pmatrix} \right\|^{-1} N\phi, \quad v = \left\| \begin{pmatrix} N\phi \\ M\phi \end{pmatrix} \right\|^{-1} M\phi.$$

This ends the proof of Claim 3.5. □
## 4 Model parameterization and quantitative variable

The values and ranges of all constant parameters are provided in Table 1. Other functional parameters, such as the proportion of eggs laid $\gamma_j$'s and the rate at which eggs become unviable $\mu_j$'s, are also assumed to be constant, with their values specified in Table 1.

**The mutation kernel.** The mutation probability $m_j$ is defined as

$$m_j(x, y) = \mathcal{G}(0, \text{var}_j)(x - y), \quad (4.1)$$

where $\mathcal{G}(0, \text{var}_j)$ stands for the normal probability density distribution with mean 0 and standard deviation $\text{var}_j$. The mutational variance $\text{var}_0$ within the unexposed mosquito population is considered variable, with values provided in Table 1. For mosquito populations in insecticide state $j \in \{1, 2\}$, the mutational variance $\text{var}_j$ is defined as $\text{var}_j = \text{var} \times \text{var}_0$, where the mutational variance ratio $\text{var}$ is also assumed to be variable, with its values specified in Table 1. Moreover, the Gaussian kernel is used here only as an illustrative example; other kernels may be considered as long as Assumption 2.2 is satisfied.

**Egg-laying rate.** Using a similar framework as in Djidjou-Demasse et al. (2023), the egg-laying rate $r(a, x)$ of an AFM (adult female mosquito) of age $a$ and resistance level $x$ is defined by:

$$r(a, x) = r_m \left[ 1 + \left( \frac{r_m - r_0}{r_0} \right) \left( \frac{r_0}{r_1} \cdot \frac{r_m - r_1}{r_m - r_0} \right)^x \right]^{-1} \times \mathbf{1}_{a>a_L}, \quad (4.2)$$

where $r_m$ is a constant due to physiological constraints and $a_L$ is the average age at which AFMs start laying eggs. The constants $r_0$ and $r_1$ are egg-laying rates of the reference 'sensitive' and 'resistant' mosquitoes to insecticide with $0 < r_1 < r_0 < r_m$. The constants $r_i$'s are given in Table 1.

**Adult mosquitoes death rate.** The death rate $d_0(a, x)$ of unexposed adult mosquitoes aged $a$ with resistant level $x$ is such that

$$d_0(a, x) = \kappa_0 \mathbf{1}_{a>a_{LS}}, \quad (4.3)$$

where $a_{LS}$ is the average lifespan of mosquitoes, and $\kappa_0$ is a positive constant. Note that the average lifespan of unexposed mosquitoes is given by $\int_0^\infty \pi_0(0, a, x)da = a_{LS} + \frac{1}{\kappa_0}$. Since $a_{LS}$ is fixed around 21 days, therefore, as soon as the constant $\kappa_0$ is sufficiently large, e.g., $\kappa_0 = 10$, we have $\int_0^\infty \pi_0(0, a, x)da = a_{LS} + \frac{1}{\kappa_0} \approx a_{LS}$.

For adult mosquitoes in the $j$-insecticide state ($j = 1, 2$), their death rate $d_j(a, x)$ accounts for the insecticide pressure such that

$$d_j(a, x) = d_0(a, x) + \bar{d}_{0,j} \left( \frac{\bar{d}_j}{\bar{d}_{0,j}} \right)^x, \quad (4.4)$$
<figure><img src="image_5.png" /><figcaption>**Fig. 1** (A) Eggs-laying rate, *r*(*x*), of AFMs with resistance level *x*. (B) Insecticide activity, $\bar{d}_{0,j} \left( \frac{\bar{d}_j}{d_{0,j}} \right)^x$, on AFMs population with resistance level *x*</figcaption></figure>

where $d_0(a, x)$ represents the baseline mortality rate in the absence of insecticide exposure.

The second term models the additional mortality induced by insecticide $I_j$. Here, $\bar{d}_{0,j}$ and $\bar{d}_j$ correspond respectively to the insecticide-induced death rates of a reference sensitive mosquito (labeled as resistance level $x=0$) and a reference resistant mosquito (labeled resistance level $x=1$). The expression $\bar{d}_{0,j} \left( \frac{\bar{d}_j}{\bar{d}_{0,j}} \right)^x$ provides an exponential interpolation between these two extremes, allowing the insecticide-induced mortality to decrease smoothly and nonlinearly as resistance increases.

Finally, note that, in general, the function *r* is assumed to be bounded (Assumption 2.1), which ensures control of egg-laying rate as a function of the resistance phenotype *x*. Moreover, the qualitative shapes of the functions *r*(*x*, *x*) and $\bar{d}_{0,j} \left( \frac{\bar{d}_j}{\bar{d}_{0,j}} \right)^x$ are shown in Figure 1. Importantly, the specific functional forms of *r* and $d_j$ are not required for the formulation or analysis of the model; the essential assumption is that both functions monotonically decrease with respect to the resistance level *x*, capturing the fitness cost of resistance through the function *r* and the fitness advantage of resistance through the function $d_j$.

**Initial conditions before deployment of insecticide I₂.** We assume that, after a prolonged period of exclusive use of insecticide $I_1$, with an exposure rate denoted by $p_0$, the mosquito population reaches a stationary state. The parameter $p_0$ represents the overall initial insecticide pressure in the environment.

At time $t=0$, insecticide $I_2$ is introduced into the insecticide landscape. At the stationary state prior to the introduction of $I_2$, the fitness function $\Theta_{[p_0,0]}$, defined in (3.12), is maximized at the phenotypic value $\bar{x}_1 \in \Theta_{[p_0,0]}^{\max}$ within a landscape composed solely of insecticide $I_1$ at proportion $p_1 = p_0$ (Figure 2A), with $\Theta_{[p_0,0]}^{\max}$ defined in (3.13). At evolutionary equilibrium, the mosquito population concentrates around this dominant phenotypic value (Figures 2B,C).
<figure><img src="image_6.png" /><figcaption>Fig. 2 Evolutionary dynamics under exclusive use of insecticide I₁ at a proportion p₀ = 0.5. (A) The fitness functions. (B) The total mosquito population at evolutionary equilibrium is largely concentrated around the dominant phenotype x̄₁, which maximizes the fitness function Θₜₚ₀. (C) In the insecticide landscape composed only of I₁ in proportion p₁ = p₀, the dynamics of total mosquito population size and distribution across phenotypic values.</figcaption></figure>

The new insecticide I₂ is characterized such that the fitness Γ₂(x₁) of mosquitoes with resistance level x₁ under its selective pressure satisfies

$$ \Gamma_2(x_1) = (1 - r_{\text{eff}}) \Gamma_1(x_1), \quad (4.5) $$

where x₁ is the optimal phenotype with respect to the fitness function Γ₁. The parameter rₐ₇₅ ∈ (0, 1) quantifies the relative effectiveness of insecticide I₂ compared to I₁ against the initial dominant phenotype established at evolutionary equilibrium under exclusive use of I₁ (Figure 3).

Let p₁S ∈ (0, 1) be given and fixed, denoting the probability of surviving I₁-insecticide exposure during one day. We set d̃₀₁ = −log(p₁S) and d̃₁ = −log(1 − p₁S). From here, the fitness function Γ₁ defined by (3.10) is fully parameterized. For the parameterization of the fitness function Γ₂, the probability p₂S of surviving I₂-insecticide exposure during one day is then determined such that (4.5) holds, i.e., for a given value of rₐ₇₅ ∈ (0, 1), p₂S satisfies

$$ \int_{0}^{\infty} r(a, x_1) e^{-\int_{0}^{a} d_0(\sigma, x_1) d\sigma} \exp \left( a \cdot \log(p_2^S) \left( \frac{\log(1 - p_2^S)}{\log(p_2^S)} \right)^{x_1} \right) da = (1 - r_{\text{eff}}) \vartheta_1, \quad (4.6) $$

with ϑ₁ = Γ₁(x₁) (τ(x₁)γ₂(x₁) / (μ₂(x₁)+γ₂(x₁)))⁻¹. While (4.6) can be solved numerically, obtaining an explicit expression for p₂S from (4.6) is generally not straightforward.

**Deployment strategies.** Insecticide I₂ may be deployed in combination with insecticide I₁ during the T years insecticide implementation period. The composition of the insecticide landscape is updated N times at regular intervals, occurring at times T₀, T₁, ..., Tₙ₋₁ introduced previously. During each interval [Tₙ₋₁, Tₙ), with n ∈ {1, ..., N}, the landscape is characterized by proportions p₁ⁿ and p₂ⁿ, corresponding to the proportions of insecticides I₁ and I₂, respectively, with p₁ⁿ + p₂ⁿ ≤ 1.
<figure><img src="image_7.png" /><figcaption>**Fig. 3 The fitness landscape.** Prior to the introduction of insecticide I₂, the dominant mosquito phenotype is denoted by x₁. The relative effectiveness parameter r<sub>eff</sub> measures the fitness reduction of phenotype x₁ under I₂ compared to I₁, and is defined such that Γ₂(x₁) = (1 - r<sub>eff</sub>)Γ₁(x₁). When insecticide I₂ is used alone, the dominant phenotype shifts to x₂. In a mixed insecticide landscape with proportions p₁ and p₂ of I₁ and I₂, respectively, the fitness function Θ<sub>[p₁, p₂]</sub> reaches its maximum at a new dominant phenotype, denoted by x̄₂. In the scenario illustrated here, p₀ = 0.5, p₁ = 0.2, p₂ = 0.4, and r<sub>eff</sub> = 0.75.</figcaption></figure>

A deployment strategy, denoted as D<sub>epS</sub>, over a period of T̄ years is defined by the temporal composition of the insecticide landscape, given by:

$$D_{epS} = \begin{pmatrix} p_1^1 & p_1^2 & \cdots & p_1^N \\ p_2^1 & p_2^2 & \cdots & p_2^N \end{pmatrix}$$

In the sequel, we set $p_j = (p_j^1, p_j^2, \cdots, p_j^N)$ to denote the corresponding proportions of insecticides $I_j$, with $j = 1, 2$.

We consider three deployment strategies: (i) a simple rotation strategy (denoted Rot), which alternates the use of insecticides $I_1$ and $I_2$ in successive update periods; (ii) a simple mosaic strategy (denoted Mos), in which each update period consists of an equal split of the insecticide landscape – 50% $I_1$ and 50% $I_2$; and (iii) an optimal mosaic strategy (denoted Mos*), which identifies the best combination of $I_1$ and $I_2$ proportions across update periods to minimize the mosquito population size over the entire deployment period.

**Relative gain.** The relative gain, $R_{gain} = R_{gain}(D_{epS}, p_0)$, from using insecticide I₂ over a period of T̄ years is calculated as the relative reduction in mosquito population size compared to the baseline population size observed under continuous usage of
insecticide I₁ alone with exposure rate p₀. Formally, this is defined as:

$$R_{\text{gain}}(D_{\text{epS}}, p_0) = 1 - \frac{\sum_{j=0}^{2} \int_{0}^{\bar{T}} \int_{\Omega} \int_{0}^{\infty} A_j(t, a, x) \mathrm{d}a \mathrm{d}x \mathrm{d}t \Big|_{D_{\text{epS}}}}{\bar{T} \int_{\Omega} \int_{0}^{\infty} (\bar{A}_0(a, x) + \bar{A}_1(a, x)) \mathrm{d}a \mathrm{d}x \Big|_{p_0}},$$

where the numerator of the fraction denotes the size of the mosquito population during the $\bar{T}$ years deployment of insecticide I₂ with landscape composition strategy $D_{epS}$, and the denominator denotes the mosquitoes population size at the stationary state under the baseline usage of I₁ only at the constant exposure rate p₀. For example, a relative gain $R_{gain}(D_{epS}, p_0)$ of 5% indicates that strategy p reduces the mosquito population size over the $\bar{T}$ years deployment period by 5%, whereas a negative value implies an increase in population size compared to the reference strategy where only I₁ is used at a constant exposure rate p₀.

**Time to resistance emergence.** The evolutionary output considered is the emergence time, denoted $T_{emerg}$, of the adapted strain $x_2$, ie., the resistance level that maximizes the fitness function $\Gamma_2$ under the selective pressure of insecticide I₂. Formally, $T_{emerg}$ is defined as: (i) the first time at which the proportion of mosquitoes with resistance level $x_2$ (denoted as $Q_{x_2}$) exceeds 10% of the total population, and (ii) remains persistently higher than that of the previously dominant strain $\bar{x}_1$. More precisely, $T_{emerg}$ is the time from which, for all $t \ge T_{emerg}$, we have:

$$Q_{x_2}(t) := \frac{\sum_{j=0}^{2} \int_{0}^{\infty} A_j(t, a, x_2) \mathrm{d}a}{\sum_{j=0}^{2} \int_{0}^{\infty} (A_j(t, a, \bar{x}_0) + A_j(t, a, x_2)) \mathrm{d}a} > 10\%.$$

The first condition (i) signals the initial emergence of strain $x_2$ from rarity to a significant frequency threshold. This threshold is chosen as a practical indicator that strain $x_2$ is no longer negligible in the population. The second condition (ii) ensures that this strain $x_2$ is not only established but that its increase is sustained over time.

# 5 Results

## 5.1 Typical dynamics simulated with the model

Introducing insecticide I₂ into a landscape initially composed solely of insecticide I₁ at a proportion $p_0 = 0.6$ results in a mixed insecticide environment with I₁ and I₂. The proportions of I₁ and I₂ are updated at regular intervals of duration $D_{ur} = 3$ years over a total deployment period of $\bar{T} = 12$ years. We consider a deployment strategy consisting of: $p_1 = (0.24, 0.2, 0.16, 0.1)$ and $p_2 = (0.48, 0.6, 0.72, 0.8)$, corresponding to the proportions of I₁ and I₂, respectively, during each update period. In this new context (at least during the first three years), the fitness function $\Theta_{[p_1, p_2]}$ attains its maximum at a dominant phenotype denoted by $\bar{x}_2$ (Figure 4A). Under this deployment scenario, the relative gain achieved by introducing I₂ is approximately $R_{gain}(D_{epS}, p_0) \approx 38\%$, while the time to resistance emergence is estimated
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Evolutionary dynamics with a deployment strategy $p_1$ = (0.24, 0.2, 0.16, 0.1), $p_2$ = (0.48, 0.6, 0.72, 0.8). Here, we have $p_0$ = 0.6, $r_{\text{eff}}$ = 0.75, and $\text{var}_0$ = 0.002. (A) The fitness functions. (B) Dynamics of total population and population with phenotypic value $x_2$. The vertical dot line represents the time of resistance emergence $T_{\text{emerg}}$. (C-E) The evolutionary dynamics for unexposed and exposed adult mosquitoes</figcaption></figure>

at $T_{\text{emerg}} \approx 4.9$ years (Figure 4B). These quantities are captured through the total mosquito population size and the proportion of mosquitoes with the phenotypic value $x_2$, which maximizes the fitness function $\Gamma_2$ under the selective pressure of insecticide $I_2$ (Figure 4B). The model captures the evolutionary dynamics of the mosquito population following the introduction of the insecticide regimen by simultaneously tracking changes in population size and resistance trait distribution (Figures 4C-E). Initially, i.e., during the first cycle corresponding to the first 3 years of deployment of insecticide $I_2$ at a proportion of 0.48, the regimen induces a sharp decline in the mosquito population size, followed by a great rebound in this cycle (Figure 4B). In the second and third cycles, increasing the proportion of $I_2$ to 0.6 and 0.72, respectively, results in a slight decline in mosquito population size over a relatively short period, followed by a rapid rebound (Figure 4B). In the final cycle, increasing the proportion of $I_2$ to 0.8 leads to a slight and temporary decline in mosquito population size. However, the subpopulation carrying the phenotypic trait $x_2$, associated with resistance to insecticide $I_2$, begins to expand one year after the beginning of the second cycle. Ultimately, as this resistant phenotype becomes dominant, the mosquito population escapes control under the current composition of the insecticide landscape, i.e., under the current deployment strategy [$p_1, p_2$] (Figure 4B).

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 5 Global sensitivity analysis.** Sensitivity of the relative gain ($R_{gain}$) and the time to resistance emergence ($T_{emerg}$) to key parameters: initial insecticide $I_1$ exposure rate before the deployment of $I_2$ ($p_0$), landscape composition defined by the proportion $p_1$ of insecticide $I_1$, relative effectiveness of insecticide $I_2$ compared to $I_1$ ($r_{eff}$), the mutational variance in the unexposed adult population ($var_0$), and the mutational variance ratio (var). The hatched portions of the bars represent the main indices (the effect of each factor alone), while the complete bars, including both the hatched and unshaded parts, represent the total effect of each factor, accounting for interactions with all other factors</figcaption></figure>

## 5.2 Sensitivity analysis

We assess the sensitivity of the relative gain ($R_{gain}$) and the time to resistance emergence ($T_{emerg}$) to five key parameters: the insecticide exposure rate prior to the deployment strategy ($p_0$), the insecticide landscape composition defined by the proportion $p_1$ of insecticide $I_1$, the relative effectiveness of insecticide $I_2$ compared to $I_1$ ($r_{eff}$), the mutational variance in the unexposed adult population ($var_0$), and the mutational variance ratio (var). The parameter ranges explored are listed in Table 1.

Sensitivity analysis identified the insecticide landscape composition and the relative effectiveness of insecticide $I_2$ compared to $I_1$ as the primary drivers of the relative gain $R_{gain}$, accounting for 29% and 27% of its variance, respectively (Figure 5). These are followed by the mutational variance in the unexposed AFM population (12%) and the insecticide exposure rate (11%). In contrast, the main determinant of the time to resistance emergence $T_{emerg}$ is the mutational variance in the unexposed AFM population, explaining 31% of its variance, followed by the insecticide landscape composition (19%) and the insecticide exposure rate (17%). The mutational variance ratio has a negligible influence on both $R_{gain}$ and $T_{emerg}$, contributing less than 0.4% to the variation in either outcome. Consequently, a potential mutagenesis component—captured by a possible mutational difference between unexposed and insecticide-exposed mosquitoes—plays only a marginal overall role.

While the relative effectiveness of insecticide $I_2$ compared to $I_1$ ($r_{eff}$) accounts for 27% of the variation in the relative gain $R_{gain}$, it explains only 0.4% of the variation in the emergence time $T_{emerg}$ (Figure 5). This result highlights a trade-off between achieving a significant reduction in mosquito population size through insecticide pressure and ensuring the long-term sustainability of insecticide effectiveness by limiting the pace of resistance emergence within mosquito populations.
## 5.3 Optimal deployment strategies

We determined the optimal deployment strategy (Mos*) through a constrained nonlinear multivariable optimization, implemented with the fmincon function in MATLAB (version 2023a). Under the same parameter configuration as in Figure 4, the optimal solution prescribes proportions of $I_1$ and $I_2$ over successive 3-year update periods as $p_1 \approx (0.0023, 0.0024, 0.0027, 0.0031)$ and $p_2 \approx (0.4898, 0.9931, 0.9943, 0.9938)$, respectively. This outcome indicates that the Mos* strategy essentially phases out the inefficient insecticide $I_1$ while progressively introducing the more effective insecticide $I_2$. Implementing this strategy yields a relative gain of approximately $R_{gain}(D_{ep}S, p_0) \approx 41.2\%$ (Figure 6B), corresponding approximately to a two-fifths reduction in mosquito population size. The time to resistance emergence is estimated at $\text{Temerg} \approx 3.9$ years (Figure 6B).

The evolutionary dynamics reveal that during the first 3 years of $I_2$ deployment, mosquito density undergoes a sharp decline followed by a rapid rebound (Figure 6B). In the second cycle, when the proportion of $I_2$ is increased to 99%, a slight temporary decrease in mosquito population size is again observed, quickly followed by resurgence. During the third and fourth cycles, maintaining $I_2$ at 99% results in sustained population growth. This pattern arises because the subpopulation carrying trait $x_2$, associated with resistance to $I_2$, begins to expand the first year of the second cycle and becomes dominant thereafter, ultimately driving the failure of the vector control strategy under the current insecticide landscape (Figure 6B-E).

Overall, the optimal mosaic strategy (Mos*) fully excluded insecticide $I_1$ regardless of the mutational variance within unexposed mosquito populations ($var_0$), the relative efficacy of $I_2$ compared to $I_1$ (reff) or the initial exposure rate to $I_1$ ($p_0$) prior to $I_2$ deployment (Figure 7). This means that, under the hypothesis that the first insecticide is ineffective against mosquitoes, simultaneous use of both insecticides is rarely optimal.

For high mutational variance within unexposed mosquito populations, the optimal mosaic strategy (Mos*) consists exclusively of insecticide $I_2$ at full (100%) coverage, regardless of the values of reff or $p_0$ (Figure 7A-C). By contrast, when $var_0$ is relatively small (Figure 7D-F), insecticide $I_1$ is fully excluded from the optimal strategy. However, the proportion of $I_2$ during the first deployment period is only moderate, except when reff is low and $p_0$ is high. During the second deployment period, the proportion of $I_2$ rapidly approaches 100% once $p_0$ exceeds a certain threshold, and from the third period onward, $I_2$ is maintained at nearly 100% coverage for all values of $p_0$. Finally, for very small values of $var_0$ (Figure 7G-I), the qualitative pattern is similar to that observed in panels D-F. However, the threshold value of $p_0$ at which $I_2$ reaches full deployment in the second period is shifted to higher values. Moreover, in this regime, the proportion of $I_2$ during the third deployment period does not consistently remain close to 100%, particularly when reff is high and $p_0$ is moderate (Figure 7H-I).
<figure><img src="image_11.png" /><figcaption>**Fig. 6** Evolutionary dynamics with the optimal deployment strategy for a prior exposure rate of $p_0 = 0.6$, reff = 0.75, and var0 = 0.002. We established that the optimal mosaic strategy is obtained for $p_1 \approx (0.0023, 0.0024, 0.0027, 0.0031)$ and $p_2 \approx (0.4898, 0.9931, 0.9943, 0.9938)$. (A) the fitness functions. (B) Dynamics of total population and population with phenotypic value $x_2$. The vertical dot line represents the time of resistance emergence $T_{emerg}$. (C-E) The evolutionary dynamics for unexposed and exposed adult mosquitoes</figcaption></figure>

## 5.4 Influence of key parameters on the effectiveness of deployment strategies

Regardless of the deployment strategy— whether the simple mosaic (Mos), the simple rotation (Rot), or the optimal mosaic (Mos*)— the mutational variance (var0) and the relative effectiveness of insecticide $I_2$ compared to $I_1$ (reff) have consistent effects on outcomes: the relative gain $R_{gain}$ systematically decreases as the initial insecticide exposure rate ($p_0$), prior to the introduction of $I_2$, increases (Figure 8). This highlights that strong insecticide pressure before the deployment phase diminishes the achievable gains in mosquito population suppression. Furthermore, the Mos strategy consistently outperforms Rot, particularly when var0 is low and reff is high (Figure 8G-I). This result is intuitive, as the Rot strategy requires periodically reintroducing the already inefficient insecticide $I_1$ at 100% coverage during certain update periods, an approach that is inherently inconsistent with maintaining effective population suppression. Finally, the Mos* strategy provides a marked performance advantage over the simple Mos strategy when var0 is low, and this advantage remains largely robust to variations in reff (Figure 8G-I).
<figure><img src="image_12.png" /><figcaption>**Fig. 7** The insecticide landscape with the mosaic optimal strategy. The mutational variance ($var_0$) and the relative effectiveness (reff) of insecticide $I_2$ compared to $I_1$ are variables. With $n \in \{1, \dots, 4\}$, $p_1^n$ and $p_2^n$ correspond to the proportions of insecticides $I_1$ and $I_2$ during the $n$-th updating period</figcaption></figure>

# 6 Conclusion and discussion

The extensive use of pyrethroid-based insecticides has led to the emergence and pro-liferation of insecticide resistance (IR) in mosquito populations, posing a significant threat to malaria control and elimination efforts, which rely heavily on vector control through insecticide interventions. Given the limited arsenal of insecticides available for public health applications, there is an urgent need to develop and implement effec-tive insecticide resistance management (IRM) strategies to ensure the sustainable use of these compounds. Building upon previous research (Kezeta-Bondja et al. 2025), we model insecticide resistance as a quantitative trait, reflecting its polygenic nature rather than a monogenic basis. This approach enabled a more accurate representation of the transient evolutionary dynamics associated with resistance development. We integrate this trait-based resistance model into an age-structured framework governed by integro-differential equations, wherein the continuous resistance trait modulates key life-history parameters, including oviposition and mortality. We assess the effi-cacy of three insecticides deployment strategies –simple mosaic (Mos), simple rotation (Rot), and optimal mosaic (Mos*)– in their ability of reducing mosquito population size, thereby contributing novel insights to the existing body of literature on IRM.
<figure><img src="image_13.png" /><figcaption>**Fig. 8** The relative gain with strategies Rot, Mos, and Mos*. The mutational variance (var₀) and the relative effectiveness (reff) of insecticide I₂ compared to I₁ are variables</figcaption></figure>

The model’s mathematical properties were rigorously analyzed, including the existence of a unique maximal bounded semiflow and conditions for the existence of steady states. Through parameterization and numerical simulations, we explored both epidemiological outcomes— quantified by the relative gain of introducing insecticide I₂ into the insecticides landscape— and evolutionary dynamics, characterized by the time of resistance emergence. Comparative analysis of the three strategies revealed their differential impacts on mosquito population reduction following the deployment of I₂. These findings offer valuable guidance for the design of sustainable and adaptive vector control strategies.

While the Mos consistently outperforms the Rot strategy, regardless of the mutational variance (var₀) and the relative effectiveness of insecticide I₂ compared to I₁ (r<sub>eff</sub>) (Figure 8), the difference in performance between Mos and Rot strategies – denoted by Δ(Mos, Rot) = R<sub>gain</sub>(Mos) – R<sub>gain</sub>(Rot) – is influenced by the insecticide landscape updating interval (D<sub>ur</sub>). Specifically, values of D<sub>ur</sub> ranging from 1 to 3 years have a negligible impact on the performance difference Δ(Mos, Rot). In contrast, when D<sub>ur</sub> = 4 years, the performance advantage of the Mos strategy over Rot is reduced with Mos still performing at least as well as Rot overall (Figure 9). This result is intuitive, since the Rot strategy requires periodically reintroducing the already
ineffective insecticide $I_1$ at full (100%) coverage during certain update periods. This structural feature explains why Rot consistently underperforms compared to Mos in our framework.

This configuration stands in contrast to those of Levick et al. (2017), who found that sequential deployment (Rot) outperformed mixtures (Mos) in some configurations. The apparent discrepancy arises from at least two main differences in initial assumptions. The first key distinction from Levick et al. (2017) lies in the initial efficacy assumptions regarding insecticides at the time of deployment. While authors in Levick et al. (2017) considered scenarios in which both insecticides retained effectiveness against mosquito populations, our framework assumes that insecticide $I_1$ has already lost efficacy due to widespread resistance, leaving $I_2$ as the sole insecticide with residual activity at deployment. The second difference pertains to the representation of the resistance trait: here, it is modeled as a scalar continuous variable $x \in \mathbb{R}$, rather than as a two-dimensional vector $x = (x_1, x_2) \in \mathbb{R}^2$, where $x_1$ and $x_2$ correspond to resistance levels against $I_1$ and $I_2$, respectively. This simplification reflects the notion that multi-insecticide resistance phenotypes often collapse onto low-dimensional trade-off axes or manifolds, a perspective that is plausible under at least two scenarios: when $I_1$ and $I_2$ represent different formulations of the same chemical compound, or when the resistance trait encapsulates a generalized mechanism effective against both insecticides.

Although evidence supporting such dimensionality reduction exists in the context of multidrug resistance (e.g., Barbosa et al. (2019); Imamovic and Sommer (2013); Lázár et al. (2014); Podnecky et al. (2018); Sakenova et al. (2024)), this may not universally apply to multi-insecticide resistance. Importantly, extending the trait space to $x \in \mathbb{R}^n$ does not fundamentally alter the main theoretical results established here, as demonstrated in Djidjou-Demasse et al. (2017). However, a significant challenge arises in specifying the age- and trait-dependent eggs-laying rate function $r(a, x)$. For instance, when considering two insecticides, this function generalizes to $r(a, x_1, x_2)$, and formulating an explicit expression without imposing a collateral insecticide resistance score $s$, which maps the multidimensional trait $(x_1, x_2) \in \mathbb{R}^2$ onto a unidimensional axis $s(x_1, x_2) \in \mathbb{R}$, remains a critical obstacle when modeling resistance as a multidimensional trait.

More than two insecticides may be available at the time of deployment. The model developed here naturally accommodates the introduction of additional insecticides into the landscape. For example, in Appendix A, we consider the case where two insecticides, $I_2$ and $I_3$, are introduced after a prolonged period of exclusive use of $I_1$. Assuming insecticide $I_3$ is more effective than $I_2$ at the time of deployment, the optimal mosaic strategy (Mos*) consists primarily of the exclusive use of $I_3$, followed by its gradual incorporation into the insecticide landscape (Figure S1). This finding suggests that the simultaneous deployment of two insecticides— differing only in their relative efficacy against the target vector population— may not represent the most effective strategy. Instead, a sequential approach, in which a first insecticide is used optimally until resistance reaches a critical threshold, followed by the introduction of the second insecticide, appears to yield greater benefits.

We also compared the simple mosaic (Mos) and simple rotation (Rot) strategies using insecticides $I_2$ and $I_3$. When the Rot strategy begins with $I_2$, and for high muta-
tional variance (var₀), Mos consistently outperforms Rot regardless of the mutational variance ratio (var), with the performance gap widening as the insecticide update interval (Dur) increases (Figure S2A, D). For var₀ = 0.002, Mos performs worse than Rot when Dur = 2 or 3 years, but better when D ur = 1 year (Figure S2B, E). A similar qualitative pattern is observed for var₀ = 0.001, though with larger performance differences (Figure S2C, F). Notably, increasing D ur slightly amplifies the performance gap (Δ(Mos, Rot)) when var₀ = 0.01, but reduces it at lower mutational variances. When the Rot strategy begins with I₃ (Figure S3), Mos and Rot achieve nearly identical relative gains at the highest mutational variance (var₀ = 0.01), irrespective of D ur or var. For lower var₀ values, however, Mos consistently outperforms Rot, and increasing D ur further amplifies the performance gap (Δ(Mos, Rot)) across all values of var.

Our results also reinforce a well-established principle in the epi-evolutionary man-
agement of insecticide resistance: the optimal dose for slowing resistance evolution
often corresponds to the lowest dose that remains epidemiological effective. Simila-
lar patterns have been documented in antimicrobial stewardship, pesticide use, and
pathogen evolution, where overly aggressive interventions intensify selection and
accelerate the spread of resistance (e.g., Blanquart et al. (2017); Day and Read (2016);
South and Hastings (2018)). In practice, however, deliberately “holding back” from
maximal pressure is difficult to implement. Such strategies may conflict with short-
term control objectives, require coordination across programs, and rely on robust
surveillance and compliance. Our findings therefore highlight a persistent tension:
although evolutionary theory shows that moderated intervention can outperform max-
imal pressure in the long run, translating this insight into operational policy remains
challenging. This underscores the need for decision-support tools and governance
frameworks that jointly optimize immediate disease-control goals and long-term resis-
tance management.

Although the model developed here could be extended to explicitly incorporate additional mosquito life stages (e.g., larval and pupal stages) and thereby accommodate interventions such as larviciding or pupaciding, our analysis focuses on control measures targeting adult mosquitoes, notably insecticide-treated nets (ITNs) and indoor residual spraying (IRS). According to WHO guidelines, ITNs and IRS are classified as core malaria vector-control tools because they provide strong individual protection and substantially reduce adult mosquito survival, generating major community-level impact (Guidelines for Malaria Vector Control 2019). By contrast, larval source management (including larvicides and pupacides) is considered only a supplementary intervention, recommended solely in settings where larval habitats are few, fixed, and easily identifiable, and not as a replacement for ITNs or IRS. This reflects WHO's conclusion that larviciding produces weaker and less consistent effects on malaria transmission (Guidelines for Malaria Vector Control 2019).

Another limitation of the model proposed here is the absence of climatic
drivers–such as temperature, humidity, rainfall, and wind–which are known to
strongly influence mosquito life-cycle traits and, consequently, population abundance,
e.g., Abdelrazec and Gumel (2017); Djidjou-Demasse et al. (2020). Incorporating
such environmental forcing would require coupling non-local interactions with a
non-autonomous dynamical system, which raises substantial analytical challenges.
Addressing this limitation therefore constitutes an important but technically demanding perspective for future work.

## A Model with 3 insecticides

Here we assume that after the failure of $I_1$, we introduced two insecticides, denoted $I_2$ and $I_3$. Insecticide $I_2$ is characterized by a relative effectiveness compared to $I_1$ with values ranging in {0.5, 0.6, 0.7, 0.8}, while insecticide $I_3$ has a relative effectiveness compared to insecticide $I_1$ fixed to 0.9. During the $T$ years period of insecticide deployment, the insecticide landscape composed of $I_1, I_2$, and $I_3$, is revised $N$ times at times $T_0, T_1, \dots, T_{N-1}$, at regular intervals of length $D_{ur}$. Moreover, during each $n^{th}$ period, the insecticide landscape composition is defined by the proportions $p_1^n, p_2^n$ and $p_3^n$, representing the proportion of insecticides $I_1, I_2$ and $I_3$, respectively, with $p_1^n + p_2^n + p_3^n \le 1$. The model with three insecticides is given for $n \in \{1, \dots, N\}$ and $t \in [T_{n-1}, T_n)$, by:

<figure><img src="image_14.png" /><figcaption>**Fig. 9** The difference in performance between Mos and Rot strategies. Here $\Delta$(Mos, Rot) = $R_{gain}$(Mos) - $R_{gain}$(Rot). The mutational variance ($var_0$) and the relative effectiveness (reff) of insecticide $I_2$ compared to $I_1$ are variables.</figcaption></figure>
$$
\left\{
\begin{aligned}
A_0^n(t, a = 0, x) &= (1 - p_1^n - p_2^n - p_3^n) \tau(x) \sum_{j=0}^{3} \gamma_j(x) E_j^n(t, x), \\
A_1^n(t, a = 0, x) &= p_1^n \tau(x) \sum_{j=0}^{3} \gamma_j(x) E_j^n(t, x), \\
A_2^n(t, a = 0, x) &= p_2^n \tau(x) \sum_{j=0}^{3} \gamma_j(x) E_j^n(t, x), \\
A_3^n(t, a = 0, x) &= p_3^n \tau(x) \sum_{j=0}^{3} \gamma_j(x) E_j^n(t, x),
\end{aligned}
\right.
$$

and for $j \in \{0, 1, 2, 3\}$:

$$
\left\{
\begin{aligned}
\partial_t E_j^n(t, x) &= H(E_j^n(t)) \int_{\Omega} \int_0^{\infty} m_j(x, y) r(a, y) A_j^n(t, a, y) \, \mathrm{d}ady - (\mu_j(x) + \gamma_j(x)) E_j^n(t, x), \\
(\partial_t + \partial_a) A_j^n(t, a, x) &= -d_j(a, x) A_j^n(t, a, x),
\end{aligned}
\right.
$$

together with the initial conditions $E_j^1(0, x) = E_{j,0}(x)$, $A_j^1(0, a, x) = A_{j,0}(a, x)$ and $E_j^n(T_{n-1}, x) = E_j^{n-1}(T_{n-1}, x)$, $A_j^n(T_{n-1}, a, x) = A_j^{n-1}(T_{n-1}, a, x)$ when $n \in \{2, 3, \dots, N\}$.

Here, $E_j^n$ represents the total number of egg laid by AFM in $j$-insecticide state and $A_j^n$ denotes the total number of adult mosquitoes in $j$-insecticide state, at time period $n$ and for $j \in \{0, 1, 2, 3\}$. The rest of parameters are unchanged.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00285-026-02343-z.

**Acknowledgements** This work received support from the project ACoMVeC (INV-047049). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Data availability** All the data used are included in this paper.

## Declarations

**Conflicts of interest** The authors of this article declare that they have no financial conflict of interest with the content of this article.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
## References

* Guidelines for Malaria Vector Control. World Health Organization, Geneva, (2019)

* Abdelrazec A, Gumel AB (2017) Mathematical assessment of the role of temperature and rainfall on mosquito population dynamics. J Math Biol 74(6):1351–1395

* Alout H, Roche B, Dabiré RK, Cohuet A (2017) Consequences of insecticide resistance on malaria transmission. PLoS Pathog 13(9):e1006499

* Barbosa C, Römhild R, Rosenstiel P, Schulenburg H (2019) Evolutionary stability of collateral sensitivity to antibiotics in the model pathogen Pseudomonas aeruginosa. Elife 8:e51481

* Bayoh M, Lindsay S (2003) Effect of temperature on the development of the aquatic stages of Anopheles gambiae sensu stricto (Diptera: Culicidae). Bull Entomol Res 93(5):375–381

* Bhatt S, Weiss DJ, Cameron E, Bisanzio D, Mappin B, Dalrymple U, Battle KE, Moyes CL, Henry A, Eckhoff PA, Wenger EA, Briët O, Penny MA, Smith TA, Bennett A, Yukich J, Eisele TP, Griffin JT, Fergus CA, Lynch M, Lindgren F, Cohen JM, Murray CLJ, Smith DL, Hay SI, Cibulskis RE, Gething PW (2015) The effect of malaria control on Plasmodium falciparum in Africa between 2000 and 2015. Nature 526(7572):207–211

* Birget PLG, Koella JC (2015) A genetic model of the effects of insecticide-treated bed nets on the evolution of insecticide-resistance. Evolution Medicine and Public Health 2015(1):205–215

* Blanquart F, Lehtinen S, Fraser C (2017) An evolutionary model to predict the frequency of antibiotic resistance under seasonal antibiotic use, and an application to Streptococcus pneumoniae. Proceedings of the Royal Society B Biological Sciences 284(1855):20170679

* Carnevale P, Robert V, editors. (2009) Les anophèles: Biologie, transmission du Plasmodium et lutte antivectorielle. IRD Éditions

* CDC. Life Cycle of Anopheles Mosquitoes, May (2024)

* Churcher TS, Lissenden N, Griffin JT, Worrall E, Ranson H (2016) The impact of pyrethroid resistance on the efficacy and effectiveness of bednets for malaria control in Africa. Elife 5:e16090

* Day T, Read AF (2016) Does High-Dose Antimicrobial Chemotherapy Prevent the Evolution of Resistance? PLoS Comput Biol 12(1):e1004689

* Djidjou-Demasse R, Abiodun GJ, Adeola AM, Botai JO (2020) Development and analysis of a malaria transmission mathematical model with seasonal mosquito life-history traits. Stud Appl Math 144(4):389–411

* Djidjou-Demasse R, Ducrot A, Fabre F (2017) Steady state concentration for a phenotypic structured problem modeling the evolutionary epidemiology of spore producing pathogens. Math Models Methods Appl Sci 27(02):385–426

* Djidjou-Demasse R, Sofonea MT, Choisy M, Alizon S (2023) Within-host evolutionary dynamics of antimicrobial quantitative resistance. Mathematical Modelling of Natural Phenomena 18:24

* Gabrielkuniyoshi ML, Piodossantos FL (2017) Mathematical modelling of vector-borne diseases and insecticide resistance evolution. J Venomous Anim Toxins Incl Trop Dis 23(1):34

* Hale JK (2014) Asymptotic Behavior of Dissipative Systems. Number v.25 in Mathematical Surveys and Monographs. American Mathematical Society, Providence

* Helps JC, Paveley ND, White S, Van Den Bosch F (2020) Determinants of optimal insecticide resistance management strategies. J Theor Biol 503:110383

* Hemingway J, Ranson H, Magill A, Kolaczinski J, Fornadel C, Gimnig J, Coetzee M, Simard F, Roch DK, Hinzoumbe CK, Pickett J, Schellenberg D, Gething P, Hoppé M, Hamon N (2016) Averting a malaria disaster: Will insecticide resistance derail malaria control? The Lancet 387(10029):1785–1788

* Hobbs NP, Hastings I (2025) The impact of insecticide decay on the rate of insecticide resistance evolution for monotherapies and mixtures. Malar J 24(1):50

* Hobbs NP, Hastings I (2025) Simulating dynamic insecticide selection pressures for resistance management in mosquitoes assuming polygenic resistance. PLoS Comput Biol 21(4):e1012944

* Hobbs NP, Weetman D, Hastings IM (2023) Insecticide resistance management strategies for public health control of mosquitoes exhibiting polygenic resistance: A comparison of sequences, rotations, and mixtures. Evol Appl 16(4):936–959

* Iannelli M (1995) Mathematical Theory of Age-Structured Population Dynamics. Giardini editori e stampatori, Pisa

* Imamovic L, Sommer MOA (2013) Use of Collateral Sensitivity Networks to Design Drug Cycling Protocols That Avoid Resistance Development. Sci Transl Med 5(204) Sept
Kezeta-Bondja S, Mann-Manyombe ML, Tewa J-J, Wondji CS, Djidjou-Demasse R (2025) Modeling the evolutionary dynamics of insecticide quantitative resistance in mosquitoes populations Jan

Killeen GF, Ranson H (2018) Insecticide-resistant malaria vectors must be tackled. The Lancet 391(10130):1551–1552

Lázár V, Nagy I, Spohn R, Csörgő B, Györkei Á, Nyerges Á, Horváth B, Vörös A, Busa-Fekete R, Hrtyan M, Bogos B, Méhi O, Fekete G, Szappanos B, Kégl B, Papp B, Pál C (2014) Genome-wide analysis captures the determinants of the antibiotic cross-resistance interaction network. Nat Commun 5(1):4352

Levick B, South A, Hastings IM (2017) A Two-Locus Model of the Evolution of Insecticide Resistance to Inform and Optimise Public Health Insecticide Deployment Strategies. PLoS Comput Biol 13(1):e1005327

Madgwick PG, Kanitz R (2022) Modelling new insecticide-treated bed nets for malaria-vector control: How to strategically manage resistance? Malar J 21(1):102

Madgwick PG, Kanitz R (2024) What is the value of rotations to insecticide resistance management? Pest Manag Sci 80(4):1671–1680

Magal P, Ruan S (2018) Theory and Applications of Abstract Semilinear Cauchy Problems, vol 201. Applied Mathematical Sciences. Springer International Publishing, Cham

Mohammed-Awel J, Gumel AB (2019) Mathematics of an epidemiology-genetics model for assessing the role of insecticides resistance on malaria transmission dynamics. Math Biosci 312:33–49

Mohammed-Awel J, Gumel AB (2023) Can insecticide resistance increase malaria transmission? A genetics-epidemiology mathematical modeling approach. J Math Biol 87(2):28

Opondo KO, Alyko E, Smith S, Levine R, Donnelly J, Hughes T, Schnabel D, Jose R, Jalloh A, Lolleh U, Conteh S, Iyikirenga L, Jacob D, Mokuena M, Yamba F, Yihdego Y (2025) Characterization of insecticide resistance mechanisms in the Anopheles gambiae population of Sierra Leone. Malar J 24(1):80

Podnecky NL, Fredheim EGA, Kloos J, Sørum V, Primicerio R, Roberts AP, Rozen DE, Samuelsen Ø, Johnsen PJ (2018) Conserved collateral antibiotic susceptibility networks in diverse clinical strains of Escherichia coli. Nat Commun 9(1):3673

Ranson H, Lissenden N (2016) Insecticide Resistance in African Anopheles Mosquitoes: A Worsening Situation that Needs Urgent Action to Maintain Malaria Control. Trends Parasitol 32(3):187–196

Sakenova N, Cacace E, Orakov A, Huber F, Varik V, Kritikos G, Michiels J, Bork P, Cossart P, Goemans CV, Typas A (2024) Systematic mapping of antibiotic cross-resistance and collateral sensitivity with chemical genetics. Nat Microbiol 10(1):202–216

South A, Hastings IM (2018) Insecticide resistance evolution with mixtures and sequences: A model-based explanation. Malar J 17(1):80

Strode C, Donegan S, Garner P, Enayati AA, Hemingway J (2014) The Impact of Pyrethroid Resistance on the Efficacy of Insecticide-Treated Bed Nets against African Anopheline Mosquitoes: Systematic Review and Meta-Analysis. PLoS Med 11(3):e1001619

Sudo M, Takahashi D, Andow DA, Suzuki Y, Yamanaka T (2018) Optimal management strategy of insecticide resistance under various insect life histories: Heterogeneous timing of selection and interpatch dispersal. Evol Appl 11(2):271–283

Tatchou-Nebangwa NMT, Mugenzi LMJ, Muhammad A, Nebangwa DN, Kouamo MFM, Tagne CSD, Tekoh TA, Tchouakui M, Ghogomu SM, Ibrahim SS, Wondji CS (2024) Two highly selected mutations in the tandemly duplicated CYP6P4a and CYP6P4b genes drive pyrethroid resistance in Anopheles funestus in West Africa. BMC Biol 22(1):286

Thieme HR (1990) Semiflows generated by Lipschitz perturbations of non-densely defined operators. Differ Integr Equ 3(6) Jan

Thieme HR (2011) Global stability of the endemic equilibrium in infinite dimension: Lyapunov functions and positive operators. J Differ Equ 250(9):3772–3801

Webb GF (1985) Theory of Nonlinear Age-Dependent Population Dynamics. Number 89 in Pure and Applied Mathematics. Dekker M (ed.) New-York Basel

Weedall GD, Riveron JM, Hearn J, Irving H, Kamdem C, Fouet C, White BJ, Wondji CS (2020) An Africa-wide genomic evolution of insecticide resistance in the malaria vector Anopheles funestus involves selective sweeps, copy number variations, gene conversion and transposons. PLoS Genet 16(6):e1008822

WHO (2024) World Malaria Report 2024. World Health Organization, Geneva, 1st edition edition
World Health Organization (2018) Global Report on Insecticide Resistance in Malaria Vectors: 2010–2016. World Health Organization, Geneva

World Health Organization and UNICEF/UNDP/World Bank/WHO Special Programme for Research and Training in Tropical Diseases (2017) *Global Vector Control Response 2017-2030*. World Health Organization, Geneva

World Health Organization (2012) Global Malaria Programme. Global plan for insecticide resistance management in malaria vectors. page 131

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.