# Optimizing MDA and antimalarial treatment in the presence of drug resistance for effective malaria control

Manuela M Nimpa, Hyacinthe N Teytsa, Joseph Mbang, Charles S Wondji, Ramsès Djidjou-Demasse

► To cite this version:

Manuela M Nimpa, Hyacinthe N Teytsa, Joseph Mbang, Charles S Wondji, Ramsès Djidjou-Demasse. Optimizing MDA and antimalarial treatment in the presence of drug resistance for effective malaria control. Mathematical Biosciences and Engineering, 2025, 22 (8), pp.1898-1930. 10.3934/mbe.2025069. hal-05121237

HAL Id: hal-05121237

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

Submitted on 19 Jun 2025

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 Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Research article

# Optimizing MDA and antimalarial treatment in the presence of drug resistance for effective malaria control

Manuela M. Nimpa<fnref n="1" />,<fnref n="2" />, Hyacinthe N. Teytsa<fnref n="4" />,<fnref n="5" />, Joseph Mbang<fnref n="2" />, Charles S. Wondji<fnref n="1" />,<fnref n="3" /> and Ramsès Djidjou-Demasse<fnref n="4" />,<fnref n="5" />,*

**Abstract:** Antimalarial drugs are critical for controlling malaria, but the emergence of drug resistance poses a significant challenge to global eradication efforts. This study explores strategies to minimize resistance prevalence and improve malaria control, particularly through the use of mass drug administration (MDA) in combination with antimalarial drugs. We develop a compartmental mathematical model that incorporates asymptomatic, paucisymptomatic, and clinical states of infection and evaluates the impact of resistance mutations on transmission dynamics. The model includes both treated and untreated states among infected and recovered individuals, with a focus on optimizing control strategies through MDA and antimalarial treatment. A global sensitivity analysis identifies the critical factors that influence malaria dynamics, including MDA coverage, treatment access for different infection states, the probability of mutation from treated sensitive human infections, to treated resistant human infections and the initial prevalence of resistance. The model is extended to include optimal control strategies that provide time-dependent control interventions for treatment and MDA. Intuitively, when the mutation rate is relatively low, the optimal strategy combines the use of antimalarial drugs and MDA, with a gradual decrease in antimalarial drug use over time, ensuring sustainable malaria control. In contrast, at higher mutation rates, the strategy prioritizes broader deployment of MDA while significantly reducing reliance on antimalarial to minimize the risk of resistance developing. Numerical simulations of the optimal control problem reinforce the importance of strategic intervention in mitigating drug resistance. This study contributes to understanding the role of MDA and treatment strategies in the control of malaria, with implications for optimizing malaria control programs in endemic regions.

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

<footnote n="2">² Department of Mathematics, Faculty of Sciences, University of Yaoundé I, 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</footnote>

<footnote>* Correspondence: Email: ramses.djidjoudemasse@ird.fr.</footnote>
**Keywords:** nonlinear dynamical systems; malaria; drug resistance; mass drug administration; optimal control

# 1. Introduction

Antimalarial drugs are crucial for controlling and eliminating malaria worldwide. However, the emergence of resistance in malaria parasites has significantly reduced their effectiveness, posing a critical challenge to global malaria control efforts [1]. Resistance often arises from genetic mutations in *Plasmodium falciparum*, enabling the parasite to withstand specific treatments. The spread and impact of these mutations on malaria endemicity have been extensively documented [1–6]. These developments raise concerns about the continued efficacy of key interventions, such as intermittent preventive treatment in pregnancy with sulfadoxine-pyrimethamine and seasonal malaria chemoprevention (SMC), particularly in affected regions. Consequently, it is essential to explore strategies to minimize the prevalence of resistance and reduce malaria-related mortality. One promising approach endorsed by the World Health Organization (WHO) is the mass drug administration (MDA) [7].

MDA involves administering antimalarial drugs to all individuals in a defined geographic area, regardless of the presence of symptoms or infection [7]. When implemented effectively, MDA can eliminate both drug-resistant and drug-sensitive parasites, providing the long-term benefit of prolonging the therapeutic efficacy of first-line treatments across entire populations, not just within the targeted region [8]. Given these advantages, incorporating MDA into models of malaria transmission dynamics is essential for understanding its role in controlling drug resistance and optimizing intervention strategies.

Malaria's pathogenesis involves multiple interacting factors, and significant research on malaria's transmission dynamics in the presence of drug resistance has been carried out over the past few decades. These studies include models with and without the MDA strategy, as well as models that incorporate genetic mutations in both the mosquito and human populations (e.g., [3, 9–11]). In [9], an epidemiological model for the spread of antimalarial resistance, excluding mutations, was developed, considering the treated and untreated statuses of individuals. In [11], a mathematical model for malaria was formulated with treatment and three levels of resistance in humans, where the results showed the limited benefits of treatments in a population with resistant strains, especially in high-transmission areas. The prevalence of mutations related to antimalarial resistance in *Plasmodium falciparum* was studied using a Bayesian malaria model in [3]. A simple compartmental malaria model including MDA was developed in [10] to compute the probability of extinction. The model analytically calculated the expected time to extinction in small populations following MDA for various scenarios.

While some human malaria infections initially progress through an asymptomatic phase before developing into symptomatic or clinical cases, others remain paucisymptomatic throughout the entire duration of the infection, meaning they never advance to symptomatic or clinical stages, (e.g., [12, 13]). Despite their mild symptoms, paucisymptomatic infections contribute to malaria transmission similarly to asymptomatic cases. Due to the overlap between malaria symptoms and those of other diseases, paucisymptomatic individuals often fail to realize they are infected. As a result, the parasites can persist in their bodies for extended periods, facilitating the spread of the disease [12, 13]. However,
the aforementioned studies do not incorporate the paucisymptomatic characteristic as addressed in the current work. In addition, we examine the combined effects of antimalarial drugs and the implementation of MDA within the host population for effective malaria control. Using a compartmental modeling approach, we assess the impact of resistance-associated mutations on malaria dynamics; we optimize control strategies, including MDA and the treatment proportions for individuals in asymptomatic, paucisymptomatic, and clinical states; and, finally, we investigate how the proportion of paucisymptomatic individuals affects the control strategies.

We begin by considering a model without optimal control that accounts for both the treated and untreated states among infected and recovered individuals. This approach is motivated by the challenges in accessing treatment in certain endemic regions. An in-depth asymptotic analysis of the model's steady states is performed using various threshold parameters, including the basic reproduction number, $\mathcal{T}_0$. Key findings include the local stability of the disease-free equilibrium and conditions under which the malaria model exhibits bistability through a backward bifurcation. Our findings suggest that increasing the proportion of infected individuals receiving treatment and the coverage of MDA significantly reduces the spread of resistance and malaria-related mortality. Additionally, a global sensitivity analysis is conducted to identify the key parameters influencing the prevalence of resistance and mortality. Four critical factors are identified: (i) The proportion of individuals covered by MDA, reflecting the strategy's acceptability in endemic areas; (ii) the proportion of asymptomatic individuals receiving antimalarial treatment, highlighting access to treatment for this group; (iii) the proportion of clinical cases receiving treatment; and (iv) the initial prevalence of resistance and the mutation rate.

Next, we extend the model to include optimal control, where the MDA rate and the proportions of infectious individuals receiving treatment serve as time-dependent control strategies. We analyze the impact of these control measures on the transmission dynamics of malaria associated with resistance mutations. Pontryagin's maximum principle [14] is applied to characterize the optimal control strategies. Numerical simulations of the optimal control problem reinforce the findings of the global sensitivity analysis.

The remainder of the manuscript is organized as follows: Section 2 presents the model's description and associated parameters. In Section 3, we establish some fundamental properties of the model, including the existence of a globally defined semiflow, the asymptotic stability of the model's stationary states, and an analysis of the backward bifurcation. Baseline model simulations and sensitivity analysis results are detailed in Section 4. Section 5 introduces the optimal control problem, including the derivation of the necessary optimality conditions and numerical simulations of the control strategies. The overall findings of the manuscript are discussed in Section 6. Finally, Section 7 provides the detailed proofs of the manuscript's main results.

## 2. Description of the model

We examine the interconnected dynamics of humans and mosquitoes in malaria transmission. We use the variable $k \in \{r, s\}$ to denote the type of infection, where $r$ represents a resistant infection and $s$ represents a sensitive infection. Furthermore, the variable $\ell \in \{U, T\}$ indicates the treatment status of the human host, with $U$ representing an untreated infection and $T$ representing a treated infection. The human population is classified into five states, namely susceptible ($S_h$), asymptomatic infected ($A_h^{k,\ell}$), paucisymptomatic infected ($P_h^{k,\ell}$), clinical cases ($I_h^{k,\ell}$), and recovered ($R_h^{k,\ell}$), while the mosquito
population is divided into susceptible ($S_m$) and infectious mosquitoes ($I_m^k$). The disease is transmitted from infectious mosquitoes to humans at time $t$ through the force of infection, denoted as

$$ \lambda_m^k(t) = \frac{\beta_m^k}{N_h(t)} I_m^k(t), \quad k \in \{s, r\}, $$

where $\beta_m^k$ represents the transmission probability, and $N_h$ denotes the total human population, given by

$$ N_h(t) = S_h(t) + \sum_{k \in \{r,s\}} \sum_{\ell \in \{T,U\}} (A_h^{k,\ell}(t) + P_h^{k,\ell}(t) + I_h^{k,\ell}(t) + R_h^{k,\ell}(t)). $$

The force of infections from infected humans, either with sensitive infection ($\lambda_h^s$) or with resistant infection ($\lambda_h^r$), to mosquitoes is given by

$$ \lambda_h^s(t) = \frac{\beta_h^s}{N_h(t)} [A_h^{s,U}(t) + P_h^{s,U}(t) + I_h^{s,U}(t) + \varepsilon_r (A_h^{r,T}(t) + P_h^{r,T}(t) + I_h^{r,T}(t))] \\ + \frac{\beta_h^r}{N_h(t)} [(1 - \varepsilon_s) (A_h^{s,T}(t) + P_h^{s,T}(t) + I_h^{s,T}(t))], $$

$$ \lambda_h^r(t) = \frac{\beta_h^r}{N_h(t)} [A_h^{r,U}(t) + P_h^{r,U}(t) + I_h^{r,U}(t) + \varepsilon_s (A_h^{s,T}(t) + P_h^{s,T}(t) + I_h^{s,T}(t))] \\ + \frac{\beta_h^s}{N_h(t)} [(1 - \varepsilon_r) (A_h^{r,T}(t) + P_h^{r,T}(t) + I_h^{r,T}(t))]. $$

Here, $\varepsilon_r$ represents the probability of mutation from treated resistant human infections to treated sensitive human infections, while $\varepsilon_s$ represents the probability of mutation from treated sensitive human infections to treated resistant human infections.

The model we shall consider then reads

$$
\left\{
\begin{aligned}
\dot{S}_h &= \Lambda_h - \mu_h S_h - (\lambda_m^s + \lambda_m^r) S_h + \rho_h \sum_{k \in \{r,s\}} \sum_{\ell \in \{U,T\}} R_h^{k,\ell}, \\
\dot{A}_h^{k,T} &= -\mu_h A_h^{k,T} + (1-p)\phi(t)\lambda_m^k S_h - \nu_h A_h^{k,T} - g(t)A_h^{k,T}, \quad k \in \{s, r\}, \\
\dot{A}_h^{k,U} &= -\mu_h A_h^{k,U} + (1-p)(1-\phi(t))\lambda_m^k S_h - \nu_h A_h^{k,U} - g(t)A_h^{k,U}, \quad k \in \{s, r\}, \\
\dot{P}_h^{k,T} &= -\mu_h P_h^{k,T} + p\phi(t)\lambda_m^k S_h - g(t)P_h^{k,T}, \quad k \in \{s, r\}, \\
\dot{P}_h^{k,U} &= -\mu_h P_h^{k,U} + p(1-\phi(t))\lambda_m^k S_h - g(t)P_h^{k,U}, \quad k \in \{s, r\}, \\
\dot{I}_h^{k,T} &= -\mu_h I_h^{k,T} + \theta(t)\nu_h A_h^{k,T} + \theta(t)\nu_h A_h^{k,U} - \gamma_h^{k,T} I_h^{k,T} - \delta_h^{k,T} I_h^{k,T}, \quad k \in \{s, r\}, \\
\dot{I}_h^{k,U} &= -\mu_h I_h^{k,U} + (1-\theta(t))\nu_h A_h^{k,U} + (1-\theta(t))\nu_h A_h^{k,U} - \gamma_h^{k,U} I_h^{k,U} - \delta_h^{k,U} I_h^{k,U}, \quad k \in \{s, r\}, \\
\dot{R}_h^{k,\ell} &= -(\mu_h + \rho_h)R_h^{k,\ell} + \gamma_h^{k,\ell} I_h^{k,\ell} + g(t)A_h^{k,\ell} + g(t)P_h^{k,\ell}, \quad k \in \{s, r\}, \ell \in \{U, T\},
\end{aligned}
\right.
\quad (2.1)
$$

and

$$
\left\{
\begin{aligned}
\dot{S}_m &= \Lambda_m - \mu_m S_m - (\lambda_h^s + \lambda_h^r) S_m, \\
\dot{I}_m^k &= \lambda_h^k S_m - \delta_m^k I_m^k, \quad k \in \{s, r\}.
\end{aligned}
\right.
\qquad (2.2)
$$

The dynamic of the total number of disease-induced deaths is defined by

$$ \dot{D}_h = \sum_{k \in \{s,r\}, \ell \in \{T,U\}} \delta_h^{k,\ell} I_h^{k,\ell}. $$
<figure><img src="image_3.png" /><figcaption>Figure 1. The model flow diagram.</figcaption></figure>

In System (2.1), the variable $\phi(t)$ denotes the rate of exposure of asymptomatic infections to antimalarial treatment without a curative intention at time $t$, whereas $\theta(t)$ represents the rate of exposure of clinical cases to antimalarial treatment with a curative intention. In addition, $g(t)$ describes the implementation of MDA efforts. Note that

$$g(t) = \frac{q_{\text{mda}}(t)}{D_{\text{mda}}},$$

where $q_{\text{mda}}(t)$ denotes the proportion of the target population participating in MDA, and $D_{\text{mda}}$ represents the delay required for an individual to be considered no longer a carrier after taking MDA treatment. Systems (2.1) and (2.2) will be considered under the following natural assumptions:

$$v_h \geq 0; \rho_h \geq 0; g \geq 0; p, \theta, \phi \in (0, 1); \Lambda_k \geq 0 \ (\text{for } k = h, m);$$
$$ \delta_m^k \ge 0; \delta_h^{k,l} \ge 0; \gamma_h^{k,l} \ge 0; \text{ (for } k=s,r; \text{ and } l=T,U). $$

The model is shown in Figure 1, and the remaining parameters of Models (2.1) and (2.2) are summarized in Tables 1 and 2.

Newborn humans and mosquitoes appear with the flux $\Lambda_h$ and $\Lambda_m$, respectively. Susceptible humans ($S_h$) become infected with a force of infection denoted by $\lambda_m^s$ for sensitive infections and $\lambda_m^r$ for resistant infections. Similarly, susceptible mosquitoes ($S_m$) become infected with a force of infection denoted by $\lambda_h^s$ for sensitive infections and $\lambda_h^r$ for resistant infections. A fraction $p$ of infected individuals will never develop clinical symptoms and will enter the group of paucisymptomatic infections ($P_h$). The rest will initially progress through an asymptomatic phase ($A_h$) before developing into symptomatic or clinical cases ($I_h$). Paucisymptomatic and asymptomatic infections are treated with antimalarial drugs without a curative intention at a rate $\phi$, leading to the partitioning of compartments $P_h$ and $A_h$ into untreated ($U$) and treated ($T$) sub compartments. Similarly, symptomatic infections are treated with a curative intention at a rate $\theta$, resulting in the partitioning of the symptomatic compartment into untreated ($U$) and treated ($T$) sub compartments. The implementation of MDA efforts is represented by $g$. $A_h$ infections progress to $I_h$ infections at a rate $v_h$, while $I_h$ infections recover at a rate $\gamma_h$. Recovered individuals lose their temporary immunity at a rate $\rho_h$. The death rates are either natural $\mu_k$ (for $k = h, m$) or due to infections for humans $\delta_h$. Infectious mosquitoes die at a rate $\delta_m^k$ (for $k = s, r$). The Notation is shown in Tables 1 and 2.

<table><caption>**Table 1.** Variable parameters, their biological meaning, and value.</caption><thead><tr><th>Parameters</th><th>Description (unit)</th><th>Value [Ref]</th></tr></thead><tbody><tr><td>φ</td><td>Rate of exposure of asymptomatic and paucisymptomatic infections to antimalarial treatment without a curative intention</td><td>(0; 1)</td></tr><tr><td>θ</td><td>Rate of exposure of clinical cases to antimalarial treatment with a curative intention</td><td>(0; 1)</td></tr><tr><td>q<sub>mda</sub></td><td>Proportion of humans under MDA</td><td>(0; 1)</td></tr><tr><td>ε<sub>s</sub></td><td>Probability of mutation from treated sensitive human infections to treated resistant human infections</td><td>10<sup>-1</sup> – 10<sup>-3</sup> (assumed)</td></tr><tr><td>γ&#773;<sub>h</sub></td><td>Ratio of the recovered rate</td><td>(0; 1)</td></tr><tr><td>δ&#773;<sub>h</sub></td><td>Ratio of the additional mortality rate of humans due to disease</td><td>(0; 1)</td></tr><tr><td>ε&#773;</td><td>Ratio of the mutation rate</td><td>(0; 1)</td></tr><tr><td>β&#773;<sub>m</sub></td><td>Ratio of the infection rate of human by infectious mosquitoes</td><td>(0; 1)</td></tr><tr><td>δ&#773;<sub>m</sub></td><td>Ratio of the mortality rate of infectious mosquitoes</td><td>(0; 1)</td></tr><tr><td>InitPrevR</td><td>Initial prevalence of drug resistance within the host population</td><td>{10<sup>-1</sup>, 10<sup>-3</sup>} (assumed)</td></tr><tr><td>ε<sub>r</sub></td><td>Probability of mutation from treated resistant human infections to treated sensitive human infections (day<sup>-1</sup>)</td><td>ε&#773;ε<sub>s</sub></td></tr><tr><td>δ<sup>r,T</sup><sub>h</sub></td><td>Additional death rate of treated infectious humans to resistant infection (day<sup>-1</sup>)</td><td>δ&#773;<sub>h</sub>δ<sup>s,U</sup><sub>h</sub></td></tr><tr><td>β<sup>r</sup><sub>m</sub></td><td>Resistant infection rate of human by infectious mosquitoes (day<sup>-1</sup> bites<sup>-1</sup>)</td><td>β&#773;<sub>m</sub>β<sup>s</sup><sub>m</sub></td></tr><tr><td>δ<sup>r</sup><sub>m</sub></td><td>Mortality rate of mosquitoes with resistant infection (day<sup>-1</sup>)</td><td>δ&#773;<sub>m</sub>δ<sup>s</sup><sub>m</sub></td></tr><tr><td>γ<sup>r,T</sup><sub>h</sub></td><td>Recovery rate of treated humans to resistant infection (day<sup>-1</sup>)</td><td>γ&#773;<sub>h</sub>γ<sup>s,T</sup><sub>h</sub></td></tr></tbody></table>
<table><caption>**Table 2.** Fixed parameters, their biological meaning and value.</caption><thead><tr><th>Parameters</th><th>Description (unit)</th><th>Value [Ref]</th></tr></thead><tbody><tr><td>T<sub>0,0</sub></td><td>Basic reproduction number for φ = 0, θ = 0, and g = 0</td><td>5 (assumed)</td></tr><tr><td>γ<sub>h</sub><sup>s,T</sup></td><td>Recovery rate of treated humans to sensitive infection (day<sup>-1</sup>)</td><td>0.0035 [15]</td></tr><tr><td>γ<sub>h</sub><sup>s,U</sup></td><td>Recovery rate of untreated humans to sensitive infection (day<sup>-1</sup>)</td><td>0 (assumed)</td></tr><tr><td>γ<sub>h</sub><sup>r,U</sup></td><td>Recovery rate of untreated humans to resistant infection (day<sup>-1</sup>)</td><td>0 (assumed)</td></tr><tr><td>β<sub>h</sub><sup>s</sup></td><td>Sensitive infection rate of mosquitoes by infected humans (day<sup>-1</sup> bites<sup>-1</sup>)</td><td>β<sub>h</sub> defined by (4.1)</td></tr><tr><td>β<sub>h</sub><sup>r</sup></td><td>Resistant infection rate of mosquitoes by infected humans (day<sup>-1</sup> bites<sup>-1</sup>)</td><td>β<sub>h</sub> defined by (4.1)</td></tr><tr><td>β<sub>m</sub><sup>s</sup></td><td>Sensitive infection rate of humans by infectious mosquitoes (day<sup>-1</sup> bites<sup>-1</sup>)</td><td>0.022 [15]</td></tr><tr><td>δ<sub>h</sub><sup>s,T</sup></td><td>Additional death rate of treated infectious humans from sensitive infection (day<sup>-1</sup>)</td><td>0 (assumed)</td></tr><tr><td>δ<sub>h</sub><sup>s,U</sup></td><td>Additional death rate of untreated infectious humans from sensitive infection (day<sup>-1</sup>)</td><td>1/5 [15]</td></tr><tr><td>δ<sub>h</sub><sup>r,U</sup></td><td>Additional death rate of untreated infectious humans from resistant infection (day<sup>-1</sup>)</td><td>9 × 10<sup>-5</sup> [15]</td></tr><tr><td>Λ<sub>h</sub></td><td>Recruitment rate of susceptible humans (day<sup>-1</sup> humans)</td><td>102 (assumed)</td></tr><tr><td>μ<sub>h</sub></td><td>Natural mortality rate of humans (day<sup>-1</sup>)</td><td>0.00224 (assumed)</td></tr><tr><td>μ<sub>m</sub></td><td>Natural mortality rate of mosquitoes (day<sup>-1</sup>)</td><td>1/20 (assumed)</td></tr><tr><td>Λ<sub>m</sub></td><td>Recruitment rate of susceptible mosquitoes (day<sup>-1</sup> mosquitoes)</td><td>100 (assumed)</td></tr><tr><td>δ<sub>m</sub><sup>s</sup></td><td>Mortality rate of mosquitoes with sensitive infection (day<sup>-1</sup>)</td><td>0.033 [15]</td></tr><tr><td>ρ<sub>h</sub></td><td>Loss of immunity rate of recovered humans (day<sup>-1</sup>)</td><td>5.5 × 10<sup>-4</sup> [15]</td></tr><tr><td>ν<sub>h</sub></td><td>The rate at which asymptomatic humans become infectious (day<sup>-1</sup>)</td><td>0.1 [15]</td></tr><tr><td>p</td><td>The proportion of paucisymptomatic infections (day<sup>-1</sup>)</td><td>0.7 (assumed)</td></tr><tr><td>D<sub>mda</sub></td><td>Delay to being no longer a carrier after MDA (days)</td><td>1 [16]</td></tr><tr><td>c<sub>h</sub></td><td>Balancing coefficients in the objective function ($ per person)</td><td>1 (assumed)</td></tr><tr><td>b<sub>h</sub></td><td>Balancing coefficients in the objective function ($ per person)</td><td>1 (assumed)</td></tr><tr><td>T<sub>f</sub></td><td>End time of controls (year)</td><td>10 (assumed)</td></tr><tr><td colspan="3" style="text-align: center;">Cost of control</td></tr><tr><td>c<sub>φ</sub></td><td>Cost of φ ($ per person)</td><td>10 [17–19]</td></tr><tr><td>c<sub>θ</sub></td><td>Cost of θ ($ per person)</td><td>3 – 20 [17–19]</td></tr><tr><td>c<sub>g</sub></td><td>Cost of MDA ($ per person)</td><td>3 – 20 [17–19]</td></tr></tbody></table>

### 3. Some basic properties of the model

The aim of this section is to provide some preliminary remarks to Systems (2.1) and (2.2). Let us introduce the following vector notations:

$$ 
\begin{aligned}
A_h &= (A_h^{s,T}, A_h^{r,T}, A_h^{s,U}, A_h^{r,U})', & P_h &= (P_h^{s,T}, P_h^{r,T}, P_h^{s,U}, P_h^{r,U})', & I_h &= (I_h^{s,T}, I_h^{r,T}, I_h^{s,U}, I_h^{r,U})', \\
R_h &= (R_h^{s,T}, R_h^{r,T}, R_h^{s,U}, R_h^{r,U})', & I_m &= (I_m^s, I_m^r),
\end{aligned}
$$
as well as the matrices $\beta_h = \text{diag}(\beta_h^k)_{k \in \{s,r\}}$, $\beta_m = \text{diag}(\beta_m^k)_{k \in \{s,r\}}$, $\gamma_h = \text{diag}(\gamma_h^{k,l})_{k \in \{s,r\}, l \in \{T,U\}}$, $\delta_h = \text{diag}(\delta_h^{k,l})_{k \in \{s,r\}}$, $\delta_m = \text{diag}(\delta_m^k)_{k \in \{s,r\}}$,

$$ \Theta = \begin{pmatrix} \theta & 0 \\ 0 & \theta \\ 1 - \theta & 0 \\ 0 & 1 - \theta \end{pmatrix}, \quad \Phi = \begin{pmatrix} \phi & 0 \\ 0 & \phi \\ 1 - \phi & 0 \\ 0 & 1 - \phi \end{pmatrix}, \quad e = \begin{pmatrix} 1 & 0 & 1 & 0 \\ 0 & 1 & 0 & 1 \end{pmatrix}, \quad \text{and} \quad \varepsilon = \begin{pmatrix} 1 - \varepsilon_s & \varepsilon_r & 1 & 0 \\ \varepsilon_s & 1 - \varepsilon_r & 0 & 1 \end{pmatrix}. $$

Then Models (2.1) and (2.2) are rewritten as

$$ \left\{ \begin{aligned} \dot{S}_h &= \Lambda_h - \mu_h S_h - \mathbf{1}_2 \lambda_m S_h + \rho_h \mathbf{1}_4 R_h, \\ \dot{A}_h &= (1-p)\Phi(t)\lambda_m S_h - (\nu_h + \mu_h + g(t))A_h, \\ \dot{P}_h &= -\mu_h P_h + p\Phi(t)\lambda_m S_h - g(t)P_h, \\ \dot{I}_h &= -\mu_h I_h + \Theta(t)ev_h A_h - \gamma_h I_h - \delta_h I_h, \\ \dot{R}_h &= \gamma_h I_h + g(t)A_h + g(t)P_h - (\mu_h + \rho_h)R_h, \\ \dot{S}_m &= \Lambda_m - \mu_m S_m - \mathbf{1}_2 \lambda_h S_m, \\ \dot{I}_m &= \lambda_h S_m - \delta_m I_m, \end{aligned} \right. \quad (3.1) $$

with $\mathbf{1}_2 = (1, 1)$, $\mathbf{1}_4 = (1, 1, 1, 1)$, and

$$ \lambda_h = \lambda_h(S_h, A_h, P_h, I_h, R_h) = \frac{\beta_h}{N_h} \varepsilon(A_h + P_h + I_h), $$

$$ \lambda_m = \lambda_m(S_h, A_h, P_h, I_h, R_h, I_m) = \frac{1}{N_h} \beta_m I_m. $$

Let us introduce the matrix $\mathcal{B} = -\text{diag}(\mu_h, (\mu_h + g(t) + \nu_h)\mathbf{I}_d, (\mu_h + g(t))\mathbf{I}_d, (\mu_h + \gamma_h + \delta_h)\mathbf{I}_d, (\mu_h + \rho_h)\mathbf{I}_d, \mu_m, \delta_m)$, where $\mathbf{I}_d$ is the identity matrix, as well as the nonlinear map

$$ \mathcal{M}[t](S_h, A_h, P_h, I_h, R_h, S_m, I_m) = \begin{pmatrix} \Lambda_h - \mathbf{1}_2 \lambda_m S_h + \rho_h \mathbf{1}_4 R_h \\ (1-p)\Phi(t)\lambda_m S_h \\ p\Phi(t)\lambda_m S_h \\ \Theta(t)ev_h A_h \\ g(t)(A_h + P_h) + \gamma_h I_h \\ \Lambda_m - \mathbf{1}_2 \lambda_h S_m \\ \lambda_h S_m \end{pmatrix}. $$

By identifying $v$ together with $(S_h, A_h, P_h, I_h, R_h, S_m, I_m)$ and by setting $v_0 = (S_h(0), A_h(0), P_h(0), I_h(0), R_h(0), S_m(0), I_m(0))$, Systems (2.1) and (2.2) are rewritten as the following Cauchy problem:

$$ \left\{ \begin{aligned} \frac{dv(t)}{dt} &= \mathcal{B}v(t) + \mathcal{M}[t](v(t)), \\ v(0) &= v_0. \end{aligned} \right. \quad (3.2) $$

### 3.1. Existence and uniqueness of the global solution

The existence and uniqueness of the bounded solutions of (3.2) reads as follows:
**Theorem 3.1.** Let $v_0 \in \mathbb{R}_+^{20}$. Assume that $\theta, \phi$, and $g$ are constant functions. Then the following hold.

(1) A unique global solution $\mathcal{V}(\cdot, v_0) : [0, \infty) \to \mathbb{R}^{20}$ exists such that $\mathcal{V}(t, v_0) = v(t)$ and $\mathcal{V}(0, v_0) = v_0$. Furthermore,

$$ \mathcal{V}(t, v_0) = e^{\mathcal{B}t} v_0 + \int_0^t e^{\mathcal{B}(t-\eta)} \mathcal{M}(v(\eta)) d\eta. $$

(2) The semiflow defined by $\{\mathcal{V}(t, v_0)\}_t$ is bounded, dissipative, and asymptotically smooth, and hence it admits a global attractor in $\mathbb{R}_+^{20}$.

We recall the following definition.

**Bounded dissipative:** A bounded set $\Omega \subset \mathbb{R}^{20}$ exists such that for any bounded set $U \subset \mathbb{R}^{20}$, $\tau = \tau(U, \Omega) \ge 0$ exists such that $\mathcal{V}(t, U) \subset \Omega$ for $t \ge \tau$.

**Asymptotically smooth:** For any nonempty, closed, bounded set $U \subset \mathbb{R}^{20}$, a nonempty compact set $J = J(U)$ exists such that $J$ attracts $\{\varphi \in U : \mathcal{V}(t, \varphi) \in U, \forall t \ge 0\}$.

Furthermore, we also recall that a nonempty set $J \subset \mathbb{R}^{20}$ is said to attract a nonempty set $U \subset \mathbb{R}^{20}$ if $\delta(\mathcal{V}(t, U), J) \to 0$ as $t \to +\infty$, where $\delta(U, J) = \sup_{u \in U} \inf_{v \in J} \|u - v\|$.

Note that the term “asymptotically smooth” is used here following the terminology of Hale, LaSalle, and Slemrod [20]. However, this is strictly equivalent to the notion of *asymptotically compact* as defined by Ladyzhenskaya [21]. Finally, refer to Section 7.1 for the detailed proof of Theorem 3.1.

## 3.2. The disease-free equilibrium and disease invasion process

To determine the threshold of System (3.2), which govern the invasion dynamics of the disease in a host population that is completely free of infections, we linearize (3.2) around the disease-free equilibrium (DFE) $E_0$. Specifically, the DFE is written as

$$ E_0 = (S_h^0, 0_{\mathbb{R}^4}, 0_{\mathbb{R}^4}, 0_{\mathbb{R}^4}, S_m^0, 0_{\mathbb{R}^2}), \quad S_h^0 = \frac{\Lambda_h}{\mu_h}, \quad S_m^0 = \frac{\Lambda_m}{\mu_m}, $$

and we linearize the infective compartments of (3.2) at the disease-free equilibrium $E_0$ to obtain the following system:

$$ \frac{dw(t)}{dt} = \hat{B}w(t) + \hat{M}w(t), $$

where $w(t) = (A_h(t), P_h(t), I_h(t), I_m(t))$, and

$$ \hat{B} = \begin{pmatrix} -(\mu_h + v_h + g)\mathbf{I}_d & 0_{4 \times 4} & 0_{4 \times 4} & 0_{4 \times 2} \\ 0_{4 \times 4} & -(\mu_h + g)\mathbf{I}_d & 0_{4 \times 4} & 0_{4 \times 2} \\ \mathbf{\Theta}e\nu_h & 0_{4 \times 4} & -(\mu_h \mathbf{I}_d + \gamma_h + \delta_h) & 0_{4 \times 2} \\ 0_{2 \times 4} & 0_{2 \times 4} & 0_{2 \times 4} & -\delta_m \end{pmatrix}, $$

$$ \hat{M} = \begin{pmatrix} 0_{4 \times 4} & 0_{4 \times 4} & 0_{4 \times 4} & (1-p)\Phi\beta_m \\ 0_{4 \times 4} & 0_{4 \times 4} & 0_{4 \times 4} & p\Phi\beta_m \\ 0_{4 \times 4} & 0_{4 \times 4} & 0_{4 \times 4} & 0_{4 \times 2} \\ \frac{\beta_h \varepsilon S_m^0}{S_h^0} & \frac{\beta_h \varepsilon S_m^0}{S_h^0} & \frac{\beta_h \varepsilon S_m^0}{S_h^0} & 0_{2 \times 2} \end{pmatrix}. $$

**Lemma 3.2.** Assume that $\theta, \phi$, and $g$ are constant functions. The following then hold.
(1) The square matrix $\hat{M}$ is positive, i.e., $\hat{M}$ is not the zero matrix, and all entries of $\hat{M}$ are non-negative.

(2) The spectral bound $\mathfrak{s}(\hat{B})$, i.e., $\mathfrak{s}(\hat{B}) = \sup \{\mathfrak{R}(\lambda) : \lambda \in \sigma(\hat{B})\}$, satisfies $\mathfrak{s}(\hat{B}) < 0$.

The proof of Lemma 3.2 is straightforward and thus the details are omitted.
We now define

$$ \mathcal{T}_0 := r(\hat{M}(-\hat{B})^{-1}), $$

where $r(\hat{M}(-\hat{B})^{-1})$ represents the spectral radius of the square matrix $\hat{M}(-\hat{B})^{-1}$.

The subsequent lemma provides a connection between the sign of $\mathcal{T}_0 - 1$ and the growth bound of the $C_0$-semigroup generated by $(\hat{B} + \hat{M})$.

**Lemma 3.3.** Assume that $\theta$, $\phi$, and $g$ are constant functions. In this case

$$ \operatorname{sgn}(\mathfrak{s}(\hat{B} + \hat{M})) = \operatorname{sgn}(\mathcal{T}_0 - 1). $$

Furthermore, we have

$$ \mathcal{T}_0 = \sqrt{\frac{S_m^0}{S_h^0} r(\mathcal{T}_h \mathcal{T}_m)}. $$

Here, the 2×2 matrices $\mathcal{T}_h$ and $\mathcal{T}_m$ describe the infectiousness of humans and mosquitoes, respectively, and are defined by

$$ \mathcal{T}_h = \beta_h \varepsilon \left\{ \left( (v_h + \mu_h + g)^{-1} I_d (1-p) + (g + \mu_h)^{-1} I_d p + (\mu_h I_d + \gamma_h + \delta_h)^{-1} \Theta e v_h \right) (v_h + \mu_h + g)^{-1} I_d (1-p) \right\} \Phi, \quad (3.3) $$

$$ \mathcal{T}_m = \beta_m \delta_m^{-1} = \operatorname{diag}(\beta_m^s \delta_m^{s-1}, \beta_m^r \delta_m^{r-1}). $$

Through a straightforward computation, an explicit formula of $\mathcal{T}_0$ is given by

$$ \mathcal{T}_0^2 = \frac{(\mathcal{T}_0^s + \mathcal{T}_0^r) + \sqrt{(\mathcal{T}_0^s - \mathcal{T}_0^r)^2 + 4Q_{sr}}}{2}, \quad (3.4) $$

where

$$ \mathcal{T}_0^s = \frac{S_m^0 \beta_m^s}{S_h^0 \delta_m^s} \beta_h^s \mathcal{K}_h^s, \quad \mathcal{T}_0^r = \frac{S_m^0 \beta_m^r}{S_h^0 \delta_m^r} \beta_h^r \mathcal{K}_h^r, $$

$$ \mathcal{K}_h^s = (1 - \varepsilon_s) \underbrace{\left( (1 - p)\theta \frac{v_h}{(v_h + \mu_h + g)(\mu_h + \gamma_h^{s,T} + \delta_h^{s,T})} + (1 - p)\phi \frac{1}{(v_h + \mu_h + g)} + p\phi \frac{1}{(g + \mu_h)} \right)}_{\text{lifespan of s-infections exposed to antimalarial drugs}} + (1 - p)(1 - \theta) \underbrace{\frac{v_h}{(v_h + \mu_h + g)(\mu_h + \gamma_h^{s,U} + \delta_h^{s,U})} + (1 - p)(1 - \phi) \frac{1}{(v_h + \mu_h + g)} + p(1 - \phi) \frac{1}{(g + \mu_h)}}_{\text{lifespan of s-infections unexposed to antimalarial drugs}} $$

$$ \mathcal{K}_h^r = (1 - \varepsilon_r) \underbrace{\left( (1 - p)\theta \frac{v_h}{(v_h + \mu_h + g)(\mu_h + \gamma_h^{r,T} + \delta_h^{r,T})} + (1 - p)\phi \frac{1}{(v_h + \mu_h + g)} + p\phi \frac{1}{(g + \mu_h)} \right)}_{\text{lifespan of r-infections exposed to antimalarial drugs}} $$
$$ + (1-p)(1-\theta) \frac{v_h}{(v_h + \mu_h + g)(\mu_h + \gamma_h^{r,U} + \delta_h^{r,U})} + (1-p)(1-\phi) \frac{1}{(v_h + \mu_h + g)} + p(1-\phi) \frac{1}{(g + \mu_h)}, $$

lifespan of r-infections unexposed to antimalarial drugs

and

$$ Q_{sr} = \varepsilon_r \beta_h^r \frac{\beta_m^r}{\delta_m^r} \left( \frac{S_m^0}{S_h^0} \right) \left( \theta(1-p) \frac{v_h}{(\mu_h + \gamma_h^{r,T} + \delta_h^{r,T})(v_h + \mu_h + g)} + (1-p)\phi \frac{1}{(v_h + \mu_h + g)} + p\phi \frac{1}{g + \mu_h} \right) $$

average rate of new s-infection from resistant infections due to mutation

$$ \times \varepsilon_s \beta_h^s \frac{\beta_m^s}{\delta_m^s} \left( \frac{S_m^0}{S_h^0} \right) \left( \theta(1-p) \frac{v_h}{(\mu_h + \gamma_h^{s,T} + \delta_h^{s,T})(v_h + \mu_h + g)} + (1-p)\phi \frac{1}{(v_h + \mu_h + g)} + p\phi \frac{1}{g + \mu_h} \right). $$

average rate of new r-infection from sensitive infections due to mutation

Furthermore, employing the next-generation matrix approach introduced in Lemma 6.12 of [22], we can determine that $\mathcal{T}_0$ represents the basic reproduction number of the entire system, while $\mathcal{T}_0^s$ and $\mathcal{T}_0^r$ denote the basic reproduction numbers induced only by s- and r-infections, respectively. Consequently, by Lemma 3.3 and the results on the asymptotic stability of equilibria for the vector fields given in Theorem 3.16 of [22], we find that the DFE $E_0$ is locally asymptotically stable if and only if $\mathcal{T}_0 < 1$.

The proof of Lemma 3.3 is appended in Section 7.2.

**Remark 1.** Parameters $\mathcal{K}_h^s$ and $\mathcal{K}_h^r$ can be interpreted as lifespans of s- and r-infections, respectively. For instance, $\mathcal{K}_h^s$ is rewritten as

$$ \mathcal{K}_h^s = (1 - \varepsilon_s)(1 - p)\theta \frac{v_h}{(v_h + \mu_h + g)} D_{s,T} + (1 - \varepsilon_s)(1 - p)\phi D_A + (1 - \varepsilon_s)p\phi D_P \\ + (1 - p)(1 - \theta) \frac{v_h}{(v_h + \mu_h + g)} D_{s,U} + (1 - p)(1 - \phi) D_A + p(1 - \phi) D_P, $$

where $D_P = (\mu_h + g)^{-1}$, $D_A = (\mu_h + v_h + g)^{-1}$, $D_{s,l} = (\mu_h + \gamma_h^{s,l} + \delta_h^{s,l})^{-1}$, and $l = T, U$ are the lifespans of paucisymptomatic, asymptomatic, and clinical infections, respectively. From this, since $(1 - \varepsilon_s)(1 - p)\theta \frac{v_h}{(v_h + \mu_h + g)}$, $(1 - \varepsilon_s)(1 - p)\phi$, and $(1 - \varepsilon_s)p\phi$ represent the proportions of s-clinical, s-asymptomatic, and s-paucisymptomatic infections exposed to antimalarial drugs, respectively, we can easily interpret the term

$$ (1 - \varepsilon_s)(1 - p)\theta \frac{v_h}{(v_h + \mu_h + g)} D_{s,T} + (1 - \varepsilon_s)(1 - p)\phi D_A + (1 - \varepsilon_s)p\phi D_P, $$

as the lifespan of s-infections exposed to antimalarial drugs. By the same way, the term

$$ (1 - p)(1 - \theta) \frac{v_h}{(v_h + \mu_h + g)} D_{s,U} + (1 - p)(1 - \phi) D_A + p(1 - \phi) D_P, $$

represents the lifespan of s-infections unexposed to antimalarial drugs. Consequently, $\mathcal{K}_h^s$ is then the lifespan of the overall s-infections.

A similar interpretation applies to $\mathcal{K}_h^r$. As for $Q_{sr}$, it can be understood as the product of the average rate of new s-infections arising from resistant infections and the average rate of new r-infections emerging from sensitive infections due to mutations.
### 3.3. The endemic equilibrium and existence of a backward bifurcation

Assume that $\theta, \phi$, and $g$ are constant functions. Let

$$E_h^* = (S_h^*, A_h^*, P_h^*, I_h^*, R_h^*, S_m^*, I_m^*)' \in \mathbb{R}_+^{20},$$

be an endemic equilibrium of System (3.1), and

$$\lambda_h^{s*} = \frac{\beta_h^s}{N_h^*} [A_h^{s,U^*} + P_h^{s,U^*} + I_h^{s,U^*} + \varepsilon_r (A_h^{r,T^*} + P_h^{r,T^*} + I_h^{r,T^*}) + (1 - \varepsilon_s) (A_h^{s,T^*} + P_h^{s,T^*} + I_h^{s,T^*})],$$

$$\lambda_h^{r*} = \frac{\beta_h^r}{N_h^*} [A_h^{r,U^*} + P_h^{r,U^*} + I_h^{r,U^*} + (1 - \varepsilon_r) (A_h^{r,T^*} + P_h^{r,T^*} + I_h^{r,T^*}) + \varepsilon_s (A_h^{s,T^*} + P_h^{s,T^*} + I_h^{s,T^*})].$$

The next theorem provides the necessary conditions on which we have the existence of a unique endemic equilibrium and the demonstration of a backward bifurcation. We examine the occurrence of a bifurcation of the endemic equilibrium in Model (3.1) at $\mathcal{T}_0 = 1$. We find that an endemic equilibrium of (3.1) exists if and only if a positive constant $\kappa$ exists such that

$$L(\mathcal{T}_0^2, \kappa) = 1,$$

where $L$ is defined by

$$L(\mathcal{T}_0^2, \kappa) = \Lambda_m \frac{S_h^0}{S_m^0} \mathcal{T}_0^2 \frac{(\mu_h + \chi_1 \kappa)}{(1 + \chi_2 \kappa) (\Lambda_h \mu_m + (\chi_2 + \Lambda_h \mu_m \omega^s + \chi_3) \kappa)},$$

where $\chi_1 = \omega^r (1 - a_s) + \omega^s (1 - a_r)$, $\chi_2 = \omega^r (b_s + a_s) + \omega^s (b_r + a_r)$, and $\chi_3 = \omega^s (c_s + d_s) + \omega^r (c_r + d_r)$, with $\omega = (\omega^s, \omega^r)$ being the principal eigenvector of the matrix $\mathcal{T}_m \mathcal{T}_h$, and with $\|\omega\| = 1$; $a_j, b_j, c_j$, and $d_j$'s are positive constants given later.

Next, we introduce the following bifurcation constant:

$$K = \Lambda_h \mu_m \chi_1 - (2\chi_2 + \Lambda_h \mu_m \omega^s + \chi_3) \mu_h.$$

We then obtain the following result on the existence and bifurcation of the endemic equilibrium.

**Theorem 3.4.** Assume that $\theta, \phi$, and $g$ are constant functions. We then have the following.

(1) If $K > 0$, then a backward bifurcation occurs at $\mathcal{T}_0 = 1$, meaning that for $\mathcal{T}_0$ values slightly less than 1, two endemic equilibria exist for Model (3.1).

(2) If $K < 0$, then a forward bifurcation occurs at $\mathcal{T}_0 = 1$, meaning that for $\mathcal{T}_0$ values slightly greater than 1, a unique endemic equilibrium exists, while for $\mathcal{T}_0$ values slightly less than 1, no endemic equilibrium exists.

(3) If $\mathcal{T}_0 > 1$, then at least one endemic equilibrium exists for (3.1). Moreover, if $K > 0$, then when $\mathcal{T}_0 = 1$, an endemic equilibrium also exists.

(4) Assume that the map $\kappa \mapsto L(\mathcal{T}_0^2, \kappa)$ is strictly decreasing on $\mathbb{R}_+$. Thus, an endemic equilibrium exists if and only if $\mathcal{T}_0 > 1$; in this case, the equilibrium is unique.

We refer to Section 7.3 for the proof of Theorem 3.4.
# 4. Model parameters, baseline model simulation, and sensitivity analysis

For the illustrative dynamics in this section, we assume that $\theta$, $\phi$, and $g$ are constant functions.

## 4.1. Model parameters and initial conditions

The parameter values used in the model are summarized in Tables 1 and 2. In general, we assume that the parameters for sensitive infections are well-established for both humans and vectors. To estimate parameters for resistant infections, we introduce variable scaling parameters. For example, the mutation probability $\varepsilon_s$ for transitions from sensitive to resistant infections is a known value. The mutation probability $\varepsilon_r$ for transitions from resistant to sensitive infections is then estimated as $\varepsilon_r = \bar{\varepsilon}\varepsilon_s$, where $\bar{\varepsilon} \in (0, 1)$ is a scaling parameter that represents the ratio of the mutation rates. This approach ensures a consistent framework for parameter estimation across sensitive and resistant infections.

The transmission rate $\beta_h$, which represents the rate of transmission from infected humans to mosquitoes and is assumed to be identical for both sensitive and resistant infections, is determined using the basic reproduction number $T_0$ as given by (3.4). Specifically, $\beta_h$ is calibrated on the basis of the value of $T_{0,0}$ under the baseline scenario where $\phi = 0$, $\theta = 0$, and $g = 0$. This ensures that the transmission dynamics align with the expected epidemiological behavior under the given conditions. More precisely, we have

$$ \beta_h = \frac{T_{0,0}^2 S_h^0}{S_m^0 \max\{C_0^s, C_0^r\}}, \quad (4.1) $$

with

$$ C_0^s = \frac{\beta_m^s}{\delta_m^s} \left( (1-p) \frac{\nu_h}{(\nu_h + \mu_h)(\mu_h + \gamma_h^{s,U} + \delta_h^{s,U})} + (1-p) \frac{1}{(\nu_h + \mu_h)} + p \frac{1}{\mu_h} \right), $$

$$ C_0^r = \frac{\beta_m^r}{\delta_m^r} \left( (1-p) \frac{\nu_h}{(\nu_h + \mu_h)(\mu_h + \gamma_h^{r,U} + \delta_h^{r,U})} + (1-p) \frac{1}{(\nu_h + \mu_h)} + p \frac{1}{\mu_h} \right). $$

The initial conditions are specified using parameters such as the initial prevalence of drug resistance (InitPrevR) within the host population. We assume an initial infection prevalence of Prev = 50% with PrevI = 10% of these infections being clinical cases. In addition, we suppose that a proportion ($p = 0.7$) of infected humans will develop paucisymptomatic infections. The initial susceptible human population, asymptomatic infected, paucisymptomatic infected, and infectious case are, respectively, given by $S_h(0) = (1 - \text{Prev})S_h^0$; $A_h(0) = (1 - p)S_h^0(1 - \text{PrevI})\text{Prev}(1 - \text{InitPrevR}, \text{InitPrevR}, 0, 0)$; $P_h(0) = pS_h^0\text{Prev}(1 - \text{InitPrevR}, \text{InitPrevR}, 0, 0)$; $I_h(0) = (1 - p)S_h^0\text{Prev}\text{PrevI}(1 - \text{InitPrevR}, \text{InitPrevR}, 0, 0)$; and $R_h(0) = (0, 0, 0, 0)$. For mosquitoes, we assume that $S_m(0) = 0.9S_m^0$ and $I_m(0) = 0.1S_m^0(1 - \text{InitPrevR}, \text{InitPrevR})$.

## 4.2. Typical dynamics simulated by the model

Model (3.1) is employed to characterize the malaria infection burden within the human host population and to capture the dynamics of resistant clinical cases.

When the probability of mutation from treated sensitive human infections to treated resistant human infections ($\varepsilon_s$) is relatively low (Figure 2), the exposure rate $\phi$ of asymptomatic and paucisymptomatic infections to antimalarial treatment has a marginal impact on both mortality and the number of resistant clinical cases (Figure 2, Column 1 vs. Column 2 vs. Column 3). However, while increasing the antimalarial treatment rate $\theta$ for clinical cases has a relatively marginal effect on the cumulative number of
deaths due to infection when $q_{mda}$ is moderate (Figure 2A–C), such an increase significantly reduces infection-related mortality when $q_{mda}$ is relatively high (Figure 2D–I). Likewise, increasing the proportion of humans covered under MDA has a pronounced effect in reducing both mortality and resistant clinical cases, with a particularly strong impact on the latter (Figure 2). Importantly, increasing the treatment rate for clinical cases is associated with a rise in the number of resistant clinical cases (Figure 2A–F), except in the scenario where the proportion $q_{mda}$ is sufficiently strong (Figure 2G–I). These contrasting behaviors highlight the beneficial effect of combining antimalarial treatment and MDA to simultaneously reduce infection-related mortality and the number of resistant clinical cases.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Dynamics of the model when the mutation rate $\varepsilon_s = 0.001$. Parameter values are the defaults in Tables 1 and 2.</figcaption></figure>

The configuration above remains similar under high mutation rates ($\varepsilon_s$) (Figure 3). However, a high
mutation rate $\varepsilon_s$ leads to a higher number of deaths within the host population (Figure 3) compared with the case of a moderate mutation rate (Figure 2). Furthermore, high values of $\varepsilon_s$ accelerate the emergence of antimalarial resistance within the host population (Figure 3A–F) compared with scenarios with low mutation rates (Figure 2A–F). Although extensive use of MDA helps reduce both deaths and resistant cases, the reduction in resistant cases occurs more slowly under high mutation rates (Figure 3G–I) than in the case of low mutation rates (Figure 2G–I).

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Dynamics of the model when the mutation rate $\varepsilon_s = 0.1$. Parameter values are the defaults in Tables 1 and 2.</figcaption></figure>
Finally, for Figures 2D–I and 3D–I, one may observe a peak in resistant clinical cases at the initial time $t = 0$. This peak is due to our chosen initial data. However, in Figure 3G–I, we observe an increase in resistant clinical cases during the first year before they start to decline. This initial increase can be explained by the high mutation rate $\varepsilon_s$. Indeed, although the proportion of humans covered under MDA ($q_{mda}$) is sufficiently high, the strong mutation capability initially leads to an increase in the number of resistant clinical cases before they begin to decline (Figure 3G–I). In contrast, this pattern is not observed when the mutation rate is low (Figure 2G–I).

### 4.3. Global sensitivity analysis

We analyze the sensitivity of cumulative deaths and resistant cases (including paucisymptomatic, asymptomatic, and clinical cases) to five key parameters: the initial prevalence of drug resistance within the host population (InitPrevR), the probability of mutation from treated sensitive human infections to treated resistant human infections ($\varepsilon_s$), the rate of exposure of asymptomatic and paucisymptomatic infections to antimalarial treatment without a curative intention ($\phi$), the rate of exposure of clinical cases to antimalarial treatment with a curative intention ($\theta$), and the proportion of humans covered under MDA ($q_{mda}$). The variation range for these parameters is specified in Table 1. Sensitivity indices are calculated by fitting an analysis of variance (ANOVA) linear model, including up to third-order interactions, to simulation-generated data. This ANOVA model demonstrates a high level of accuracy, explaining more than 99% of the variance.

Overall, the sensitivity analysis identifies $q_{mda}$ as the primary parameter influencing both the number of deaths and the prevalence of resistant cases (Figure A1). It accounts for 89% of the variance in the number of deaths (Figure A1, left) and 53% in the prevalence of resistant cases (Figure A1, right). When interactions with other parameters are included, these indices increase to 94% and 83%, respectively. The second most influential parameter is $\theta$, which explains 5% of the variance in the number of deaths (rising to 10% when accounting for interactions with other parameters) and 14% of the variance in the prevalence of resistance (rising to 44% with interactions). The parameters InitPrevR and $\varepsilon_s$ have a comparable and modest impact, each contributing 0.2%–0.4% to the variance in both outcomes. Finally, $\phi$ has a negligible effect on both the number of deaths and the prevalence of resistance, explaining less than 0.09% of the variance, even when accounting for interactions.

## 5. Optimal control

Although previous analyses were performed under the assumption of constant parameter values for the rate of exposure of asymptomatic and paucisymptomatic infections to antimalarial treatment without curative intent ($\phi$), the rate of exposure of clinical cases to antimalarial treatment with curative intent ($\theta$), and the proportion of humans covered under MDA ($g$). This section considers these parameters as time-dependent functions and assigns reasonable upper and lower bounds to account for potential temporal variations.

The control scheme ($\theta, \phi, g$) lies on the set

$$ \Omega = \{(\theta, \phi, g) \in (L^{\infty}(0, T_f))^{3} : 0 \le \theta \le 1; \ 0 \le \phi \le 1; \ 0 \le g \le g_{max}\}, $$

with $g_{max}$ being a positive constant.
Although our objective is not to provide a detailed analysis of the controlled system, i.e., System (3.1) where the control scheme $(\theta, \phi, g)$ is considered as a time-dependent function, we present the following well-posedness result.

**Theorem 5.1.** Assume that the time-dependent functions $\theta, \phi$, and $g$ are continuous. Then the controlled system, that is, the System (3.1), where the control scheme $(\theta(t), \phi(t), g(t))$ is considered as a time-dependent function, generating a unique, globally defined, positive, and bounded non autonomous semiflow.

The proof of this result relies on a standard methodology (e.g., [23, 24]).

## 5.1. The optimal control problem

We assume that a successful control scheme is one that effectively reduces both the total number of deaths and the total number of resistant cases. Consequently, the control scheme is considered optimal if it minimizes the following objective functional:

$$ \mathcal{J}(\theta, \phi, g) = \int_{0}^{T_f} \left( c_h D_h(t) + b_h \sum_{\ell \in \{T, U\}} \left( A_h^{r,\ell}(t) + P_h^{r,\ell}(t) + I_h^{r,\ell}(t) \right) + \frac{c_\theta}{2}\theta^2(t) + \frac{c_\phi}{2}\phi^2(t) + \frac{c_g}{2}g^2(t) \right) dt. \quad (5.1) $$

Here, $b_h$ and $c_h$ are balancing coefficients that transform the integral into the cost expended over a finite period of $T_f$ years (see Table 2). The first integral, weighted by $b_h$ and $c_h$, captures the economic losses associated with deaths and resistant infections, respectively. The second integral represents the costs of implementing the three control measures, with quadratic terms for the controls included to account for the nonlinear costs.

We find that $(\theta^*, \phi^*, g^*)$ satisfies

$$ \mathcal{J}(\theta^*, \phi^*, g^*) = \min_{(\theta, \phi, g) \in \Omega} \mathcal{J}(\theta, \phi, g). \quad (5.2) $$

## 5.2. Necessary optimality conditions

To address the necessary optimality conditions, we follow the framework provided by [25–27]. Given a solution $v = (S_h, A_h, P_h, I_h, R_h, S_m, I_m)$ of System (3.1), where the control scheme $(\theta(t), \phi(t), g(t))$ is considered as a time-dependent function, we define the corresponding adjoint functions as $\vartheta = (\vartheta_{S_h}, \vartheta_{A_h}, \vartheta_{P_h}, \vartheta_{I_h}, \vartheta_{R_h}, \vartheta_{S_m}, \vartheta_{I_m}, \vartheta_{D_h})$. Theorem 5.1 gives the existence and boundedness of the solutions of System (3.1) with time-dependent control functions. According to the results in [25, 27], the set of controls $\Omega$ is nonempty, convex, and closed. Moreover, the right-hand side of System (3.1) is continuous and can be expressed as a linear function of the controls $\theta, \phi$, and $g$, with coefficients that depend on both time and the state variables.

If we choose $a_1 \in \mathbb{R}_+^*$, $a_2 = \min\left\{\frac{c_\theta}{2}, \frac{c_\phi}{2}, \frac{c_g}{2}\right\}$, and $\alpha = 2$, the integrant $h(v, \theta, \phi, g)(t)$ of $\mathcal{J}$ verifies

$$ \begin{aligned} h(v, \theta, \phi, g)(t) = &c_h D_h(t) + b_h \sum_{\ell \in \{T, U\}} \left( A_h^{r,\ell}(t) + P_h^{r,\ell}(t) + I_h^{r,\ell}(t) \right) + \frac{c_\theta}{2}\theta^2 + \frac{c_\phi}{2}\phi^2 + \frac{c_g}{2}g^2 \\ &\geq \frac{c_\theta}{2}\theta^2 + \frac{c_\phi}{2}\phi^2 + \frac{c_g}{2}g^2 \end{aligned} $$
$$ \geq -a_1 + \min \left\{ \frac{c_{\theta}}{2}, \frac{c_{\phi}}{2}, \frac{c_{g}}{2} \right\} (\theta^2 + \phi^2 + g^2). $$

Therefore, the existence of the optimal control is given by the following result.

**Theorem 5.2.** Assume that the time-dependent functions $\theta$, $\phi$, and $g$ are continuous. Let $v = (S_h, A_h, P_h, I_h, R_h, S_m, I_m)$ be the solution of System (3.1) where the control scheme $(\theta, \phi, g)$ is considered as a time-dependent function. Here, $(\theta^*, \phi^*, g^*) \in \Omega$ exists such that

$$ \min_{(\theta, \phi, g) \in \Omega} \mathcal{J}(\theta, \phi, g) = \mathcal{J}(\theta^*, \phi^*, g^*). $$

The existence of such a control is well-established using standard techniques, as detailed in [25, 27].

Pontryagin's maximum principle [14] allows us to reformulate the problem (5.2) to the problem of minimizing the Hamiltonian defined as follows:

$$ \mathcal{H} = c_h D_h + b_h \sum_{\ell \in \{T, U\}} (A_h^{r, \ell} + P_h^{r, \ell} + I_h^{r, \ell}) + \frac{c_{\theta}}{2}\theta^2 + \frac{c_{\phi}}{2}\phi^2 + \frac{c_{g}}{2}g^2 + \sum_{z \in \mathcal{Z}} \vartheta_z f_z, $$

where $\mathcal{Z} = \{S_h, A_h, P_h, I_h, R_h, D_h, S_m, I_m\}$, $(\vartheta_z)_{z \in \mathcal{Z}}$ are adjoint functions and $f_z$ is given by the right-hand side of Models (2.1) and (2.2) for the z-compartment.

The necessary conditions for the existence of the solution to the problem (5.2) are

$$ \dot{\vartheta}_z = -\nabla_z \mathcal{H}, \quad \text{for } z \in \mathcal{Z}, \qquad (5.3) $$

and

$$ \frac{\partial \mathcal{H}}{\partial w} = 0, \quad \text{for } w \in \{\theta, \phi, g\}. \qquad (5.4) $$

From (5.3), we find that the adjoint system $(\vartheta_k)_{k \in \mathcal{Z}}$ verifies:

$$ \left\{ \begin{aligned} \dot{\vartheta}_{S_h} &= \left(\mu_h + 1_2 \lambda_m \left(1 - \frac{S_h}{N_h}\right)\right) \vartheta_{S_h} - \left(1 - \frac{S_h}{N_h}\right) (1-p)(\Phi \lambda_m)^t \vartheta_{A_h} \\ &\quad - \left(1 - \frac{S_h}{N_h}\right) p (\Phi \lambda_m)^t \vartheta_{P_h} - 1_2 \lambda_h \frac{S_m}{N_h} \vartheta_{S_m} + \frac{S_m}{N_h} \langle \lambda_h, \vartheta_{I_m} \rangle, \\ \dot{\vartheta}_{A_h} &= b - \frac{S_h}{N_h} (1_2 \lambda_m) 1_4^t \vartheta_{S_h} + (\mu_h + v_h + g(t)) \vartheta_{A_h} + \frac{S_h}{N_h} (1-p) (1_4 (\Phi \lambda_m)^t \vartheta_{A_h})^t + \frac{S_h}{N_h} p (1_4 (\Phi \lambda_m)^t \vartheta_{P_h})^t \\ &\quad - (\Theta e v_h)^t \vartheta_{I_h} - g(t) \vartheta_{R_h} + \frac{S_m}{N_h} \left( (1_2 \beta_h \varepsilon)^t - (1_2 \lambda_h) 1_4^t \right) \vartheta_{S_m} - \frac{S_m}{N_h} \left( (\beta_h \varepsilon)^t \vartheta_{I_m} - \langle \lambda_h, \vartheta_{I_m} \rangle 1_4^t \right), \\ \dot{\vartheta}_{P_h} &= b - \frac{S_h}{N_h} (1_2 \lambda_m) 1_4^t \vartheta_{S_h} + (\mu_h + g(t)) \vartheta_{P_h} + \frac{S_h}{N_h} p (1_4 (\Phi \lambda_m)^t \vartheta_{P_h})^t + \frac{S_h}{N_h} (1-p) (1_4 (\Phi \lambda_m)^t \vartheta_{A_h})^t \\ &\quad - g(t) \vartheta_{R_h} + \frac{S_m}{N_h} \left( (1_2 \beta_h \varepsilon)^t - (1_2 \lambda_h) 1_4^t \right) \vartheta_{S_m} - \frac{S_m}{N_h} \left( (\beta_h \varepsilon)^t \vartheta_{I_m} - \langle \lambda_h, \vartheta_{I_m} \rangle 1_4^t \right), \\ \dot{\vartheta}_{I_h} &= b - \frac{S_h}{N_h} (1_2 \lambda_m) 1_4^t \vartheta_{S_h} + \frac{S_h}{N_h} (1-p) (1_4 (\Phi \lambda_m)^t \vartheta_{A_h})^t + \frac{S_h}{N_h} p (1_4 (\Phi \lambda_m)^t \vartheta_{P_h})^t - \gamma_h \vartheta_{R_h} \\ &\quad + (\mu_h I_d + \gamma_h + \delta_h) \vartheta_{I_h} + \delta_h 1_4^t \vartheta_{D_h} + \frac{S_m}{N_h} \left( (1_2 \beta_h \varepsilon)^t - (1_2 \lambda_h) 1_4^t \right) \vartheta_{S_m} - \frac{S_m}{N_h} \left( (\beta_h \varepsilon)^t \vartheta_{I_m} - \langle \lambda_h, \vartheta_{I_m} \rangle 1_4^t \right), \\ \dot{\vartheta}_{R_h} &= -\frac{S_h}{N_h} (1_2 \lambda_m) 1_4^t \vartheta_{S_h} + \frac{S_h}{N_h} (1-p) (1_4 (\Phi \lambda_m)^t \vartheta_{A_h})^t + \frac{S_h}{N_h} p (1_4 (\Phi \lambda_m)^t \vartheta_{P_h})^t - \rho_h 1_4^t \vartheta_{S_h} + (\mu_h + \rho_h) \vartheta_{R_h} \\ &\quad + \frac{S_m}{N_h} \left( (1_2 \lambda_h) 1_4^t \vartheta_{S_m} - \langle \lambda_h, \vartheta_{I_m} \rangle 1_4^t \right), \\ \dot{\vartheta}_{D_h} &= -c_h, \end{aligned} \right. \qquad (5.5) $$
and

$$
\left\{
\begin{aligned}
\dot{\vartheta}_{S_m} &= (\mu_m + 1_2 \lambda_h) \vartheta_{S_m} - \langle \lambda_h, \vartheta_{I_m} \rangle, \\
\dot{\vartheta}_{I_m} &= \delta_m \vartheta_{I_m} + \frac{S_h}{N_h} (1_2 \beta_m)^t \vartheta_{S_h} - \frac{S_h}{N_h} (\Phi \beta_m)^t ((1-p)\vartheta_{A_h} + p\vartheta_{P_h}),
\end{aligned}
\right.
\qquad (5.6)
$$

with $b = (0, -b_h, 0, -b_h)$. Systems (5.5) and (5.6) are associated with the boundary conditions or
transversality conditions, $\vartheta_z(T_f) = 0$, for all $z \in \mathcal{Z}$.

Furthermore, by (5.4), we have:

$$
\frac{\partial \mathcal{H}}{\partial \theta} = c_{\theta} \theta - \nu_{h} \sum_{k} A_{h}^{k,T} (\vartheta_{I_{h}}^{k,T} - \vartheta_{I_{h}}^{k,U}) - \nu_{h} \sum_{k} A_{h}^{k,U} (\vartheta_{I_{h}}^{k,T} - \vartheta_{I_{h}}^{k,U}) = 0,
$$

$$
\frac{\partial \mathcal{H}}{\partial \phi} = c_{\phi} \phi - \sum_{k} (1-p) \lambda_{m}^{k} S_{h} \vartheta_{A_{h}}^{k,U} - \sum_{k} p \lambda_{m}^{k} S_{h} \vartheta_{P_{h}}^{k,U} + \sum_{k} (1-p) \lambda_{m}^{k} S_{h} \vartheta_{A_{h}}^{k,T} + \sum_{k} p \lambda_{m}^{k} S_{h} \vartheta_{P_{h}}^{k,T} = 0,
$$

$$
\frac{\partial \mathcal{H}}{\partial g} = c_g g - b_h \left( (A_h^{r,T} + A_h^{r,U}) + (P_h^{r,T} + P_h^{r,U}) \right) - \sum_{k,l} (A_h^{k,l})^t (\vartheta_{A_h}^{k,l} - \vartheta_{R_h}^{k,l}) - \sum_{k,l} (P_h^{k,l})^t (\vartheta_{P_h}^{k,l} - \vartheta_{R_h}^{k,l}) = 0. \tag{5.7}
$$

Therefore, if $(\theta^*, \phi^*, g^*) \in \Omega$ is a solution of (5.2), then by (5.7), it is characterized by:

$$
\theta = \max \left\{ 0, \min \left\{ 1; \frac{\nu_h \sum_k (A_h^{k,T} + A_h^{k,U}) (\vartheta_{I_h}^{k,T} - \vartheta_{I_h}^{k,U})}{c_\theta} \right\} \right\},
$$

$$
g = \max \left\{ 0, \min \left\{ g_{\max}; \frac{b_h \left( (A_h^{r,T} + A_h^{r,U}) + (P_h^{r,T} + P_h^{r,U}) \right) + \sum_{k,l} [(A_h^{k,l})^t (\vartheta_{A_h}^{k,l} - \vartheta_{R_h}^{k,l}) + (P_h^{k,l})^t (\vartheta_{P_h}^{k,l} - \vartheta_{R_h}^{k,l})] }{c_g} \right\} \right\}.
$$

The state system (3.1), the adjoint systems (5.5) and (5.6), and the control characterization system (5.8) together form the optimality system that must be solved numerically. These systems operate as a feedback system in the following sense.

(1) State system input: A given controller serves as the input for the state system.

(2) Adjoint system input: The controller, along with the output of the state system, becomes the input for the adjoint system.

(3) Control update: The outputs of the state and adjoint systems are then used to update the controller via the control characterization system.

(4) Loop closure: The updated controller is fed back as the input for the state system, forming a “closed loop” system.

Since the state equations are governed by the initial conditions and the adjoint equations rely on the
final-time conditions, the optimality system cannot be solved straightforwardly by sweeping forward
in time alone. Instead, an iterative approach known as the forward–backward sweep method is applied
(see [26]). This method alternates between solving the state equations forward in time and the adjoint
equations backward in time, iteratively updating the control until convergence is achieved.
Figure 4 demonstrates the impact of an optimal intervention strategy on the number of deaths and clinical cases. The optimal control strategy strongly discourages the exposure of asymptomatic and paucisymptomatic infections to antimalarial treatment without curative intent (Figure 4A). This aligns with the findings of the global sensitivity analysis, which indicates that this parameter has a minimal impact on both mortality and the prevalence of resistant strains. Antimalarial treatment with curative intention is employed at relatively high levels during the initial years, followed by a gradual reduction toward the end of the control period. In contrast, the use of MDA is strongly recommended at a consistently high rate throughout the control period (Figure 4A). Furthermore, the reduction in antimalarial treatment is effectively offset by the increased use of MDA within the host population. This approach is likely driven by the assumption that the costs of both control strategies are equivalent in this scenario. The optimal intervention strategies markedly diminish the severity of the projected outbreaks, reducing both the total number of deaths (Figure 4B) and the prevalence of clinical cases (Figure 4C).

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Effect of intervention strategies on the model's dynamics. (A) The implementation of intervention strategies $\phi$, $\theta$, and $q_{mda}$ as functions of time $t$. (B) The total number of deaths over time. (C) The total number of clinical cases over time. Here, $\varepsilon_s = 10^{-3}$ and the control costs are $c_\phi = 10$, $c_\theta = 20$, and $c_g = 20$. Parameter values are the defaults in Table 2.</figcaption></figure>

**6. Discussion**

Antimalarial treatments and mass drug administration (MDA) are key strategies to alleviate the burden of malaria. However, the rise of parasite resistance to antimalarial drugs presents a significant challenge to effective malaria control and treatment. To assess the impact of combining antimalarial treatments and MDA for effective disease control, the human population is categorized into five distinct states: susceptible, asymptomatic infected (individuals with very mild or no symptoms who may later develop symptomatic infections), paucisymptomatic infected (individuals with very mild or no symptoms who do not progress to symptomatic infections), clinical cases, and recovered individuals. An infected human, whether symptomatic, paucisymptomatic, or asymptomatic, can be treated with antimalarial drugs either with curative intent or without. Meanwhile, the mosquito population is divided into two states: susceptible and infectious mosquitoes.
We use optimal control theory [26] to characterize an optimal strategy that significantly reduces the number of deaths and resistant clinical cases. Our formulation assumes a quadratic cost for the control effort. Overall, for a relatively low mutation rate from treated sensitive human infections to treated resistant human infections ($ε_s = 10^{-3}$), the optimal strategy strongly limits antimalarial treatments for asymptomatic and paucisymptomatic infections without curative intent. Curative treatments are applied heavily in the initial years, then taper off while MDA is consistently maintained at high levels throughout the control period (Figure 4).

In the scenario of a relatively high mutation rate ($ε_s = 10^{-1}$) (see Figure 5), the configuration closely resembles that of a low mutation rate (Figure 4). However, a higher mutation rate leads to a reduced deployment of antimalarial drugs, which is compensated by an increased reliance on MDA within the host population for effective control (Figures 4A and 5A). Intuitively, when the mutation rate is relatively low, the optimal strategy combines the use of both antimalarial treatment with curative intent and MDA, with a gradual decrease in antimalarial use over time, ensuring sustainable malaria control. In contrast, at higher mutation rates, the strategy prioritizes broader deployment of MDA while significantly reducing reliance on antimalarial drugs to minimize the risk of resistance developing.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Effect of intervention strategies on the model dynamics as in Figure 4 but with $ε_s = 10^{-1}$.</figcaption></figure>

Figure 6 assesses the impact of control costs on the model's dynamics. For a high cost of curative treatment and a low cost of MDA (as shown in Figure 6A–C), the control strategies favor higher reliance on MDA and reduced use of antimalarial treatments. This pattern is evident in society, where high treatment costs often prevent individuals from seeking care until they are already ill. As a result, implementing MDA on a large scale is crucial to effectively reduce both deaths and clinical cases over time.

In contrast, when the cost of curative treatments is lower than that of MDA (see Figure 6D–F), the control strategy advocates for a combination of both approaches: treating clinical cases with curative intent while also implementing a significant level of MDA to more effectively reduce deaths and clinical cases. However, the use of MDA is less extensive in this scenario compared with the previous one (see Figure 6A–C), likely due to the higher cost of MDA. It is also clear that the number of deaths and
clinical cases is lower when MDA is more affordable (Figure A2). It is essential to note that treating
infected individuals without curative intent in both scenarios is strongly discouraged to prevent the
emergence of resistance over time.

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

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Effect of the costs of intervention strategies on the model’s dynamics.</figcaption></figure>

Human malaria infections are characterized by the actual presence of parasites within an individual, but with the absence or very mild malaria-related symptoms, this remains a challenge for malaria control programs due to their significant influence on transmission dynamics [12, 13]. While some infected individuals never develop clinical symptoms and transit to the group of paucisymptomatic infections, others initially progress through an asymptomatic phase (i.e., without symptoms) before developing into symptomatic or clinical cases. Treating the infectious parasite reservoir in paucisymptomatic and asymptomatic individuals is a critical intervention strategy with a positive impact on reducing the malaria burden (Figures 4 and 5). The proportion of paucisymptomatic infections (p) within the host population significantly affects control strategies that combine MDA and antimalarial treatments (Figure 7). Specifically, a low proportion of paucisymptomatic infections results in a greater reliance on
MDA for effective control (Figure 7, Row 1). This is because a lower prevalence of paucisymptomatic infections corresponds to a higher proportion of asymptomatic infections that eventually progress to symptomatic cases. Consequently, antimalarial use remains relatively high during the initial years of the control period and continues, albeit at a lower rate, over the duration of the control period (Figure 7A). In contrast, when the proportion of paucisymptomatic infections is moderate (Figure 7, Row 2) or high (Figure 7, Row 3), the deployment of MDA is less intense compared with scenarios with low paucisymptomatic prevalence. Additionally, the use of antimalarials is robust during the initial years of the control period but decreases substantially afterward (Figure 7D,G).

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Effect of paucisymptomatic infections on the optimal control strategy.</figcaption></figure>

The proposed model assumes heterogeneity within both the human and mosquito populations. However, this presents a potential limitation, as incorporating factors such as the chronological age of human and mosquito populations and the time since infection can have significant epidemiological implications for malaria control strategies [13, 28–30]. According to previous studies, e.g., [28, 29], an age-structured model formulation within this context may indeed be interesting to investigate in detail but would deserve a dedicated study.

Another potential limitation of this study is the arbitrary choice of the balancing coefficients in the
objective function. Specifically, the value of the balancing cost lacks precise justification, introducing an element of uncertainty into the outcome of the optimal control function. Although we have deliberately chosen its value, we acknowledge that this parameter can significantly influence the resulting optimal control function.

# 7. Proof of main results

## 7.1. Proof of Theorem 3.1

Note that the linear operator $\mathcal{B}$ is the infinitesimal generator of a strongly continuous semigroup $\{e^{\mathcal{B}t}\}_t$ which preserves $\mathbb{R}_+^{20}$. Furthermore, the map $\mathcal{M}: \mathbb{R}^{20} \to \mathbb{R}^{20}$ is locally Lipschitz; therefore, classical results (e.g., [24, 31]) shows that for any $v_0 \in \mathbb{R}^{20}$, we can find a time $T = T(v_0) > 0$ such that (3.2) has a unique solution $v \in C([0, T)) \cap C^1((0, T))$ such that

$$ \mathcal{V}(t, v_0) := v(t) = e^{\mathcal{B}t}v_0 + \int_0^t e^{\mathcal{B}(t-\eta)} \mathcal{M}(v(\eta))d\eta, \forall t \in [0, T). $$

Moreover, statements below are valid.

(i) If $T < +\infty$, then $\lim_{t \to T^-} \|\mathcal{V}(t, v_0)\|_{\mathbb{R}^{20}} = +\infty$.

(ii) Let $t \in (0, T)$. If a sequence $\{u_n\}_n \subset \mathbb{R}^{20}$ such that $u_n \to v_0$ in $\mathbb{R}^{20}$ as $n \to \infty$, then $\mathcal{V}(\eta, u_n)$ exists for all $\eta \in [0, t]$, where $n$ is sufficiently large and $\sup_{\eta \in [0, t]} \|\mathcal{V}(\eta, u_n) - \mathcal{V}(\eta, v_0)\|_{\mathbb{R}^{20}} \to 0$, as $n \to \infty$.

We are now interested in the positiveness of the solutions. Note that $\mathcal{B}$ is resolvent positive, i.e., for all $\zeta$, with $\zeta > -\omega := \min\{\mu_h, \mu_m, \delta_m\} > 0$, $(\zeta I_d - \mathcal{B})$ is invertible and $(\zeta I_d - \mathcal{B})^{-1}\mathbb{R}_+^{20} \subset \mathbb{R}_+^{20}$. Moreover, for all $\eta > 0$, we can find $\zeta_\eta > 0$ such that

$$ \mathcal{M}(v(t)) + \zeta_\eta v(t) \geq 0_{\mathbb{R}^{20}}, \quad \forall v \in \mathbb{R}_+^{20} \cap \bar{B}(0, \eta), \forall t \in [0, T), \quad (7.1) $$

with $B(0, \eta)$ being the open ball in $\mathbb{R}^{20}$ centered at 0 with a radius $\eta$. Therefore, for all $v_0 \in \mathbb{R}_+^{20}$, $\mathcal{V}(t, v_0) \in \mathbb{R}_+^{20}$ for all $t \in [0, T)$, i.e.,

$$ \mathcal{V}(t, \mathbb{R}_+^{20}) \subset \mathbb{R}_+^{20}, \quad \forall t \in [0, T). $$

The statement above follows from well-known results, e.g., [34, Theorem 5.2.7, p. 226]. Inequality (7.1) should be understood in the following sense. We recall that the mapping $\mathcal{M}: \mathbb{R}^{20} \to \mathbb{R}^{20}$ is such that $\mathcal{M} = (\mathcal{M}_\ell)_{\ell=1,\dots,20}$. Thus, (7.1) expresses that, for all, $\ell = 1, \dots, 20$, $\mathcal{M}_\ell(v) + \zeta_\eta v_\ell \ge 0$, where $v = (v_\ell)_{\ell=1,\dots,20} \in \mathbb{R}^{20}$.

Next, for $v_0 \in \mathbb{R}_+^{20}$, let $\mathcal{V}(t, v_0) \in \mathbb{R}_+^{20}$, for all $t \in [0, T(v_0))$, be the solution of (3.2). Then

$$ \dot{N}_h(t) \le \Lambda_h - \mu_h N_h(t), \quad \dot{N}_m(t) \le \Lambda_m - \min(\mu_m, \delta_m^*) N_m(t), \quad \text{with} \quad \delta_m^* = \min(\delta_m^s, \delta_m^r), $$

i.e.,

$$ N_h(t) \le N_h(0)e^{-\mu_h t} + \frac{\Lambda_h}{\mu_h}(1 - e^{-\mu_h t}), $$
$$N_m(t) \leq N_m(0)e^{-\min(\mu_m, \delta_m^*)t} + \frac{\Lambda_m}{\min(\mu_m, \delta_m^*)}(1 - e^{-\min(\mu_m, \delta_m^*)t}).$$

From the case where $\forall t \in [0, T)$,

$$N_h(t) \leq \max\left(N_h(0), \frac{\Lambda_h}{\mu_h}\right), \quad N_m(t) \leq \max\left(N_m(0), \frac{\Lambda_h}{\min(\mu_m, \delta_m^*)}\right).$$

Consequently

$$\sup_{t \in [0, T(v_0))} \|\mathcal{V}(t, v_0)\| < \infty,$$

which ends the global well-posedness in $\mathbb{R}_+^{20}$.

It remains to show the bounded dissipativity, that is, a bounded subset $Q \subset \mathbb{R}_+^{20}$ exists such that, for each bounded subset $Q \subset \mathbb{R}_+^{20}$, $t_0 = t_0(Q) > 0$ exists such that $\mathcal{V}(t, Q) \subset Q$, for all $t \ge t_0$. Let $v_0 \in \mathbb{R}_+^{20}$, and let $\mathcal{V}(t, v_0) = (S(t), A_h(t), P_h(t), I_h(t), R_h(t), S_m(t), I_m(t))$ be the corresponding semiflow. Setting

$$Q = \left\{ \begin{array}{l} (S_h, A_h, P_h, I_h, R_h, S_m, I_m) \in \mathbb{R}_+^{20} : \\ S_h + S_m + \|A_h\| + \|P_h\| + \|I_h\| + \|R_h\| + \|I_m\| \le \frac{\Lambda_h}{\mu_h} + \frac{\Lambda_m}{\min(\mu_m, \delta_m^*)} \end{array} \right\},$$

we easily find the bounded dissipativity of the semiflow.

To end with the proof of Theorem 3.1, it remains to show that the semiflow $\{\mathcal{V}(t, \cdot)\}_t$ is asymptotically smooth, i.e., for any closed, bounded, and positively invariant set $Q \subset \mathbb{R}_+^{20}$, a compact set $K \subset \mathbb{R}_+^{20}$ exists such that $d(\mathcal{V}(t, Q), K) \to 0$, as $t \to \infty$, where $d$ is the Hausdorff semi distance, see e.g., [32]. Note that

$$\mathcal{V}(t, v_0) = e^{\beta t} v_0 + \int_0^t e^{\beta(t-\eta)} \mathcal{M}(\mathcal{V}(\eta, v_0)) d\eta.$$

Since $\{e^{\beta t}\}_t$ is strictly contractive and $\mathcal{M}$ is compact, the lemma follows from the results established in [33].

## 7.2. Proof of Lemma 3.3

By Lemma 3.2, we can employ the results in [22] to establish the first part of the lemma.

For the second part of the lemma, note that

$$(-\hat{B})^{-1} = \begin{pmatrix} (\mu_h + g + v_h)^{-1}\mathbf{I}_d & 0_{4 \times 4} & 0_{4 \times 4} & 0_{4 \times 2} \\ 0_{4 \times 4} & (\mu_h + g)^{-1}\mathbf{I}_d & 0_{4 \times 4} & 0_{4 \times 2} \\ v_h(\mu_h + g + v_h)^{-1}(\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1}\mathbf{\Theta}\mathbf{e} & 0_{4 \times 4} & (\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1} & 0_{4 \times 2} \\ 0_{2 \times 4} & 0_{4 \times 4} & 0_{4 \times 4} & \delta_m^{-1} \end{pmatrix}$$

From this, $\hat{M}(-\hat{B})^{-1}$ can be written as

$$\hat{M}(-\hat{B})^{-1} = \begin{pmatrix} 0_{12 \times 12} & V \\ U & 0_{2 \times 2} \end{pmatrix},$$

where $U$ and $V$ are, respectively, 2 × 12 and 12 × 2 matrices such that

$$U = \frac{S_m^0}{S_h^0} \left( (\mu_h + g + v_h)^{-1} \beta_h \varepsilon [\mathbf{I}_d + v_h(\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1} \mathbf{\Theta} \mathbf{e}] \beta_h \varepsilon (\mu_h + g_h)^{-1} \mathbf{I}_d \quad \beta_h \varepsilon (\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1} \right),$$
$$V = \begin{pmatrix} (1-p)\Phi\beta_m\delta_m^{-1} \\ p\Phi\beta_m\delta_m^{-1} \\ 0_{4\times2} \end{pmatrix}.$$

Therefore,

$$(\hat{M}(-\hat{B})^{-1})^2 = \operatorname{diag}(UV, VU),$$

where

$$(r(\hat{M}(-\hat{B})^{-1}))^2 = r(\hat{M}(-\hat{B})^{-1})^2 = \max\{r(UV), r(VU)\}.$$

Finally, we have

$$r(UV) = r(VU) = \frac{S_h^0}{S_h^0} r(\mathcal{T}_h \mathcal{T}_m),$$

with $\mathcal{T}_h$ and $\mathcal{T}_m$ introduced by (3.3).

### 7.3. Proof of Theorem 3.4

The different components of $E_h^*$ are given by

$$ \begin{gathered} \begin{pmatrix} S_h^* \\ R_h^* \end{pmatrix} = \begin{pmatrix} \mu_h + \mathbf{1}_2 \lambda_m^* & -\rho_h \mathbf{1}_4 \\ -[\gamma_h(\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1} \mathbf{\Theta} \mathbf{e} v_h (v_h + \mu_h + g)^{-1}(1-p) + g(g + \mu_h)^{-1} \mathbf{I}_d p \\ +g(v_h + \mu_h + g)^{-1}(1-p)] \Phi \lambda_m^* & (\mu_h + \rho_h) \mathbf{I}_d \end{pmatrix}^{-1} \begin{pmatrix} \Lambda_h \\ 0_4 \end{pmatrix} \\ := f(\lambda_m^*). \end{gathered} $$

Set

$$Z = [\gamma_h(\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1} \mathbf{\Theta} \mathbf{e} v_h (v_h + \mu_h + g)^{-1}(1-p) + g(g + \mu_h \mathbf{I}_d)^{-1} p + g(v_h + \mu_h + g)^{-1}(1-p)] \Phi,$$

then

$$\begin{pmatrix} S_h^* \\ R_h^* \end{pmatrix} := f(\lambda_m^*) = \begin{pmatrix} \Lambda_h[\mu_h + (\mathbf{1}_2 - \rho_h \mathbf{1}_4(\mu_h + \rho_h)^{-1} \mathbf{I}_d Z)\lambda_m^*]^{-1} \\ \Lambda_h(\mu_h + \rho_h)^{-1} \mathbf{I}_d Z[\mu_h + (\mathbf{1}_2 - \rho_h \mathbf{1}_4(\mu_h + \rho_h)^{-1} \mathbf{I}_d Z)\lambda_m^*]^{-1} \lambda_m^* \end{pmatrix}.$$

Consequently

$$ \begin{pmatrix} A_h^* \\ P_h^* \\ I_h^* \end{pmatrix} = \begin{pmatrix} (v_h + \mu_h + g)^{-1} \mathbf{I}_d (1-p) \Phi \lambda_m^* & 0_{4 \times 4} \\ (g + \mu_h)^{-1} \mathbf{I}_d p \Phi \lambda_m^* & 0_{4 \times 4} \\ (\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1} \mathbf{\Theta} \mathbf{e} v_h (v_h + \mu_h + g)^{-1} \mathbf{I}_d (1-p) \Phi \lambda_m & 0_{4 \times 4} \end{pmatrix} f(\lambda_m^*). $$

i.e.,

$$ \begin{pmatrix} A_h^* \\ P_h^* \\ I_h^* \end{pmatrix} = \begin{pmatrix} (v_h + \mu_h + g)^{-1}(1-p)\Phi\Lambda_h[\mu_h + (\mathbf{1}_2 - \rho_h \mathbf{1}_4(\mu_h + \rho_h)^{-1}\mathbf{I}_d Z)\lambda_m^*]^{-1}\lambda_m^* \\ (g + \mu_h)^{-1}\mathbf{I}_d p\Phi\Lambda_h[\mu_h + (\mathbf{1}_2 - \rho_h \mathbf{1}_4(\mu_h + \rho_h)^{-1}\mathbf{I}_d Z)\lambda_m^*]^{-1}\lambda_m^* \\ (\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1} \mathbf{\Theta} \mathbf{e} v_h (v_h + \mu_h + g)^{-1} \mathbf{I}_d (1-p)\Phi\Lambda_h[\mu_h + (\mathbf{1}_2 - \rho_h \mathbf{1}_4(\mu_h + \rho_h)^{-1} \mathbf{I}_d Z)\lambda_m^*]^{-1} \lambda_m^* \end{pmatrix}. $$
At the equilibrium, the total population of human is written as

$$N_h^* = (1 \ 1_4) \begin{pmatrix} S_h^* \\ R_h^* \end{pmatrix} + (1_4 \ 1_4 \ 1_4) \begin{pmatrix} A_h^* \\ P_h^* \\ I_h^* \end{pmatrix}.$$

We have

$$N_h^* = (1 \ 1_4) f(\lambda_m^*) + (1_4 \ 1_4 \ 1_4) \begin{pmatrix} (v_h + \mu_h + g)^{-1} I_d (1-p) \Phi \lambda_m^* & 0_{4 \times 4} \\ (g + \mu_h)^{-1} I_d p \Phi \lambda_m^* & 0_{4 \times 4} \\ (\mu_h I_d + \gamma_h + \delta_h)^{-1} \Theta e v_h (v_h + \mu_h + g)^{-1} I_d (1-p) \Phi \lambda_m^* & 0_{4 \times 4} \end{pmatrix} f(\lambda_m^*) \\ := g(\lambda_m^*) f(\lambda_m^*).$$

With this, the two forces of infection are written as follows

$$\lambda_h^* = \frac{\beta_h}{g(\lambda_m^*) f(\lambda_m^*)} \varepsilon \begin{pmatrix} 1 & 1 & 1 \end{pmatrix} \begin{pmatrix} (v_h + \mu_h + g)^{-1} I_d (1-p) \Phi \lambda_m^* & 0_{4 \times 4} \\ (g + \mu_h)^{-1} I_d p \Phi \lambda_m^* & 0_{4 \times 4} \\ (\mu_h I_d + \gamma_h + \delta_h)^{-1} \Theta e v_h (v_h + \mu_h + g)^{-1} I_d (1-p) \Phi \lambda_m^* & 0_{4 \times 4} \end{pmatrix} f(\lambda_m^*) \\ := J(\lambda_m^*), \\ \lambda_m^* = \frac{1}{g(\lambda_m^*) f(\lambda_m^*)} \beta_m I_m^*.$$

$$J(\lambda_m^*) = \frac{\beta_h}{g(\lambda_m^*) f(\lambda_m^*)} \varepsilon \begin{pmatrix} 1 & 1 & 1 \end{pmatrix} \begin{pmatrix} (v_h + \mu_h + g)^{-1} (1-p) \Phi \Lambda_h [\mu_h + (\mathbf{1}_2 - \rho_h \mathbf{1}_4 (\mu_h + \rho_h)^{-1} \mathbf{I}_d Z) \lambda_m^*]^{-1} \lambda_m^* \\ (g + \mu_h)^{-1} \mathbf{I}_d p \Phi \Lambda_h [\mu_h + (\mathbf{1}_2 - \rho_h \mathbf{1}_4 (\mu_h + \rho_h)^{-1} \mathbf{I}_d Z) \lambda_m^*]^{-1} \lambda_m^* \\ (\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1} \Theta e v_h (v_h + \mu_h + g)^{-1} \mathbf{I}_d (1-p) \Phi \Lambda_h [\mu_h + (\mathbf{1}_2 - \rho_h \mathbf{1}_4 (\mu_h + \rho_h)^{-1} \mathbf{I}_d Z) \lambda_m^*]^{-1} \lambda_m^* \end{pmatrix}.$$

Furthermore, replacing $\lambda_h^*$ in the two last equations of System (3.1) yields

$$\left\{ \begin{array}{l} S_m^* = (\mu_m + \mathbf{1}_2 \lambda_h^*)^{-1} \Lambda_m = (\mu_m + \mathbf{1}_2 J(\lambda_m^*))^{-1} \Lambda_m, \\ I_m^* = \delta_m^{-1} \lambda_h^* S_m = \delta_m^{-1} J(\lambda_m^*) (\mu_m + \mathbf{1}_2 J(\lambda_m^*))^{-1} \Lambda_m. \end{array} \right.$$

Consequently, the existence of an endemic equilibrium is given by the solution of the following fixed-point equation

$$\mathbb{R}^2 \ni \lambda_m^* = \frac{1}{g(\lambda_m^*) f(\lambda_m^*)} \beta_m \delta_m^{-1} J(\lambda_m^*) (\mu_m + \mathbf{1}_2 J(\lambda_m^*))^{-1} \Lambda_m^* := F(\lambda_m^*) \in \mathbb{R}^2,$$

The system above is then rewritten as

$$\lambda_m^* = G(\lambda_m^*) \beta_m \delta_m^{-1} \beta_h \varepsilon \{(v_h + \mu_h + g)^{-1}(1-p) + (g+\mu_h)^{-1}p + (\mu_h I_d + \gamma_h + \delta_h)^{-1} \Theta e v_h (v_h + \mu_h + g)^{-1}(1-p)\} \Phi \lambda_m^*,$$
that is,

$$\mathcal{T}_m \mathcal{T}_h \lambda_m^* = \frac{1}{G(\lambda_m^*)} \lambda_m^*, \qquad (7.2)$$
where

$$G : \lambda_m^* \in \mathbb{R}^2 \mapsto G(\lambda_m^*) = \frac{\Lambda_m G_1(\lambda_m^*)}{G_2(\lambda_m^*) G_3(\lambda_m^*)} \in \mathbb{R},$$

with

$$ 
\begin{aligned}
G_1(\lambda_m^*) &= \mu_h + (\mathbf{1}_2 - \rho_h \mathbf{1}_4(\mu_h + \rho_h)^{-1} \mathbf{I}_d Z) \lambda_m^*, \\
G_2(\lambda_m^*) &= 1 + \{\mathbf{1}_4[(v_h + \mu_h + g)^{-1}(1-p) + (g+\mu_h)^{-1} \mathbf{I}_d p]\Phi\} \lambda_m^* \\
&\quad + \{\mathbf{1}_4(\mu_h \mathbf{I}_d + \gamma_h + \delta_h)^{-1} \Theta e v_h (v_h + \mu_h + g)^{-1} \mathbf{I}_d (1-p) \Phi + \mathbf{1}_4(\mu_h + \rho_h)^{-1} \mathbf{I}_d Z\} \lambda_m^*, \\
G_3(\lambda_m^*) &= \Lambda_h \mu_m G_2(\lambda_m^*) + \mathbf{1}_2 \beta_h \varepsilon \{(v_h + \mu_h + g)^{-1}(1-p) + (g+\mu_h)^{-1} p + (\mu_h \mathbf{I}_d \\
&\quad + \gamma_h + \delta_h)^{-1} \Theta e v_h (v_h + \mu_h + g)^{-1}(1-p)\} \Phi \lambda_m^*.
\end{aligned}
$$

Since $\mathcal{T}_m\mathcal{T}_h$ is a positive and irreducible matrix, by the Perron-Frobenius theorem, we have

$$ 
\begin{aligned}
\frac{1}{G(\lambda_m^*)} &= r(\mathcal{T}_m\mathcal{T}_h), \\
\lambda_m^* &= \kappa\omega,
\end{aligned}
$$

where $\omega = (\omega^s, \omega^r)$ is the eigenvector, with $\|\omega\| = 1$, associated with the principal eigenvalue $r(\mathcal{T}_m\mathcal{T}_h)$ of the matrix $\mathcal{T}_m\mathcal{T}_h$, and $\kappa$ is a positive constant to be determined.

The constant $\kappa$ is then the solution of the equation $\frac{1}{G(\lambda_m^*)} = r(\mathcal{T}_m\mathcal{T}_h)$, i.e., $\frac{1}{G(\kappa\omega)} = r(\mathcal{T}_m\mathcal{T}_h)$, which is rewritten as

$$ \Lambda_m r(\mathcal{T}_m\mathcal{T}_h) \frac{G_1(\kappa\omega)}{G_2(\kappa\omega)G_3(\kappa\omega)} = 1. $$

The equality above is equivalent to

$$ L(\kappa) := \Lambda_m r(\mathcal{T}_m\mathcal{T}_h) \frac{(\mu_h + \chi_1\kappa)}{(1 + \chi_2\kappa)(\Lambda_h\mu_m + (\chi_2 + \Lambda_h\mu_m\omega^s + \chi_3)\kappa)} = 1, \quad (7.3) $$

with

$$ 
\begin{aligned}
\chi_1 &= \omega^r (1 - a_s) + \omega^s (1 - a_r), \\
\chi_2 &= \omega^r (b_s + a_s) + \omega^s (b_r + a_r), \\
\chi_3 &= \omega^s (c_s + d_s) + \omega^r (c_r + d_r),
\end{aligned}
$$

and

$$ 
\begin{aligned}
a_j &= \frac{v_h \theta \gamma_h^{j,T} \rho_h}{(\mu_h + \gamma_h^{j,T} + \delta_h^{j,T})(v_h + \mu_h)(\mu_h + \rho_h)} + \frac{v_h \gamma_h^{j,U} \rho_h (1 - \theta)}{(\mu_h + \gamma_h^{j,U} + \delta_h^{j,U})(v_h + \mu_h)(\mu_h + \rho_h)}, \\
b_j &= \frac{\theta v_h (1 - p)}{(\mu_h + \gamma_h^{j,T} + \delta_h^{j,T})(v_h + \mu_h + g)} + \frac{(1 - \theta) v_h (1 - p)}{(\mu_h + \gamma_h^{j,U} + \delta_h^{j,U})(v_h + \mu_h + g)} + \frac{(1 - p)}{(v_h + \mu_h + g)} + \frac{p}{(g + \mu_h)}, \\
c_j &= \left( \frac{(1 - p) \phi}{(v_h + \mu_h + g)} + \frac{p \phi}{(g + \mu_h)} + \frac{\theta (1 - p) v_h}{(\mu_h + \gamma_h^{j,T} + \delta_h^{j,T})(v_h + \mu_h + g)} \right) ((1 - \varepsilon_j) \beta_h^s + \varepsilon_j \beta_h^r), \\
d_j &= \beta_h^j \left( \frac{(1 - p) v_h (1 - \theta)}{(\mu_h + \gamma_h^{j,U} + \delta_h^{j,U})(v_h + \mu_h + g)} + \frac{(1 - p)(1 - \phi)}{(v_h + \mu_h + g)} + \frac{p(1 - \phi)}{(g + \mu_h)} \right).
\end{aligned}
$$
Since $L(0) = \Lambda_m r(\mathcal{T}_m \mathcal{T}_h) \frac{\mu_h}{\Lambda_h \mu_m} = \frac{S_m^0}{S_h^0} r(\mathcal{T}_m \mathcal{T}_h) = \mathcal{T}_0^2$, equality (7.3) is rewritten as $\mathcal{T}_0^2$-parametric equation

$$L(\mathcal{T}_0^2, \kappa) = 1, \qquad (7.4)$$

such that $L(1, 0) = 1$.

Then, using the implicit function theorem, we get

$$\left. \frac{d\kappa}{d\mathcal{T}_0^2} \right|_{(\mathcal{T}_0^2=1, \kappa=0)} = -\frac{\partial_{\mathcal{T}_0^2} L(1,0)}{\partial_{\kappa} L(1,0)}.$$

We have

$$\partial_{\mathcal{T}_0^2} L(1, 0) = 1,$$

and

$$\partial_{\kappa} L(1, 0) = \frac{\Lambda_m S_h^0}{S_m^0 (\Lambda_h \mu_m)^2} K,$$

with

$$K = \Lambda_h \mu_m \chi_1 - (2\chi_2 + \Lambda_h \mu_m \omega^s + \chi_3) \mu_h.$$

A backward bifurcation occurs at $\mathcal{T}_0^2 = 1$ if and only if

$$\left. \frac{d\kappa}{d\mathcal{T}_0^2} \right|_{(\mathcal{T}_0^2=1, \kappa=0)} < 0$$

which is equivalent to $\partial_{\kappa} L(1, 0) > 0$, that is, $K > 0$.

Similarly, a forward bifurcation occurs at $\mathcal{T}_0^2 = 1$ if and only if

$$\left. \frac{d\kappa}{d\mathcal{T}_0^2} \right|_{(\mathcal{T}_0^2=1, \kappa=0)} > 0$$

which is equivalent to $\partial_{\kappa} L(1, 0) < 0$, that is, $K < 0$.

Furthermore, for all $\mathcal{T}_0$, we find that $L(\mathcal{T}_0^2, 0) = \mathcal{T}_0^2$ and $\lim_{\kappa \to \infty} L(\mathcal{T}_0^2, \kappa) = 0$. Therefore, if $\mathcal{T}_0^2 > 1$, then at least one solution $\kappa$ to the equation $L(\mathcal{T}_0^2, \kappa) = 1$ exists. Next, assume that $\mathcal{T}_0^2 = 1$ and $K > 0$. Then, on one hand, we have $L(1, 0) = 1$ and $\partial_{\kappa} L(1, 0) > 0$. Hence, we obtain $L(1, \kappa_0) > 1$ for any $\kappa_0 > 0$ that is sufficiently small. Consequently, at least one $\kappa > 0$ exists such that $L(1, \kappa) = 1$.

Finally, the uniqueness of a positive $\kappa$ satisfying $L(k) = 1$ is obtained if the function $\kappa \mapsto L(\mathcal{T}_0^2, \kappa)$ is strictly decreasing on $\mathbb{R}_+$. Indeed, since $L(\mathcal{T}_0^2, 0) = \mathcal{T}_0^2$, we readily see that a unique $\kappa > 0$ exists such that $L(\mathcal{T}_0^2, \kappa) = 1$. We also note that a forward bifurcation occurs in this configuration since here, we have $K < 0$.

## Use of AI tools declaration

The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.
**Acknowledgments**

This work received support from the African Consortium in Modeling for Effective Vector Control (ACoMVeC) project INV-047049. The sponsors had no role in study design, data collection and analysis, decision to publish or manuscript preparation.

The authors would like to thank the anonymous referees for the valuable comments which helped to improve the overall quality of the manuscript.

**Conflict of interest**

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

**References**

1. A. Hamilton, F. Haghpanah, M. Hasso-Agopsowicz, I. Frost, G. Lin, E. Schueller, et al., Modeling of malaria vaccine effectiveness on disease burden and drug resistance in 42 African countries, *Commun. Med.*, **3** (2023), 144. https://doi.org/10.1038/s43856-023-00373-y

2. J. Muwanguzi, G. Henriques, P. Sawa, T. Bousema, C. J. Sutherland, K. B. Beshir, Lack of K13 mutations in Plasmodium falciparum persisting after artemisinin combination therapy treatment of Kenyan children, *Malar. J.*, **15** (2016), 1–6. https://doi.org/10.1186/s12936-016-1095-y

3. P. Patel, P. K. Bharti, D. Bansal, N. A. Ali, R. K. Raman, P. K. Mohapatra, et al., Prevalence of mutations linked to antimalarial resistance in Plasmodium falciparum from Chhattisgarh, Central India: A malaria elimination point of view, *Sci. Rep.*, **7** (2017), 16690. https://doi.org/10.1038/s41598-017-16866-5

4. C. W. Tarama, H. Soré, M. Siribié, S. Débé, R. Kinda, A. Ganou, et al., Plasmodium falciparum drug resistance-associated mutations in isolates from children living in endemic areas of Burkina Faso, *Malar. J.*, **22** (2023), 213. https://doi.org/10.1186/s12936-023-04645-9

5. A. G. B. Tuedom, E. M. Sarah-Matio, C. E. E. Moukoko, B. L. Feufack-Donfack, C. N. Maffo, A. N. Bayibeki, et al., Antimalarial drug resistance in the Central and Adamawa regions of Cameroon: Prevalence of mutations in P. falciparum crt, Pfmdr1, Pfdhfr and Pfdhps genes, *PLoS One*, **16** (2021), e0256343. https://doi.org/10.1371/journal.pone.0256343

6. M. I. Veiga, S. K. Dhingra, P. P. Henrich, J. Straimer, N. Gnädig, A. C. Uhlemann, et al., Globally prevalent PfMDR1 mutations modulate Plasmodium falciparum susceptibility to artemisinin-based combination therapies, *Nat. Commun.*, **7** (2016), 11553. https://doi.org/10.1038/ncomms11553

7. World Health Organization, *Malaria Microscopy Quality Assurance Manual*, World Health Organization, Geneva, 2nd edition, 2016.

8. T. D. Nguyen, T. N.-A. Tran, D. M. Parker, N. J. White, M. F. Boni, Antimalarial mass drug administration in large populations and the evolution of drug resistance, *PLoS Global Public Health*, **3** (2023), e0002200. https://doi.org/10.1371/journal.pgph.0002200
9. J. C. Koella, R. Antia. Epidemiological models for the spread of anti-malarial resistance, *Malar. J.*, **2** (2003), 1–11. https://doi.org/10.1186/1475-2875-2-3

10. P. Pemberton-Ross, N. Chitnis, E. Pothin, T. A. Smith, A stochastic model for the probability of malaria extinction by mass drug administration, *Malar. J.*, **16** (2017), 1–9. https://doi.org/10.1186/s12936-017-2010-x

11. J. M. Tchuenche, C. Chiyaka, D. Chan, A. Matthews, G. Mayer, A mathematical model for antimalarial drug resistance, *Math. Med. Biol. J. IMA*, **28** (2011), 335–355. https://doi.org/10.1093/imammb/dqq017

12. D. D. Laishram, P. L. Sutton, N. Nanda, V. L. Sharma, R. C. Sobti, J. M. Carlton, et al., The complexities of malaria disease manifestations with a focus on asymptomatic malaria, *Malar. J.*, **11** (2012), 1–15. https://doi.org/10.1186/1475-2875-11-29

13. R. Telfils, A. Y. Dossou, A. Djènontin, E. Adimi, R. Akoho, J. Bailly, et al., Dynamics of submicroscopic and microscopic asymptomatic malaria infection and associated factors: A longitudinal study in South Benin, *PLoS One*, **19** (2024), e0311217. https://doi.org/10.1371/journal.pone.0311217

14. L. S. Pontryagin, V. G. Boltyanskii, *Mathematical Theory of Optimal Processes*, John Wiley & Sons, 1st edition, 1962. https://doi.org/10.1201/9780203749319

15. N. Chitnis, J. M. Hyman, J. M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, *Bull. Math. Biol.*, **70** (2008), 1272–1296. https://doi.org/10.1007/s11538-008-9299-0

16. A. Mahamar, M. J. Smit, K. Sanogo, Y. Sinaba, S. M. Niambele, A. Sacko, et al., Artemether-lumefantrine with or without single-dose primaquine and sulfadoxine-pyrimethamine plus amodiaquine with or without single-dose tafenoquine to reduce Plasmodium falciparum transmission: A phase 2, single-blind, randomised clinical trial in Ouelessebougou, Mali, *Lancet Microbe*, **5** (2024), 633–644. https://doi.org/10.1016/S2666-5247(24)00023-5

17. L. Cirera, B. Galatas, S. Alonso, K. Paaijmans, M. Mamuquele, H. Martí-Soler, et al., Moving towards malaria elimination in southern Mozambique: Cost and cost-effectiveness of mass drug administration combined with intensified malaria control, *PLoS One*, **15** (2020), e0235631. https://doi.org/10.1371/journal.pone.0235631

18. S. S. Kyaw, G. Delmas, T. L. Drake, O. Celhay, W. Pan-Ngum, S. Pukrittayakamee, et al., Estimating the programmatic cost of targeted mass drug administration for malaria in Myanmar, *BMC Public Health*, **21** (2021), 826. https://doi.org/10.1186/s12889-021-10842-5

19. J. O. Yukich, C. Scott, K. Silumbe, B. A. Larson, A. Bennett, T. P. Finn, et al., Cost-effectiveness of focal mass drug administration and mass drug administration with dihydroartemisinin-piperaquine for malaria prevention in Southern Province, Zambia: Results of a Community-Randomized Controlled Trial, *Am. J. Trop. Med. Hyg.*, **103**(2_Suppl) (2020), 46–53. https://doi.org/10.4269/ajtmh.19-0661
20. J. K. Hale, J. P. LaSalle, M. Slemrod, Theory of a general class of dissipative processes, *J. Math. Anal. Appl.*, **39** (1972), 177–191.

21. O. A. Ladyzhenskaya, On the determination of minimal global attractors for the Navier-Stokes and other partial differential equations, *Russ. Math. Surv.*, **42** (1987), 25–60.
https://doi.org/10.1070/RM1987v042n06ABEH001503

22. H. R. Thieme, Spectral bound and reproduction number for infinite-dimensional population structure and time heterogeneity, *SIAM J. Appl. Math.*, **70** (2009), 188–211.
https://doi.org/10.1137/080732870

23. P. Magal, Compact attractors for time-periodic age-structured population models, *Electron. J. Dif-fer. Equations (EJDE)*, **2001** (2001), 1–35.

24. A. Pazy, *Semigroups of Linear Operators and Applications to Partial Differential Equations*, Springer-Verlag Berlin and Heidelberg GmbH & Co. K, New York Berlin Heidelberg, 1983.

25. W. Fleming, R. Rishel, *Deterministic and Stochastic Optimal Control*, Springer New York, New York, NY, 1975.

26. S. Lenhart, J. T. Workman, *Optimal Control Applied to Biological Models*, Chapman and Hall/CRC, New York, May 2007. https://doi.org/10.1201/9781420011418

27. J. Macki, A. Strauss, *Introduction to Optimal Control Theory*, Springer, New York Berlin Heidel-berg, corr print edition edition, 1981.

28. Q. Richard, M. Choisy, T. Lefèvre, R. Djidjou-Demasse, Human-vector malaria transmission model structured by age, time since infection and waning immunity, *Nonlinear Anal. Real World Appl.*, **63** (2022), 103393. https://doi.org/10.1016/j.nonrwa.2021.103393

29. Q. Richard, M. Choisy, T. Lefèvre, R. Djidjou-Demasse, On the necessity of accounting for age structure in human malaria transmission modeling, *Math. Biosci.*, **378** (2024), 109319.
https://doi.org/10.1016/j.mbs.2024.109319

30. B. M. Somé, E. Guissou, D. F. Da, Q. Richard, M. Choisy, K. B. Yameogo, et al., Mosquito ageing modulates the development, virulence and transmission potential of pathogens, *Proc. R. Soc. B Biol. Sci.*, **291** (2024), 20232097. https://doi.org/10.1098/rspb.2023.2097

31. T. Cazenave, A. Haraux, *An Introduction to Semilinear Evolution Equations*, Oxford University Press, 1998.

32. J. K. Hale, *Asympstotic Behavior of Dissipative Systems*, American Mathematical Society, 1988.

33. G. F. Webb, Compactness of bounded trajectories of dynamical systems in infi-nite dimensional spaces, *Proc. R. Soc. Edinburgh Sec. A: Math.*, **84** (1979), 19–33.
https://doi.org/10.1017/S0308210500016930

34. P. Magal, S. Ruan, *Theory and Applications of Abstract Semilinear Cauchy Problems*, Manhattan, NY, USA, Springer International Publishing, 2018. https://doi.org/10.1007/978-3-030-01506-0
**Appendix: Supplementary figure**

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

**Figure A1.** Sensitivity indices of disease-induced cumulative deaths (left) and resistant cases (right). The shaded portions of the bars represent the main indices (the effect of each factor alone), while the complete bars, including both the shaded and unshaded parts, represent the total effect of each factor, accounting for interactions with all other factors.

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

**Figure A2.** Disease-induced deaths and clinical cases under the control functions shown in Figure 6.

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

AIMS Press

©2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)