A mathematical modeling and optimal control analysis of the effect of treatment-seeking behaviors on the spread of malaria

Malaria, an infectious disease transmitted by mosquitoes and caused by the Plasmodium parasite, poses a significant global public health challenge, especially in areas lacking modern medical infrastructure. Traditional medicine often serves as either a primary or complementary treatment avenue. This study introduces a novel deterministic model that considers the impact of treatment seeking-behaviors on malaria transmission dynamics. Expanding upon the existing model, we incorporate distinct groups: individuals seeking treatment at health facilities and those self-treating with traditional remedies, which lack clinical validation. The study employs mathematical techniques for a comprehensive analysis of the model, including positivity, boundedness, existence and uniqueness, equilibrium, reproduction number, sensitivity, optimal control, and numerical simulations performed using MATLAB and the fourth-order Runge-Kutta method. Furthermore, we explore three time-dependent optimal control variables: antimalarial drug treatment, personal protective measures like ITNs, and promoting awareness to discourage inappropriate traditional medicine usage, all aimed at reducing disease transmission. Sensitivity analysis helps identify key parameters affecting malaria dynamics. Notably, increased utilization of health facilities for treatment significantly reduces the basic reproduction number, highlighting the importance of effective healthcare interventions. Numerical simulations underscore the vital role of treating infected individuals at health facilities in malaria eradication efforts. Optimal control analysis suggests that a combination of the three control strategies is most effective in combating malaria. This provides insights for public health policies to address the risk factors of using clinically not validated traditional medicine in malaria-endemic areas.

1 Introduction

Malaria remains a significant global public health challenge, particularly in endemic regions [1]. According to the World Malaria Report (2023), the estimated number of malaria cases worldwide reached 249 million in 2022, marking an increase of 16 million from the pre-pandemic level of 233 million in 2019 [2]. Malaria-related mortality has also remained high, with 608,000 deaths reported in 2022. The disease burden is disproportionately concentrated in the WHO African Region, accounting for approximately 94% of global cases and 96% of malaria deaths. Malaria is caused by Plasmodium parasites, which are transmitted to humans through the bites of infected female Anopheles mosquitoes [2].¹ Of the five Plasmodium species known to infect humans—Plasmodium falciparum, Plasmodium vivax, Plasmodium malariae, Plasmodium ovale, and Plasmodium knowlesi—P. falciparum and P. vivax are the most virulent and contribute to the highest morbidity and mortality rates [3].

The transmission dynamics of malaria are influenced by various factors, including human treatment-seeking behavior [4]. Effective and timely treatment plays a critical role in mitigating malaria transmission and reducing its societal impact [5, 6]. Historically, traditional medicine has played a crucial role in malaria treatment, particularly in regions with limited access to modern healthcare [7, 8]. It is estimated that between 25% and 75% of individuals in endemic areas rely on traditional remedies as either a primary or complementary treatment [9]. However, traditional medicine use presents substantial challenges, including a lack of rigorous clinical validation, unpredictable dosage variations, and concerns regarding safety and efficacy [10, 11]. The authors in [12] found that the main challenge with herbal malaria treatment is inappropriate dosage specifications. Unlike standardized pharmaceutical treatments, herbal malaria therapies are often prepared and administered based on empirical knowledge rather than precise scientific formulations. The absence of standardized dosing regimens increases the risk of treatment failure, severe complications, and, potentially, the emergence of antimalarial drug resistance, which is currently a major challenge to malaria control [9].

Mathematical modeling serves as a powerful tool for understanding infectious disease dynamics, evaluating intervention strategies, and informing public health policies [13]. Since the early work of Sir Ronald Ross in 1911 [14], as the disease burden increased, numerous malaria transmission models have been developed [15–18], incorporating factors such as latency periods, human–mosquito interactions, immunity, and heterogeneous susceptibility. The study in [16] presented a mathematical model for malaria dynamics, considering blood transfusion transmission and saturated treatment function to delay in administering anti-malaria drug treatment transmission in the population. They investigated the role of saturated treatment on the dynamics of malaria and verified that delay in treatment is responsible for the decrease in recovery from infection. The authors [17] developed and analyzed a population-level compartmental model of human–mosquito interactions that takes into account an intervention using transmission-blocking drugs. Their results show that the effect of treatment rate on reducing reproduction number depends on the efficacy of the drug. A mathematical model was developed by [19] to assess the effectiveness of antimalarial herbal therapy in controlling malaria transmission. Their study suggests a massive campaign for the use of antimalarial traditional medicine as a treatment for malaria infection. However, further empirical studies and clinical trials are needed to validate the model's predictions and assess the feasibility and effectiveness of implementing herbal therapy on a larger scale. The authors in [5] proposed a deterministic model that explains the transmission dynamics between mosquito and human populations and the impacts of control interventions. A malaria dynamic model with partial immunity and protected travelers using optimal control and cost-effectiveness analysis was developed [16]. They showed that a quadruple of control strategies could minimize infection and implementation costs, suggesting not all strategies are necessary with limited resources. In addition, the study in [20] developed a deterministic model to assess the impact of awareness campaigns on malaria transmission. They applied optimal control theory to minimize disease control costs, concluding that media-driven awareness campaigns are effective for cost-efficient malaria management. The authors in [21] formulated a model incorporating drug resistance, treatment, and mosquito net usage to study malaria dynamics in Nigeria. By fitting the model to real-world data, they identified key parameters influencing disease transmission and emphasized the need for improved control measures targeting resistant strains. Again, the study in [22] proposed a malaria model incorporating insecticide-treated bed nets, human treatment, sterile mosquito techniques, and interventions for pregnant women and newborns, demonstrating that effective control strategies can reduce malaria transmission. The study presented in [23] a nonlinear deterministic model to analyze climatic variability's impact on malaria transmission, highlighting the significance of environmental factors in disease control strategies. The researchers developed a model incorporating partial immunity and protected travelers' movement, using optimal control theory to evaluate intervention strategies for malaria management in populations with varying immunity levels and mobility patterns [24]. Moreover, the study introduced a two-group model distinguishing between symptomatic and asymptomatic carriers of malaria [25]. Optimal control strategies were analyzed to determine effective interventions, emphasizing the role of asymptomatic individuals in disease persistence and control. The study in [15] develops a nonlinear mathematical model to investigate the effects of relapse on malaria transmission dynamics in both human and mosquito populations. Their result showed that relapse occurs not only because of incomplete treatment but also when symptoms reappear after the parasites have been eliminated from blood but persist as dormant hypnozoites in liver cells. However, the model does not explicitly incorporate treatment-seeking behavior, which is a critical factor in malaria control efforts.

Despite extensive research, the disease continues to spread due to various factors, including treatment-seeking behaviors. Many individuals seek treatment from traditional medicine practitioners, whose remedies often lack clinical validation. This study aims to fill the gap in mathematical modeling by incorporating the effect of treatment-seeking behaviors on malaria transmission dynamics using optimal control approaches. Thus, this study is an extension of the model studied in [15], taking into account the treatment-seeking behavior. It extends the existing model by classifying infected individuals based on their treatment choices: those seeking treatment at healthcare facilities and those using traditional remedies without clinical validation. The model further integrates three optimal control strategies: antimalarial drug treatment, the use of insecticide-treated bed nets (ITNs), and awareness campaigns to discourage inappropriate reliance on traditional medicine.

In general, the rest of this article is organized as follows: In Section 2 of the work, the mathematical model is formulated. Section 3 discusses the basic properties of the model, and a sensitivity analysis of the basic reproduction number is conducted in Section 4, followed by numerical simulations in Section 5. In Section 6, optimal control analysis is made analytically, and numerical simulations of optimal control are discussed. Then, discussions and conclusions are provided in Section 7. Lastly, in Section 8, limitations of the study are provided.

2 Model formulation

This study is an extension of the model provided in [15], taking into account the treatment-seeking behavior by adding more compartments and assumptions. The total human population at a time t, denoted by N_h(t), is divided into six epidemiological groups in the presence of the disease as susceptible (S_h), infected (E_h), infectious (I_h), infectious individuals who are undergoing treatment at health facilities (T_h), infectious individuals who are undergoing treatment with inappropriate traditional medicine by traditional practitioners (T_s), and recovered (individuals who recover as a result of receiving treatment using traditional medicine or natural immunity R_h). Thus, at any time t ≥ 0, the total human population is given by:

$\begin{array}{l}{N}_h\left(t\right)=S_h\left(t\right)+E_h+I_h\left(t\right)+T_h\left(t\right)+T_s\left(t\right)+R_h\left(t\right).\end{array}$

Likewise, the vector population is categorized into three groups: susceptible (S_m), exposed (E_m), and infectious (I_m). The vector compartment does not include an immune class as mosquitoes never recover from the infection, that is, their infected period ends with their death due to their relatively short lifecycle. Thus, the total vector population denoted N_m(t), at any time t ≥ 0 is given by:

$\begin{array}{l}{N}_m\left(t\right)=S_m\left(t\right)+E_m\left(t\right)+I_m\left(t\right).\end{array}$

Individuals naturally enter the susceptible class S_h at a constant rate Λ_h, due to the natural loss of immunity from the human treatment class at health facilities (a constant rate ρ), by loss of immunity from the recovered class (at a constant rate ρ₁) and can be decreased by a natural death at a constant rate μ or infected after a bite from an infectious mosquito, and the sporozoites are passed on to them. The transmission rate of infections in a susceptible human population (from an infectious mosquito to a susceptible human) is assumed to be given by a rate λ_h. An exposed individual becomes infectious at a constant rate θ_h and an infectious human I_h who is seeking treatment at health facilities will move to the treatment class, T_h at a constant rate τ₁. It is assumed that individuals who are undergoing treatment at health facilities will recover successfully and then go back to the susceptible class S_h by loss of immunity at a constant rate ρ. Again, an infectious human who is seeking treatment with traditional medicine without a prescription from health professionals will move to the treatment class, T_s at a constant rate τ₂. We assumed that an individual in a treatment class, T_s with traditional medicine can recover temporarily due to natural immunity and the use of traditional medicine and will move to the recovered class at a constant rate γ; while the others go back to a health facility due to the ineffectiveness of traditional medicine treatment at a constant progression rate ω. A temporarily recovered human class can be entered into the susceptible human class by losing immunity at a constant rate ρ₁ if the merozoites of parasites clear from the blood completely, and if they don't, they go back to the infectious class I_h at a constant rate γ₁. All human population classes decreased through natural death at a constant rate μ and disease-induced death rates for I_h and T_s at a constant δ and ψ, respectively.

Similarly, new vectors are recruited at a constant rate, Λ_m into the susceptible mosquito compartment, S_m. When a susceptible mosquito bites an infectious human, I_h or human in treatment with traditional medicine, T_s, the parasite (in the form of gametocytes) enters the mosquito, and the mosquito moves from the susceptible class to the exposed class, E_m at a rate Λ_m. An infected mosquito will become infectious at a constant rate θ_m. The vector population classes decreased by natural death at a constant rate η.

The formulation of the model is primarily guided by the following assumptions:

1. When an individual is infected with malaria, we categorize the human treatment class into two groups based on community treatment-seeking behaviors: those receiving treatment at health facilities and those undergoing treatment with traditional medicines without a prescription from healthcare professionals.

2. In this study, the term ‘traditional medicines' (TMs) refers to antimalarial drugs prepared by traditional practitioners. However, it is important to note that these medicines lack clinical validity and are devoid of quality control, safety measures, standardized dosage, and potential drug interactions [26].

3. An infectious individual receiving treatment at health facilities will recover fully and subsequently lose natural immunity, rendering them susceptible to malaria again.

4. An infectious individual undergoing treatment with traditional medicine may experience temporary recovery due to their natural immunity. However, it is assumed that the merozoites of parasites are not completely eliminated from the bloodstream due to the ineffectiveness of traditional antimalarial drugs. Consequently, they may revert to the infectious human class [15].

The flow chart (Figure 1) shows the dynamics of malaria in human and vector populations.

Figure 1

Figure 1. Flow diagram of the malaria transmission dynamics.

The dynamical system of the flow chart of Figure 1 is given as follows:

$\begin{array}{ll}\{\begin{array}{l}\frac{d{S}_h}{dt}\left(t\right)={Ʌ}_h-{\lambda }_h{S}_h+\rho T_h+{\rho }_1{R}_h-\mu S_h\hfill \\ \frac{d{E}_h}{dt}\left(t\right)={\lambda }_h{S}_h-\left({\theta }_h+\mu \right)E_h\hfill \\ \frac{d{I}_h}{dt}\left(t\right)={\theta }_h{E}_h+{\gamma }_1{R}_h-\left(\delta +{\tau }_1+{\tau }_2+\mu \right)I_h\hfill \\ \frac{d{T}_h}{dt}\left(t\right)={\tau }_1{I}_h+\omega T_s-\left(\rho +\mu \right)T_h\hfill \\ \frac{d{T}_s}{dt}\left(t\right)={\tau }_2{I}_h-\left(\omega +\gamma +\mu +\psi \right)T_s\hfill \\ \frac{d{R}_h}{dt}\left(t\right)=\gamma T_s-\left({\gamma }_1+{\rho }_1+\mu \right)R_h\hfill \\ \frac{d{S}_m}{dt}\left(t\right)={Ʌ}_m-\left({\lambda }_m+\eta \right)S_m\hfill \\ \frac{d{E}_m}{dt}\left(t\right)={\lambda }_m{S}_m-\left({\theta }_m+\eta \right)E_m\hfill \\ \frac{d{I}_m}{dt}\left(t\right)={\theta }_m{E}_m-\eta I_m\hfill \end{array}& \left(1\right)\end{array}$

With an initial conditions [S_h(0), E_h(0), I_h(0), T_h(0), T_s(0), R_h(0), S_m(0), E_m(0), I_m(0)]$\in {ℝ}_{+}^9$ and where λ_h and λ_m are the force of infections from vector to human and human to vector, respectively given by

$\begin{array}{ll}{\lambda }_h\left(I_m\right)=\frac{q{\beta }_{mh}{I}_m}{N_h},{\lambda }_m\left(I_h+T_s\right)=\frac{q{\alpha }_{1m}{I}_h}{N_h}+\frac{q{\alpha }_{2m}{T}_s}{N_h}& \left(2\right)\end{array}$

3 Model analysis

The model represented by the system of differential equation (Equation 1) will be analyzed in the feasible region and since the model represents the populations all the state variables and the parameters are assumed positive [26]. The positivity and boundedness are given in lemma 1 and 2.

Lemma 1 (Positivity). [18, 27] If S_h(0) > 0, E_h(0) > 0, I_h(0) > 0, T_h(0) > 0, T_s(0) > 0, R_h(0) > 0, S_m(0) > 0, E_m(0) > 0, I_m(0) > 0, then the solution S_h(t), E_h(t), I_h(t), T_h(t), T_s(t), R_h(t), S_m(t), E_m(t) and I_m(t) of the system (Equation 1) are all positive for all time t ≥ 0.

Proof: From the first equation of system (Equation 1), we have $\frac{d{S}_h}{dt}+\left({\lambda }_h+\mu \right)S_h={\Lambda }_h+\rho T_h+{\rho }_1{R}_h$.

After some simplification and integrating both sides given as:

$\begin{array}{l}{\displaystyle {\int }_0^t}{\left(e^{Q\left(s\right)+\mu s}{S}_h\left(s\right)\right)}^{\prime }ds={\displaystyle {\int }_0^t}\left({\Lambda }_h+\rho T_h+{\rho }_1{R}_h\right)e^{Q\left(s\right)+\mu s}ds.\\ S_h\left(t\right)={S_h\left(0\right)e}^{-Q\left(t\right)+Q\left(0\right)-\mu t}+e^{-Q\left(t\right)-\mu t}\\ {\displaystyle {\int }_0^t}\left({\Lambda }_h+\rho T_h+{\rho }_1{R}_h\right)e^{Q\left(s\right)+\mu \left(s\right)}ds.& \left(3\right)\end{array}$

Since S_h(0) > 0 and f(t) > 0, for all t > 0, and also the exponential function is always positive, then the solution S_h(t) > 0 for all t > 0.

Similarly, it can be shown that the rest of the state variables, E_h(t), I_h(t), T_h(t), T_s(t), R_h(t), S_m(t), E_m(t),I_m(t) are non-negative for all time t > 0.

Thus, we can conclude that all solutions of model (Equation 1) remain positive for all t > 0 and the solution of the model is positive.

Lemma 2 (Boundedness) [28, 29]. The feasible region Ω defined by

$\begin{array}{l}{\Omega }_h=\left\{\right(S_h,E_h,I_h,T_h,T_s,R_h)\in {ℝ}_{+}^6:\\ (S_h+E_h+I_h+T_h+T_s+R_h)\le \frac{{Ʌ}_h}{\mu }\},\\ {\Omega }_m=\left\{\right(S_m,E_m,I_m)\in {ℝ}_{+}^3:(S_m+E_m+I_m)\\ \le {Ʌ}_m/\eta \},\end{array}$

With initial condition S_h(0) > 0, E_h(0) > 0, I_h(0) > 0, T_h(0) > 0, T_s(0) > 0, R_h(0) > 0, S_m(0) > 0, E_m(0) > 0,I_m(0) > 0 is positively invariant for the system (Equation 1) in ${ℝ}_{+}^9$.

Proof: Given N_h(t) = S_h(t) + E_h(t) + I_h(t) + T_h(t) + T_s(t) + R_h(t) and N_m(t) = S_m(t) + E_m(t) + I_m(t). The rate of change in total human population size N_h(t) = S_h(t) + E_h(t) + I_h(t) + T_h(t) + T_s(t) + R_h(t) can be obtained as:

$\begin{array}{l}\frac{d{N}_h}{dt}\left(t\right)=({Ʌ}_h-\mu S_h-\mu E_h-\left(\delta +\mu \right)I_h-\mu T_h\\ -(\mu +\psi )T_S-\mu R_h)\le {Ʌ}_h-\mu N_h\left(t\right).\end{array}$

After some simplification and integrating both sides, we have:

$\begin{array}{ll}{N}_h\left(t\right)\le {N_h\left(0\right)e}^{-\mu t}+\frac{{\Lambda }_h}{\mu }\left(1-e^{-\mu t}\right)& \left(4\right)\end{array}$

Now, taking the limsup as t → ∞, gives $N_h\left(t\right)\le \underset{t\to \infty }{lim}sup{N_h\left(0\right)e}^{-\mu t}+\frac{{\Lambda }_h}{\mu }\left(1-e^{-\mu t}\right)=\frac{{\Lambda }_h}{\mu }$.

$\begin{array}{ll}{N}_h\left(t\right)\le \frac{{\Lambda }_h}{\mu }& \left(5\right)\end{array}$

Similarly, the boundedness of the total vector population is given as:

$\begin{array}{ll}{N}_m\left(t\right)\le \frac{{\Lambda }_m}{\eta }& \left(6\right)\end{array}$

In particular, $N_h\left(t\right)\le \frac{{\Lambda }_h}{\mu }$ if the total human population at the initial instant of time, $N_h\left(0\right)\le \frac{{\Lambda }_h}{\mu }$ and $N_m\left(t\right)\le \frac{{\Lambda }_m}{\eta }$, if $N_m\left(0\right)\le \frac{{\Lambda }_m}{\eta }$. Therefore, the region is positively invariant for the system (Equation 1). Moreover, if $N_h\left(t\right)>\frac{{\Lambda }_h}{\mu }$ and $N_m\left(t\right)>\frac{{\Lambda }_m}{\eta }$ , then either the solution of model (Equation 1) enters Ω in a finite time, or N_h(t) approaches $\frac{{\Lambda }_h}{\mu }$ and N_m(t) approaches $\frac{{\Lambda }_m}{\eta }$ asymptotically. Hence, the region Ω is positively invariant and attracts all solutions of model (Equation 1) in ${ℝ}_{+}^6$x${ℝ}_{+}^3$.

Lemma 3 (Existence and uniqueness). In model (Equation 1) if the initial conditions S_h(0) > 0, E_h(0) > 0, I_h(0) > 0, T_h(0) > 0, T_s(0) > 0, R_h(0) > 0, S_m(0) > 0, E_m(0) > 0,I_m(0) > 0, and t₀ > 0, then for all t∈ℝ the solutions S_h(t), E_h(t), I_h(t), T_h(t), T_s(t), R_h(t), S_m(t), E_m(t) and I_m(t) exist in ${ℝ}_{+}^9$ [17, 28, 30, 31].

Proof: The model (Equation 1) can be described as:

$\begin{array}{ll}Ẋ\left(t\right)=f\left(X\left(t\right)\right),X\left(0\right)=X_0,& \left(7\right)\end{array}$

where X(t)=(S_h(t), E_h(t), I_h(t), T_h(t), T_s(t), R_h(t), S_m(t), E_m(t), $I_m\left(t\right){)}^T$ with the corresponding initial condition.

$\begin{array}{l}{V}^{-1}=\left(\begin{array}{ccccc}\frac{1}{\mu +{\theta }_h}& 0& 0& 0& 0\\ -\frac{-{\eta }^2{\theta }_h-\eta {\theta }_h{\theta }_m}{\eta \left(\mu +{\theta }_h\right)\left(\eta +{\theta }_m\right)\left(\delta +\mu +{\tau }_1+{\tau }_2\right)}& \frac{1}{\delta +\mu +{\tau }_1+{\tau }_2}& 0& 0& 0\\ 0& 0& \frac{1}{\eta +{\theta }_m}& 0& 0\\ 0& 0& \frac{{\theta }_m}{\eta \left(\eta +{\theta }_m\right)}& \frac{1}{\eta }& 0\\ \frac{{\eta }^2{\theta }_h{\tau }_2+\eta {\theta }_h{\theta }_m{\tau }_2}{\eta \left(\gamma +\mu +\psi +\omega \right)\left(\mu +{\theta }_h\right)\left(\eta +{\theta }_m\right)\left(\delta +\mu +{\tau }_1+{\tau }_2\right)}& \frac{{\tau }_2}{\left(\gamma +\mu +\psi +\omega \right)\left(\delta +\mu +{\tau }_1+{\tau }_2\right)}& 0& 0& \frac{1}{\gamma +\mu +\psi +\omega }\\ \end{array}\right)\end{array}$

$\begin{array}{ll}F.V^{-1}=\left(\begin{array}{ccccc}0& 0& \frac{q{\beta }_{hm}{\theta }_m}{\eta \left(\eta +{\theta }_m\right)}& \frac{q{\beta }_{hm}}{\eta }& 0\\ 0& 0& 0& 0& 0\\ -\frac{q\mu {\alpha }_m\left(-{\eta }^2{\theta }_h-\eta {\theta }_h{\theta }_m\right){\Lambda }_m}{{\eta }^2\left(\mu +{\theta }_h\right)\left(\eta +{\theta }_m\right){\Lambda }_h\left(\delta +\mu +{\tau }_1+{\tau }_2\right)}+\frac{q\mu {\alpha }_{2m}{\Lambda }_m\left({\eta }^2{\theta }_h{\tau }_2+\eta {\theta }_h{\theta }_m{\tau }_2\right)}{{\eta }^2\left(\gamma +\mu +\psi +\omega \right)\left(\mu +{\theta }_h\right)\left(\eta +{\theta }_m\right){\Lambda }_h\left(\delta +\mu +{\tau }_1+{\tau }_2\right)}& \frac{q\mu {\alpha }_m{\Lambda }_m}{\eta {\Lambda }_h\left(\delta +\mu +{\tau }_1+{\tau }_2\right)}+\frac{q\mu {\alpha }_{2m}{\Lambda }_m{\tau }_2}{\eta \left(\gamma +\mu +\psi +\omega \right){\Lambda }_h\left(\delta +\mu +{\tau }_1+{\tau }_2\right)}& 0& 0& \frac{q\mu {\alpha }_{2m}{\Lambda }_m}{\eta \left(\gamma +\mu +\psi +\omega \right){\Lambda }_h}\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ \end{array}\right)& \left(8\right)\end{array}$

X₀=(S_h(0), E_h(0), I_h(0), T_h(0), T_s(0), R_h(0), S_m(0), E_m(0), I_m(0)) and f(X) is the vector-valued function representing the hand-right side of the system (Equation 1) which is

$\begin{array}{ll}f\left(X\right)=\left(\begin{array}{c}\frac{{\Lambda }_h-{{\lambda }_{h}S}_h+\rho T_h+{\rho }_1{R}_h-\mu S_h}{{{\lambda }_{h}S}_h-\left({\theta }_h+\mu \right)E_h}\\ {\theta }_h{E}_h+{\gamma }_1{R}_h-\left(\delta +{\tau }_1+{\tau }_2+\mu \right)I_h\\ {\tau }_1{I}_h+\omega T_s-\left(\rho +\mu \right)T_h\\ {\tau }_2{I}_h-\left(\omega +\gamma +\mu +\psi \right)T_s\\ \gamma T_s-\left({\gamma }_1+{\rho }_1+\mu \right)R_h\\ {\Lambda }_m-\left({\lambda }_m+\eta \right)S_m\\ {\lambda }_m{S}_m-\left({\theta }_m+\eta \right)E_m\\ {\theta }_m{E}_m-\eta I_m\\ \end{array}\right)& \left(9\right)\end{array}$

Since a function f in the system (Equation 1) has a continuous first derivative in ${ℝ}_{+}^9$ and $X_0\in {ℝ}_{+}^9$, N_h≠0 and N_m≠0, f is well defined at t = 0. Consequently, according to the widely used fundamental existence and uniqueness theorem [32], lemmas 1 and 2 above, there exists a unique, positive, and bounded solution for the model (Equation 1) in ${ℝ}_{+}^9$.

3.1 Malaria-free equilibrium point, E₀

Malaria-free equilibrium points (MFE) are steady-state solutions where there is no malaria [33]. The malaria-free equilibrium of the system (Equation 1) is $E_0=\left(\frac{{\Lambda }_h}{\mu },0,0,0,0,0,\frac{{\Lambda }_m}{\eta },0,0\right)$.

3.2 The basic reproduction number

The reproduction number $\left({ℛ}_0\right)$ is defined as the average number of secondary infections produced by one infected individual in a completely susceptible population [13]. To obtain ${\mathcal{R}}_0$ for model (Equation 1), we use the next-generation matrix technique described in [18, 34] and is the spectral radius ρ(FV⁻¹), where

$\begin{array}{l}F=\left(\begin{array}{ccccc}0& 0& 0& q{\beta }_{mh}& 0\\ 0& 0& 0& 0& 0\\ 0& \frac{q\mu {\alpha }_{1m}{\Lambda }_m}{\eta {\Lambda }_h}& 0& 0& \frac{q\mu {\alpha }_{2m}{\Lambda }_m}{\eta {\Lambda }_h}\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ \end{array}\right)\text{and}\\ V=\left(\begin{array}{ccccc}\left({\theta }_h+\mu \right)& 0& 0& 0& 0\\ -{\theta }_h& \left(\delta +{\tau }_1+{\tau }_2+\mu \right)& 0& 0& 0\\ 0& 0& \left({\theta }_m+\eta \right)& 0& 0\\ 0& 0& -{\theta }_m& \eta & 0\\ 0& -{\tau }_2& 0& 0& \left(\omega +\gamma +\mu +\psi \right)\\ \end{array}\right)\end{array}$

By next generation operator method, the basic reproduction number of the model (Equation 1) is given:

$\begin{array}{l}\text{ }{ℛ}_0=\rho \left(F{V}^{-1}\right)\\ =\sqrt{\frac{q^2\mu {\beta }_{hm}{\theta }_h{\theta }_m{Ʌ}_m\left(\left(\gamma +\mu +\omega +\psi \right){\alpha }_{1m}+{\tau }_2{\alpha }_{2m}\right)}{{\eta }^2{Ʌ}_h\left(\mu +{\theta }_h\right)\left(\eta +{\theta }_m\right)\left(\gamma +\mu +\omega +\psi \right)\left(\delta +\mu +{\tau }_1+{\tau }_2\right)}}& \left(10\right)\end{array}$

Theorem 1. [13, 35, 36] The malaria-free equilibrium point, E₀ of the model in the system (Equation 1) is locally asymptotically stable if the reproduction number ${\mathcal{R}}_0<1$ and is unstable if${\mathcal{R}}_0>1$.

Proof: The Jacobian matrix J(E₀ ) of the model to (S_h, E_h, I_h, T_h, T_S, R_h, S_m, E_m, I_m) at the malaria-free equilibrium point is the following:

$\begin{array}{ll}J\left(E_0\right)=\left(\begin{array}{ccccccccc}\mu & 0& 0& \rho & 0& {\rho }_1& 0& 0& -A\\ 0& -B& 0& 0& 0& 0& 0& 0& A\\ 0& {\theta }_h& -C& 0& 0& {\gamma }_1& 0& 0& 0\\ 0& 0& {\tau }_1& -D& \omega & 0& 0& 0& 0\\ 0& 0& {\tau }_2& 0& -E& 0& 0& 0& 0\\ 0& 0& 0& 0& \gamma & -F& 0& 0& 0\\ 0& 0& -G& 0& -H& 0& -\eta & 0& 0\\ 0& 0& G& 0& H& 0& 0& -I& 0\\ 0& 0& 0& 0& 0& 0& 0& {\theta }_m& -\eta \\ \end{array}\right)& \left(11\right)\end{array}$

where A = qβ_mh, B = (θ_h+μ), C = (δ+τ₁+τ₂+μ), D = (ρ+μ), E = (ω+γ+μ+ψ), F = (γ₁+ρ₁+μ), $G=\frac{q\mu {\alpha }_{1m}{\Lambda }_m}{\eta {\Lambda }_h}$, $H=\frac{q\mu {\alpha }_{2m}{\Lambda }_m}{\eta {\Lambda }_h}$, I = (θ_m+η), where eigenvalues of λ₁ = −μ and λ₆ = − η.

The remaining eigenvalues are the roots of the characteristic equation for Equation 11 given by

$\begin{array}{l}\left|\begin{array}{ccccccc}-B-\lambda & 0& 0& 0& 0& 0& A\\ {\theta }_h& -C-\lambda & 0& 0& {\gamma }_1& 0& 0\\ 0& {\tau }_1& -D-\lambda & \omega & 0& 0& 0\\ 0& {\tau }_2& 0& -E-\lambda & 0& 0& 0\\ 0& 0& 0& \gamma & -F-\lambda & 0& 0\\ 0& G& 0& H& 0& -I-\lambda & 0\\ 0& 0& 0& 0& 0& {\theta }_m& -\eta -\lambda \end{array}\right|\\ =0& \left(12\right)\end{array}$

$\begin{array}{l}\left(-D-\lambda \right)(-A\left(F+\lambda \right){\theta }_h{\theta }_m\left(G\left(E+\lambda \right)+H{\tau }_2\right)\\ +\left(I+\lambda \right)\left(B+\lambda \right)\left(\eta +\lambda \right)(\left(C+\lambda \right)\left(E+\lambda \right)\left(F+\lambda \right)\\ -\gamma {\gamma }_1{\tau }_2))=0& \left(13\right)\end{array}$

Now, we have λ₃ = −(ρ+μ) and the remaining eigenvalues are obtained (Equation 14):

$\begin{array}{ll}{P_6\lambda }^6+{P_5\lambda }^5+{P_4\lambda }^4+{P_3\lambda }^3+{P_2\lambda }^2+P_1\lambda +P_0=0& \left(14\right)\end{array}$

Where, P₆ = 1 > 0, P₅ = (γ+δ+2η+4μ+ω+γ₁+θ_h+θ_m+ρ₁+τ₁+τ₂) > 0, P₄ > 0, P₃ > 0, P₂ > 0, $P_1=(IBCEF+IBCF\eta +IBEF\eta +ICEF\eta +BCEF\eta +IBCE\eta \left(1-{ℛ}_0^2\right)$ −AFGθ_hθ_m−(IB+Iη+Bη)γγ₁τ₂) and $P_0=IB\eta \left(CEF\left(1-{{\mathcal{R}}_0}^2\right)-\gamma {\gamma }_1{\tau }_2\right)>0$, where ${\mathcal{R}}_0<1$ and γγ₁τ₂<CEF.

Applying the Routh–Hurwitz stability criterion [37] and after some little algebraic manipulations, it can be shown that the eigenvalues of the block matrix have negative real parts .i.e. λ₁ = −μ < 0 , λ₃ = −(ρ+μ) < 0, λ₆ = −η < 0, if ${\mathcal{R}}_0<1$. If ${\mathcal{R}}_0>1$, then P₁ < 0, thus the matrix JE₀ has at least one eigenvalue with a positive real part. Hence, malaria-free equilibrium E₀ of the system (Equation 1) is locally asymptotically stable if ${\mathcal{R}}_0<1$ and unstable if ${\mathcal{R}}_0>1$.

3.3 Malaria-present equilibrium point, E^*

The endemic equilibrium point of the dynamical system of Equation 4 is obtained by making the right side of the system equal to zero, providing that E_h ≠ 0, I_h ≠ 0, T_s ≠ 0, R_h ≠ 0, E_m ≠ 0, I_m ≠ 0. We have supposed the endemic equilibrium point of the model is denoted by E^* $={S_h}^{*}$, ${E_h}^{*}$, ${I_h}^{*}$, ${T_h}^{*}$, ${T_S}^{*}$, ${R_h}^{*}$, ${S_m}^{*}$, ${E_m}^{*}$, ${I_m}^{*}$ and the corresponding forces of infection are:

${{\lambda }_h}^{*}$ $=\frac{q{\beta }_{mh}{I_m}^{*}{S_h}^{*}}{N_h}$ and ${{\lambda }_m}^{*}=\frac{q{\alpha }_{1m}{I_h}^{*}{S_m}^{*}}{N_h}+\frac{q{\alpha }_{2m}{T_s}^{*}{S_m}^{*}}{N_h}$. Next, after some steps of the computations, we have determined the required endemic equilibrium point for the model given as:

$\begin{array}{l}{S}_h^{\ast }=\frac{\left({\theta }_h+\mu \right)}{{\lambda }_h}(\frac{\left(\delta +{\tau }_1+{\tau }_2+\mu \right)}{{\theta }_h}\\ -\frac{{\gamma }_1\gamma {\tau }_2}{{\theta }_h\left({\gamma }_1+{\rho }_1+\mu \right)\left(\omega +\gamma +\mu +\psi \right)})I_h^{\ast },E_h^{\ast }\\ =\left(\frac{\left(\delta +{\tau }_1+{\tau }_2+\mu \right)}{{\theta }_h}-\frac{{\gamma }_1\gamma {\tau }_2}{{\theta }_h\left({\gamma }_1+{\rho }_1+\mu \right)\left(\omega +\gamma +\mu +\psi \right)} \right)I_h^{\ast }\\ T_h^{\ast }=\left(\frac{{\tau }_1}{\left(rho+\mu \right)}+\frac{\omega {\tau }_2}{\left(\omega +\gamma +\mu +\psi \right)}\right)I_h^{\ast },\\ T_S^{\ast }=\left(\frac{{\tau }_2}{\left(\omega +\gamma +\mu +\psi \right)}\right)I_h^{\ast },\\ R_h^{\ast }=\left(\frac{\gamma {\tau }_2}{\left({\gamma }_1+{\rho }_1+\mu \right)\left(\omega +\gamma +\mu +\psi \right)} \right)I_h^{\ast },\\ I_h^{\ast }=-\frac{{Ʌ}_h}{\frac{\rho {\tau }_1}{\mu +\rho }+\frac{\rho \omega {\tau }_2}{\gamma +\mu +\psi +\omega }+\frac{\gamma {\rho }_1{\tau }_2}{(\gamma +\mu +\psi +\omega )(\mu +{\gamma }_1+{\rho }_1)}+\frac{(-\mu -{\theta }_h)(\mu +{\lambda }_h)(\delta +\mu +{\tau }_1+(1-\frac{\gamma {\gamma }_1}{(\gamma +\mu +\psi +\omega )(\mu +{\gamma }_1+{\rho }_1)}){\tau }_2)}{{\theta }_h{Ʌ}_h}},\\ S_m^{\ast }=\left(\frac{{Ʌ}_m}{\left({\lambda }_m+\eta \right)}\right),E_m^{\ast }=\left(\frac{{Ʌ}_m{\lambda }_m}{\left({\lambda }_m+\eta \right)\left({\theta }_m+\eta \right)}\right)and{I}_m^{\ast }=\left(\frac{{\theta }_m{\lambda }_m{Ʌ}_m}{\left({\theta }_m+\eta \right)\left({\lambda }_m+\eta \right)} \right).\end{array}$

Then, by substituting ${I_h}^{*}$ and ${I_m}^{*}$ in the ${{\lambda }_h}^{*}$ and ${{\lambda }_m}^{*}$, respectively, we get the simplified form:

$\begin{array}{ll}{P}_1{\lambda }_h+P_0=0& \left(15\right)\end{array}$

where P₁ = (IBDEΛ_h(FC−γγ₁τ₂)+Aθ_mΛ_m(BDEF(δ+μ) +EF(BD−ρθ_h)τ₁+D(BEF−BEγγ₁−θ_h(Fρω+γρ₁))τ₂)) and $P_0=D{\eta }^2{\Lambda }_{h}IB\left(CEF\left({{\mathcal{R}}_0}^2-1\right)+\gamma {\gamma }_1{\tau }_2\right)$.

Consequently, the following result has been determined.

Theorem 2. The model of the system (Equation 1) admits precisely [17].

a. One endemic equilibrium point if ${\mathcal{R}}_0>1$ and P₁ < 0 or P₀ > 0 and P₁ < 0.

b. One endemic equilibrium point if ${\mathcal{R}}_0<1$, γγ₁τ₂ > CEF and P₁ < 0 or P₀ > 0 and P₁ < 0.

c. One endemic equilibrium point if ${\mathcal{R}}_0<1$, γγ₁τ₂<CEF and P₁ > 0 or P₀ < 0 and P₁ > 0.

d. No equilibrium point otherwise.

From theorem 2, we can observe that, the existence of one positive root in condition (b) and one positive root in condition (c) suggests that the model can exhibit backward bifurcation phenomenon, where both locally asymptotically stable malaria-free and endemic steady states co-exist, where ${\mathcal{R}}_0<1$. It suggests that while eliminating diseases from the population is dependent on the criterion ${\mathcal{R}}_0<1$, this need is no longer sufficient to ensure disease eradication. Thus, when ${\mathcal{R}}_0<{{\mathcal{R}}_0}^{*}<1$, it can be demonstrated that backward bifurcation happens for certain values of ${\mathcal{R}}_0$. The backward bifurcation phenomena in system (Equation 1) make it challenging to manage or eradicate malaria because ${\mathcal{R}}_0<{{\mathcal{R}}_0}^{*}<1$ requires lowering the model's reproduction number to a level below unity. This implies that further control efforts will be required to eradicate malaria from the population.

Theorem 3. The malaria-present equilibrium, $E^{*}=\left({S_h}^{*},{E_h}^{*},{I_h}^{*},{T_h}^{*},{T_S}^{*},{R_h}^{*},{S_m}^{*},{E_m}^{*},{I_m}^{*}\right)$ of the system (Equation 1), is locally asymptotically stable (LAS), if ${\mathcal{R}}_0>1$ and unstable if ${\mathcal{R}}_0<1$ [38].

Proof: The Jacobian matrix of the model (Equation 1) is given by

$\begin{array}{ll}J\left(E^{*}\right)=\left(\begin{array}{ccccccccc}{a}_1& 0& 0& a_4& 0& a_6& 0& 0& -a_9\\ b_1& -b_2& 0& 0& 0& 0& 0& 0& b_9\\ 0& c_2& -c_3& 0& 0& c_6& 0& 0& 0\\ 0& 0& d_3& -d_4& d_5& 0& 0& 0& 0\\ 0& 0& e_3& 0& -e_5& 0& 0& 0& 0\\ 0& 0& 0& 0& f_5& -f_6& 0& 0& 0\\ 0& 0& -g_3& 0& -g_5& 0& -g_7& 0& 0\\ 0& 0& h_3& 0& h_5& 0& h_7& -h_8& 0\\ 0& 0& 0& 0& 0& 0& 0& k_8& -k_9\\ \end{array}\right)& \left(16\right)\end{array}$

where $a_1=\mu +\frac{q{\beta }_{mh}{I_m}^{*}}{{N_h}^{*}}$, a₄ = ρ, a₆ = ρ₁, $a_9=\frac{q{\beta }_{mh}{S_h}^{*}}{{N_h}^{*}}$, $b_1=\frac{q{\beta }_{mh}{I_m}^{*}}{{N_h}^{*}}$, b₂ = θ_h+μ, $b_9=\frac{q{\beta }_{mh}{S_h}^{*}}{{N_h}^{*}}$ , c₂ = θ_h, c₃ = (δ+τ₁+τ₂+μ), c₆ = γ₁, d₃ = τ₁, d₄ = (ρ+μ), d₅ = ω, e₁ = e₂ = e₃ = τ₂, e₅ = (ω+γ+μ+ψ), f₅ = γ, f₆ = γ₁+ρ₁+μ, $g_3=\frac{q{\alpha }_{1m}}{{N_h}^{*}}{S_m}^{*}$, $g_5=\frac{q{\alpha }_{2m}}{{N_h}^{*}}{S_m}^{*}$, $g_7=\left(q{\alpha }_{1m}\frac{{I_h}^{*}}{{N_h}^{*}}+q{\alpha }_{2m}\frac{{T_S}^{*}}{{N_h}^{*}}+\eta \right)$, $h_3=\frac{q{\alpha }_{1m}}{{N_h}^{*}}{S_m}^{*}$, $h_5=\frac{q{\alpha }_{2m}}{{N_h}^{*}}{S_m}^{*}$, $h_7=\left(q{\alpha }_{1m}\frac{{I_h}^{*}}{{N_h}^{*}}+q{\alpha }_{2m}\frac{{T_S}^{*}}{{N_h}^{*}}\right)$, h₈ = θ_m+η, k₈ = θ_m, k₉ = η. The corresponding characteristic equation of the Jacobian matrix with eigenvalue λ is given by |J(E^*)−λI| = 0; that is,

$\begin{array}{l}{P_9\lambda }^9+{P_8\lambda }^8+{P_7\lambda }^7+{P_6\lambda }^6+{P_5\lambda }^5+{P_4\lambda }^4+{p_3\lambda }^3\\ +{P_2\lambda }^2+P_1\lambda +P_0=0& \left(17\right)\end{array}$

where P₀, P₁, P₂, …, P₉ (see Appendix I).

Therefore, we show that when ${\mathcal{R}}_0>1$, all the coefficients P_i of the characteristic (Equation 17), and the first column values b₁, c₁, d₁, e₁, f₁, g₁, h₁ and i₀ of the Routh array are positive, so by the Routh-Hurwitz stability criterion, all the eigenvalues of the Jacobian matrix (Equation 17) have negative real parts. Thus, the malaria-present equilibrium is locally asymptotically stable for ${\mathcal{R}}_0>1$.

4 Sensitivity analysis of the reproduction number

This section aims to perform a sensitivity analysis of the basic reproductive number, which is crucial for designing mitigation strategies to slow the spread of malaria by reducing the reproduction number. Sensitivity analysis helps [39–41] researchers, public health officials, and policymakers prioritize interventions and decisions based on the most strongly influencing factors, enabling them to understand the effects of each parameter on the reproduction number and design effective strategies.

Definition 1. Normalized forward sensitivity index of ${\mathcal{R}}_0$ which is differentiable to a given parameter (p) is defined as [42]

$\begin{array}{ll}SI\left(p\right)=\left(\frac{\partial {\mathcal{R}}_0}{\partial p}\right)*\left(\frac{p}{{\mathcal{R}}_0}\right)& \left(18\right)\end{array}$

By using the parameter values given in Table 1 and Def. 1, the respective sensitivity indices values for reproduction number are computed in Table 2 and we plot the sensitivity indices in Figure 2.

Table 1

Table 1. Parameters description of the model with their values.

Table 2

Table 2. Sensitivity indices of ${\mathcal{R}}_0$ to parameters evaluated at the parameter values given in Table 1.

Figure 2

Figure 2. Sensitivity indices of the model parameters for reproduction number, ${\mathcal{R}}_0.$

The positive parameter values increase the reproduction number when they are increased, while the negative parameter values decrease the reproduction number when they are increased. From the parameter values displayed in Table 2, we can observe that factors such as the biting rate of mosquitoes and the probability of transmission of infection from humans in treatment with traditional medicine to susceptible mosquitoes are the most sensitive factors. So, interventions should be introduced to reduce these factors. In another way, factors such as the rate of infectious humans who are undergoing treatment at health facilities (τ₁) and the progression rate (ω) due to failure of treatment with traditional medicine to treatment at health facilities are also sensitive to the basic reproduction number. Increasing these factors will decrease the reproductive number and reduce malaria. Moreover, if SI(p)=y, then a 1% increase in the parameter p results in a y% increase if y > 0 (decrease if y < 0) in reproduction number. From our sensitivity indices results, for example, a 10% increase in the probability of transmission of infection from human treatment with traditional medicine to susceptible mosquitoes results in a 3.5% increase in the value of reproduction number. A 10% increase in treatment with traditional medicine results in a 1.9% decrease in value of the reproduction number; a 10% increase in the progression rate (ω) results in a 3% decrease in the value of the reproduction number. Again, an 10% increase in treatment at health facilities results in a 3.1% decrease in the value of the reproduction number and, hence a reduction in malaria.

5 Numerical simulations

This section conducts a numerical simulation to validate the analytical results of a model (Equation 1) designed to study the dynamics of malaria in communities. Parameter values from many sources are used in the model, which employs the fourth-order Runge–Kutta method and MATLAB R2023a to examine the impact of parameters on malaria transmission.

Here, one of our main goals is to look into how treatment with inappropriate traditional medicines affects the dynamics of malaria transmission. Figure 3a shows how varying the treatment rate at health facilities affects the infectious human populations. We can observe that when the treatment rate at health facilities increases, the population of infectious humans can be reduced, and the disease can be eliminated from the community after 140 months.

Figure 3

Figure 3. Effect of treatment-seeking behaviors on infectious humans. (a) The effect of treatment at health facilities on infectious humans. (b) The effect of treatment with traditional medicines on infectious humans.

Figure 3b shows how varying the treatment rate with traditional medicines affects the infectious human populations. We can observe that when the treatment rate with traditional medicine increases, the population of infectious humans can be increased because of factors such as lack of quality control, safety, dosage, and potential drug interactions and can be reduced when traditional remedies are considered with caution and in consultation with healthcare professionals or skilled traditional practitioners. If so, it takes more than 400 months to eliminate the disease in the targeted community.

Figure 4a shows how varying the progression rate (awareness) to health facilities affects the human population's treatment with traditional medicines. We can observe that when the progression rate from T_s to T_h increases, the population treated with traditional medicines decreases, resulting in the population of infectious humans being reduced, and the disease can be eliminated from the community after 200 months, as shown in Figure 4b.

Figure 4

Figure 4. Effect of progression rate on T_s and I_h. (a) Effect of progression rate from treatment with TM to health facility on T_s. (b) Effect of progression rate from treatment with TM to health facility on I_h.

The depicted graph in Figure 5 delineates the relationship between the reproduction number and the treatment rate (τ₁) for the transition from the infectious human class to the human treatment class at a health facility. Remarkably, a notable intersection occurs at a reproduction number of one and a treatment rate of τ₁ = 0.17. This intersection point carries significant epidemiological implications, suggesting a critical threshold where the spread of infection is balanced by treatment efforts. Now, the rate of new infections matches the rate of individuals transitioning to treatment, potentially indicating a pivotal point for disease control and management strategies. Understanding and leveraging this intersection point can inform targeted interventions aimed at effectively curtailing outbreaks and minimizing their impact on public health.

Figure 5

Figure 5. Effect of treatment rate at health facility on reproduction number, ${\mathcal{R}}_0$.

The convergence of the treatment rate (τ₂) and the reproduction number on the Figure 6 denotes a significant intersection point, occurring at a value of 0.80. This intersection holds profound implications for the transition from the infectious human class to the human treatment class with traditional medicines. At this point, the treatment rate aligns with the reproduction number, indicating a pivotal balance between the rate of new infections and the effectiveness of treatment interventions. Notably, to achieve a reproduction number of less than one, it is essential for the treatment rate (τ₂) to exceed 0.80. Such equilibrium underscores a critical threshold in disease management, where treatment efforts are optimally poised to mitigate the spread of infection effectively. Recognizing and leveraging this intersection point is paramount for devising targeted strategies to control outbreaks and safeguard public health.

Figure 6

Figure 6. Effect of treatment rate with traditional medicines on reproduction number, ${\mathcal{R}}_0$.

6 Formulation and analysis of optimal control

In this part, we use Pontryagin's maximum principle [45, 46] to identify the parameters required for effective malaria control. We use three time-dependent controls in the model (Equation 1) to find the best disease control approach. The study proposes a model to minimize malaria transmission by introducing control measures like ITNs, antimalarial drug treatment, and raising awareness about avoiding inappropriate traditional medicines. The model is designed to investigate the potential impact of these interventions, with the use of insecticide-treated bed nets, antimalarial drug treatment, and awareness creation being incorporated into the model to determine the most effective approach as follows:

$\begin{array}{ll}\{\begin{array}{c}\frac{d{S}_h}{dt}\left(t\right)={Ʌ}_h-\left(\right(1-u_1){\lambda }_h+\mu )S_h+\rho T_h+{\rho }_1{R}_h\\ \frac{d{E}_h}{dt}\left(t\right)=(1-u_1){\lambda }_h{S}_h-\left({\theta }_h+\mu \right)E_h\\ \frac{d{I}_h}{dt}\left(t\right)={\theta }_h{E}_h+{\gamma }_1{R}_h-\left(\delta +u_2{\tau }_1+(1-u_3){\tau }_2+\mu \right)I_h\\ \frac{d{T}_h}{dt}\left(t\right)=u_2{\tau }_1{I}_h+\omega T_s-\left(\rho +\mu \right)T_h\\ \frac{d{T}_s}{dt}\left(t\right)=\left(1-u_3\right){\tau }_2{I}_h-\left(\omega +\gamma +\mu +\psi \right)T_s\\ \frac{d{R}_h}{dt}\left(t\right)=\gamma T_s-\left({\gamma }_1+{\rho }_1+\mu \right)R_h\\ \frac{d{S}_m}{dt}\left(t\right)={Ʌ}_m-\left((1-u_1){\lambda }_m+\eta \right)S_m\\ \frac{d{E}_m}{dt}\left(t\right)=(1-u_1){\lambda }_m{S}_m-\left({\theta }_m+\eta \right)E_m\\ \frac{d{I}_m}{dt}\left(t\right)={\theta }_m{E}_m-\eta I_m\end{array}& \left(19\right)\end{array}$

with the corresponding initial conditions

$\begin{array}{l}{S}_h\left(0\right)>0,E_h\left(0\right)\ge 0,I_h\left(0\right)\ge 0,T_h\left(0\right)\ge 0,T_s\left(0\right)\ge 0,\\ R_h\left(0\right)\ge 0,S_m\left(0\right)>0,E_m\left(0\right)\ge 0,I_m\left(0\right)\ge 0& \left(20\right)\end{array}$

and 0 ≤ u₁(t) ≤ 1 represents malaria infection personal protection measures, 0 ≤ u₂ ≤ 1 represents antimalarial drug treatment control and 0 ≤ u₃ ≤ 1 represents providing awareness on avoiding the inappropriate use of traditional medicines.

The objective function of the system (Equation 19) is to minimize the total number of exposed humans, infectious humans, and mosquitoes while minimizing the cost of control u₁(t), u₂(t), u₃(t) and to find the optimal control values $u^{*}=\left({u_1}^{*},{u_2}^{*},{u_3}^{*}\right)$ of the controls u = (u₁, u₂, u₃) such that the associated state trajectories $\left({S_h}^{*},{E_h}^{*},{I_h}^{*},{T_h}^{*},{T_s}^{*},{R_h}^{*},{S_m}^{*},{E_m}^{*},{I_m}^{*}\right)$ are solutions of the optimal control system (Equation 19) with initial conditions as given in Equation 20. The objective function is given by

$\begin{array}{l}J\left(u_1,u_2,u_3\right)={\displaystyle {\int }_0^{t_f}}{(}_D{E}_h\left(t\right)+D_2{I}_h\left(t\right)+D_3{N}_m\left(t\right)\\ +\frac{1}{2}\left(w_1{u}_1^2+w_2{u}_2^2+w_3{u}_3^2\right))dt& \left(21\right)\end{array}$

Subject to (Equation 22) where t_f is the final time, N_m = (S_m+E_m+I_m), D₁, D₂, D₃, w₁, w₂ and w₃ are positive weight constants. Our choice of control functions agrees with other literature on control of epidemics [16, 40, 46, 47]. The quadratic costs $w_1{u}_1^2\left(t\right)$, $w_2{u}_2^2\left(t\right)$ and $w_3{u}_3^2$ are the costs associated with the use of personal protection measures (ITNs), the use of antimalarial drugs for infectious humans, and awareness on avoiding the un-prescribed use of antimalarial drugs. Here, the fact that the controls are linearly in Equation 19 and quadratic in the objective function allows the Hamiltonian associated with the optimal control problem to be maximized. We seek to find, using the Pontryagin maximum principle [48], an optimal control ${u_1}^{*}\left(t\right),{u_2}^{*}\left(t\right)$ and ${u_3}^{*}\left(t\right)$ satisfying (Equation 19) such that

$\begin{array}{ll}J\left({u_1}^{*},u_2^{*},u_3^{*}\right)=min\left\{J\left(u_1,u_2,u_3\right):u_1,u_2,u_3\in U\right\},& \left(22\right)\end{array}$

Where U = {(u₁, u₂, u₃)} such that u₁, u₂, u₃ are measurable with 0 ≤ u₁ ≤ 1, 0 ≤ u₂ ≤ 1, 0 ≤ u₃ ≤ 1, for t∈[0, t_f] is the control set. The controls (u₁(t), u₂(t), u₃(t)) are bounded Lebesgue integrable functions [49].

The optimal controls (u₁(t), u₂(t), u₃(t)) satisfy necessary conditions from Pontryagin's maximum principle [48]. This principle converts Equations 19, 21 into a problem of minimizing pointwise a Hamiltonian H, to the controls (u₁(t), u₂(t), u₃(t)). The optimality conditions are obtained using the Hamiltonian formulated from the cost functional (Equation 21) and the governing dynamics (Equation 19). Thus, Hamiltonian $\mathcal{H}$ is given as:

$\begin{array}{l}\mathscr{H}=D_1{E}_h\left(t\right)+D_2{I}_h\left(t\right)+D_3{N}_m\left(t\right)+\frac{1}{2}\left(w_1{u}_1^2+w_2{u}_2^2+w_3{u}_3^2\right)\\ +{\Gamma }_1\left[{Ʌ}_h-\left(\right(1-u_1){\lambda }_h+\mu )S_h+\rho T_h+{\rho }_1{R}_h\right]\\ +{\Gamma }_2\left[\left(1-u_1\right){\lambda }_h{S}_h-\left({\theta }_h+\mu \right)E_h\right]\\ +{\Gamma }_3\left[{\theta }_h{E}_h+{\gamma }_1{R}_h-\left(\delta +u_2{\tau }_1+\left(1-u_3\right){\tau }_2+\mu \right)I_h\right]\end{array}$

$\begin{array}{l}+{\Gamma }_4\left[u_2{\tau }_1{I}_h+\omega T_s-\left(\rho +\mu \right)T_h\right]\\ +{\Gamma }_5\left[\left(1-u_3\right){\tau }_2{I}_h-\left(\omega +\gamma +\mu +\psi \right)T_s\right]\\ +{\Gamma }_6\left[\gamma T_s-\left({\gamma }_1+{\rho }_1+\mu \right)R_h\right]\\ +{\Gamma }_7\left[{\Lambda }_m-\left(\left(1-u_1\right){\lambda }_m+\eta \right)S_m\right]\\ +{\Gamma }_8\left[\left(1-u_1\right){\lambda }_m{S}_m-\left({\theta }_m+\eta \right)E_m\right]\\ +{\Gamma }_9\left[{\theta }_m{E}_m-\eta I_m\right]& \left(23\right)\end{array}$

Where Γ_i, i = 1, 2, …, 9 are the associated adjoints for the state variables.

The system of equations is found by taking the appropriate partial derivatives of the Hamiltonian (26) to the associated state variable.

Theorem 4. [45, 47] Let ${u_1}^{*},u_2^{*},u_3^{*}$ be the optimal control and solutions $\left({S_h}^{*},{E_h}^{*},{I_h}^{*},{T_h}^{*},{T_s}^{*},{R_h}^{*},{S_m}^{*},{E_m}^{*},{I_m}^{*}\right)$ of the corresponding state system (19) that minimizes J(u₁, u₂, u₃) over U. Then, there exists adjoint variables Γ_i, i = 1, 2, …, 9 satisfying the following canonical equations

$\begin{array}{l}-\frac{d{\Gamma }_1}{dt}=-{\Gamma }_1\left(\left(1-u_1\right){\lambda }_h+\mu \right)+{\Gamma }_2\left(1-u_1\right)\\ -\frac{d{\Gamma }_2}{dt}=D_1-{\Gamma }_2\left({\theta }_h+\mu \right)+{\Gamma }_3{\theta }_h\\ -\frac{d{\Gamma }_3}{dt}=D_2-{\Gamma }_3\left(\delta +u_2{\tau }_1+\left(1-u_3\right){\tau }_2+\mu \right)+{\Gamma }_4\left(u_2{\tau }_1\right)\\ +{\Gamma }_5\left(1-u_3\right){\tau }_2+\left({\Gamma }_8-{\Gamma }_7\right)\frac{\left(1-u_1\right)q{\alpha }_{1m}{S}_m}{N_h}\\ -\frac{d{\Gamma }_4}{dt}={\Gamma }_1\rho -{\Gamma }_4\left(\rho +\mu \right)\end{array}$

$\begin{array}{l}-\frac{d{\Gamma }_5}{dt}={\Gamma }_4\omega -{\Gamma }_5\left(\omega +\gamma +\mu +\psi \right)+{\Gamma }_6\gamma \\ +\left({\Gamma }_8-{\Gamma }_7\right)\frac{\left(1-u_1\right)q{\alpha }_{2m}{S}_m}{N_h}& \left(24\right)\end{array}$

$\begin{array}{l}-\frac{d{\Gamma }_6}{dt}={\Gamma }_1{\rho }_1+{\Gamma }_3{\gamma }_1-{\Gamma }_6\left({\gamma }_1+{\rho }_1+\mu \right)\\ -\frac{d{\Gamma }_7}{dt}=D_3+{\Gamma }_8\left(1-u_1\right){\lambda }_m-{\Gamma }_7\left(\left(1-u_1\right){\lambda }_m+\eta \right)\\ -\frac{d{\Gamma }_8}{dt}=D_3+{\Gamma }_9{\theta }_m-{\Gamma }_8\left({\theta }_m+\eta \right)\\ -\frac{d{\Gamma }_9}{dt}=D_3+\left({\Gamma }_2-{\Gamma }_1\right)\left(1-u_1\right)q{\beta }_{mh}-{\Gamma }_9\eta \end{array}$

With transversality conditions

$\begin{array}{l}{\Gamma }_1\left(t_f\right)={\Gamma }_2\left(t_f\right)={\Gamma }_3\left(t_f\right)={\Gamma }_4\left(t_f\right)={\Gamma }_5\left(t_f\right)\\ ={\Gamma }_6\left(t_f\right)={\Gamma }_7\left(t_f\right)={\Gamma }_8\left(t_f\right)={\Gamma }_9\left(t_f\right)=0& \left(25\right)\end{array}$

and

$\begin{array}{l}{u_1}^{*}=max\left\{0,min\left\{\frac{\left({\Gamma }_2-{\Gamma }_1\right){\lambda }_h{S}_h+\left({\Gamma }_8-{\Gamma }_7\right){\lambda }_m{S}_m}{w_1},1\right\}\right\},\\ u_2^{*}=max\left\{0,min\left\{\frac{\left({\Gamma }_3-{\Gamma }_4\right){\tau }_1{I}_h}{w_2},1\right\}\right\},\\ u_3^{*}=max\left\{0,min\left\{\frac{\left({\Gamma }_5-{\Gamma }_3\right){\tau }_2{I}_h}{w_3},1\right\}\right\}& \left(26\right)\end{array}$

Proof: By using [50], the integrand of the objective functional J(u₁, u₂, u₃) is a convex function of u₁, u₂ and u₃. Since the solution of the system (Equation 19) is bounded, hence the system satisfies the Lipshitz property with respect to the variables S_h, E_h, I_h, T_h, T_s, R_h, S_m, E_m and I_m. Therefore, there exists an optimal pair$\left({u_1}^{*},u_2^{*},u_3^{*}\right)$. The governing equations of the adjoint variables are obtained by differentiation of the Hamiltonian function, evaluated at the optimal control. Then, the adjoint system can be written as:

$\begin{array}{l}\frac{d{\Gamma }_1}{dt}=-\frac{\partial \mathcal{H}}{\partial S_h}=-{\Gamma }_1\left(\left(1-u_1\right){\lambda }_h+\mu \right)+{\Gamma }_2\left(1-u_1\right)\text{Or}\\ {\Gamma }_1^{\prime }={\Gamma }_1\left(\left(1-u_1\right){\lambda }_h+\mu \right)-{\Gamma }_2\left(1-u_1\right)\\ \frac{d{\Gamma }_2}{dt}=-\frac{\partial \mathcal{H}}{\partial E_h}=D_1-{\Gamma }_2\left({\theta }_h+\mu \right)+{\Gamma }_3{\theta }_h\\ \frac{d{\Gamma }_3}{dt}=-\frac{\partial \mathcal{H}}{\partial I_h}=D_2-{\Gamma }_3\left(\delta +u_2{\tau }_1+\left(1-u_3\right){\tau }_2+\mu \right)\\ +{\Gamma }_4\left(u_2{\tau }_1\right)+{\Gamma }_5\left(1-u_3\right){\tau }_2+\left({\Gamma }_8-{\Gamma }_7\right)\frac{\left(1-u_1\right)q{\alpha }_{1m}{S}_m}{N_h}\\ \frac{d{\Gamma }_4}{dt}=-\frac{\partial \mathcal{H}}{\partial T_h}={\Gamma }_1\rho -{\Gamma }_4\left(\rho +\mu \right)\end{array}$

$\begin{array}{l}\frac{d{\Gamma }_5}{dt}=-\frac{\partial \mathcal{H}}{\partial T_s}={\Gamma }_4\omega -{\Gamma }_5\left(\omega +\gamma +\mu +\psi \right)+{\Gamma }_6\gamma \\ +\left({\Gamma }_8-{\Gamma }_7\right)\frac{\left(1-u_1\right)q{\alpha }_{2m}{S}_m}{N_h}& \left(27\right)\end{array}$

$\begin{array}{l}\frac{d{\Gamma }_6}{dt}=-\frac{\partial \mathcal{H}}{\partial R_h}={\Gamma }_1{\rho }_1+{\Gamma }_3{\gamma }_1-{\Gamma }_6\left({\gamma }_1+{\rho }_1+\mu \right)\\ \frac{d{\Gamma }_7}{dt}=-\frac{\partial \mathcal{H}}{\partial S_m}=D_3+{\Gamma }_8\left(1-u_1\right){\lambda }_m-{\Gamma }_7\left(\left(1-u_1\right){\lambda }_m+\eta \right)\\ \frac{d{\Gamma }_8}{dt}=\frac{\partial \mathcal{H}}{\partial E_m}=D_3+{\Gamma }_9{\theta }_m-{\Gamma }_8\left({\theta }_m+\eta \right)\\ \frac{d{\Gamma }_9}{dt}=-\frac{\partial \mathcal{H}}{\partial I_m}=D_3+\left({\Gamma }_2-{\Gamma }_1\right)\left(1-u_1\right)q{\beta }_{mh}-{\Gamma }_9\eta \end{array}$

With transversality conditions Γ₁(t_f) = Γ₂(t_f) = Γ₃(t_f) = Γ₄(t_f) = Γ₅(t_f) = Γ₆(t_f) = Γ₇(t_f) = Γ₈(t_f) = Γ₉(t_f) = 0.

Then, from the optimality condition, we have:

$\begin{array}{l}\frac{\partial \mathcal{H}}{\partial u_i}=0,at{u}_1=u_i^{*}fori=1,2,3\end{array}$

Using [51], we get, $u_1=\frac{\left({\Gamma }_2-{\Gamma }_1\right){\lambda }_h{S}_h+\left({\Gamma }_8-{\Gamma }_7\right){\lambda }_m{S}_m}{w_1}$, $u_2=\frac{\left({\Gamma }_3-{\Gamma }_4\right){\tau }_1{I}_h}{w_2}$, $u_3=\frac{\left({\Gamma }_5-{\Gamma }_3\right){\tau }_2{I}_h}{w_3}$ and solving for $u_1^{*}$, $u_2^{*}$ and $u_3^{*}$ we get:

$\begin{array}{l}{u_1}^{*}=max\left\{0,min\left\{\frac{\left({\Gamma }_2-{\Gamma }_1\right){\lambda }_h{S}_h+\left({\Gamma }_8-{\Gamma }_7\right){\lambda }_m{S}_m}{w_1},1\right\}\right\},\\ u_2^{*}=max\left\{0,min\left\{\frac{\left({\Gamma }_3-{\Gamma }_4\right){\tau }_1{I}_h}{w_2},1\right\}\right\}\\ \text{and }{u}_3^{*}=max\left\{0,min\left\{\frac{\left({\Gamma }_5-{\Gamma }_3\right){\tau }_2{I}_h}{w_3},1\right\}\right\}\\ \end{array}$

Additionally, the sufficient condition (the second derivative of the Hamilton function with respect to u₁, u₂ and u₃). i.e. $\frac{{\partial }^2\mathcal{H}}{\partial u_1^2}=0\Rightarrow \frac{{\partial }^2\mathcal{H}}{\partial u_1^2}=w_1$, $\frac{{\partial }^2\mathcal{H}}{\partial u_2^2}=0\Rightarrow \frac{{\partial }^2\mathcal{H}}{\partial u_2^2}=w_2$ and $\frac{{\partial }^2\mathcal{H}}{\partial u_3^2}=0\Rightarrow \frac{{\partial }^2\mathcal{H}}{\partial u_3^2}=w_3$.

This shows that the optimal problem goes minimization at u₁, u₂ and u₃.

6.1 Numerical simulations of optimal control

The optimality system consists of the state system in the system (Equation 1), optimal control set in Equation 19, adjoint system in Equation 24, boundary conditions in Equation 25, and initial conditions S_h(0) = 85, E_h(0) = 20, I_h(0) = 20, T_h(0) = 5, T_s(0) = 10, R_h(0) = 5, S_m(0) = 1, 000, E_m(0) = 20, and I_m(0) = 30, [16], while the parameter values are given in Table 2. We carried out the numerical simulation with MATLAB R2023a by using the forward-backward fourth-order Runge-Kutta iterative method to solve the optimality system [46]. The state (Equation 19) is solved with the initial values of state variables using the fourth-order forward Runge–Kutta method and then we used backward fourth-order Runge-Kutta to solve the adjoint equations [45, 46, 51]. To determine the impact of control measures on the reduction of malaria infection, we have the following cases of optimal control strategies:

Intervention I: Optimal use of ITNs, treatment, and awareness creation toward avoiding inappropriate traditional medicines

Figures 7a–c depicts when the three control mechanisms, namely ITNs (u₁ ≠ 0), antimalarial drug treatment (u₂ ≠ 0), and awareness creation toward avoiding un-prescribed TM (u₃ ≠ 0) are implemented to optimize the objective function J. The results show that the number of exposed humans, infectious humans, and mosquitoes are decreased more rapidly when controls are in use than the case without controls. We found that while ITN use, treatment, and raising awareness against the use of un-prescribed traditional medicines can all contribute to a reduction in the malaria burden. Combining these interventions is the most effective way to get improved outcomes

Figure 7

Figure 7. Effect of intervention –I on I_h and I_m over time. (a) Simulation of optimal use of u₁, u₂ and u₃ on infectious humans. (b) Simulation of optimal use of all interventions on infectious humans when TM is not used. (c) Simulation of optimal use of u₁, u₂ and u₃ on infectious mosquitoes.

Intervention II: Optimal use of ITNs and antimalarial drug treatment

In this case, simulation is done for infectious humans and mosquitoes when there is no control strategy in place and when there are protective efforts using ITNs and treatment for infectious individuals. As Figures 8a, b illustrates, when controls are in place, the number of infectious humans and infectious mosquitoes decline more quickly than the case without controls, respectively.

Figure 8

Figure 8. Effect of intervention II on I_h and I_m over time. (a) Simulation of optimal use of u₁ and u₂ on infectious humans. (b) Simulation of optimal use of u₁ and u₂ on infectious mosquitoes.

Intervention III: Optimal use of ITNs and awareness creation toward inappropriate use of traditional medicines

In this case, simulation is done for infectious humans and mosquitoes when there is a control strategy in place and when there are combined protective efforts using ITNs and using media campaigns for awareness creation of avoiding unprescribed traditional medicines. As Figures 9a, b illustrate, when controls are in place, the number of infectious humans and mosquitoes declines more quickly than in the case without controls.

Figure 9

Figure 9. Effect of intervention III on I_h and I_m over time. (a) Simulation of optimal use of u₁ and u₃ on infectious humans. (b) Simulation of optimal use of u₁ and u₃ on infectious mosquitoes.

Intervention IV: Optimal use of treatment of Infectious Humans and Awareness creation toward inappropriate use of traditional medicines

In this case, simulation is done for infectious humans and mosquitoes when there is no control strategy in place and when there are combined effects of protective efforts using treatment for infectious individuals and using media campaigns to raise awareness of avoiding clinically not valid traditional medicines. As shown in Figures 10a, b, when controls are in place, the number of infectious humans and mosquitoes declines more rapidly than in the case without controls. Here, on comparison, we find that using strategy IV is the powerful one to mitigate malaria compared to the strategy III because it reduces the infectious individuals and mosquitoes to zero after 100 months.

Figure 10

Figure 10. Effect of intervention IV on I_h and I_m over time. (a) Simulation of optimal use of u₂ and u₃ on infectious humans. (b) Simulation of optimal use of u₂ and u₃ on infectious mosquitoes.

Intervention V: Optimal use of ITNs

Figures 11a, b illustrates how the adoption of personal protective measures using ITNs alone has a positive effect on curbing the prevalence of malaria within the community, but it is not sufficient to eliminate the disease.

Figure 11

Figure 11. Effect of intervention V on I_h and I_m over time. (a) Simulation of optimal use of u₁ on infectious humans. (b) Simulation of optimal use of u₁ on infectious mosquitoes.

Intervention VI: Optimal use of antimalarial drug treatment

The implementation of treatment alone can effectively slow the transmission of malaria within a community, but it is insufficient to eradicate the illness, as shown by Figures 12a, b.

Figure 12

Figure 12. Effect of intervention VI on I_h and I_m over time. (a) Simulation of optimal use of u₃ on infectious humans. (b) Simulation of Optimal use of u₃ on infectious mosquitoes.

Intervention VII: Optimal use of awareness creation toward inappropriate use of TMs

The implementation of awareness creation toward unprescribed traditional medicines alone can reduce the spread of malaria within a community, though it is insufficient to eradicate the illness, as shown by Figures 13a, b.

Figure 13

Figure 13. Effect of intervention VII on I_h and I_m over time. (a) Simulation of optimal use of u₃ on infectious humans. (b) Simulation of optimal use of u₃ on infectious mosquitoes.

Moreover, when there is no treatment with inappropriate TMs (τ₂ = 0), the three control mechanisms strictly decrease the infected humans initially from 30 to 0 at 20 months, and when there is use of inappropriate TMs (τ₂ > 0), the three control mechanisms decrease the infectious humans initially from 30 to 0 after 100 months, as shown in Figures 9a, b, respectively. Although the three strategies do lower the prevalence of malaria, when traditional medicines are used for treatment, the disease is still endemic in the community. Therefore, using the three control mechanisms at the same time, avoiding treatment with inappropriate traditional medicines, and seeking treatment by health professionals at health facilities are the best strategies to mitigate malaria.

7 Discussion and conclusion

This study presents a comprehensive mathematical model to analyze the impact of treatment-seeking behaviors using optimal control on the transmission dynamics of malaria. The model incorporates distinct groups of individuals seeking treatment at health facilities and those relying on traditional remedies, which lack clinical validation. By employing a new deterministic model, the study explores the effects of various control measures, including antimalarial drug treatment, personal protective measures like insecticide-treated bed nets (ITNs), and awareness campaigns to discourage the inappropriate use of traditional medicine. The findings highlight that the utilization of health facilities for malaria treatment significantly reduces the basic reproduction number, underscoring the importance of effective healthcare interventions in malaria-endemic regions. Sensitivity analysis reveals that parameters such as the biting rate of mosquitoes and the probability of transmission from individuals undergoing traditional medicine treatment to susceptible mosquitoes are highly sensitive to changes in reproduction number. This suggests that interventions targeting these parameters could be particularly effective in reducing malaria transmission. Moreover, from our sensitivity results, we can observe that the value of the transmission of infection from I_h to S_m(α_1m = 0.0990734), transmission of infection from T_s to S_m(α_2m = 0.359975) and transmission of infection from I_m to S_h (β_mh = 0.5) implies that they increase the reproduction number. This means that the positive parameter values increase the reproduction number when they are increased, while the negative parameter values decrease the reproduction number when they are increased.

Numerical simulations demonstrate that increasing the treatment rate at health facilities leads to a rapid decline in the infectious human population, potentially eliminating the disease. In contrast, reliance on traditional medicine without proper clinical validation can prolong the disease's persistence and can take a long time to eradicate. The study also explores the role of awareness campaigns in reducing the use of inappropriate traditional medicines. Increasing the progression rate from traditional medicine treatment to health facility treatment significantly reduces the population relying on traditional remedies, thereby decreasing the infectious human population. This highlights the importance of public health campaigns in guiding individuals toward clinically validated treatments. Optimal control analysis further emphasizes that a combination of control strategies—ITNs, antimalarial drug treatment, and awareness campaigns—is the most effective approach to combating malaria. The simulations show that implementing all three controls simultaneously leads to a faster reduction in the number of infected individuals and mosquitoes compared to using any single control measure. This finding aligns with previous studies that advocate for integrated control strategies to achieve malaria eradication.

The study highlights the importance of healthcare facilities in malaria control, highlighting the need for increased treatment rates to reduce the basic reproduction number and accelerate the decline in the infectious human population. Traditional medicine, widely used in malaria-endemic regions, faces challenges due to its lack of clinical validation and standardized dosage. Relying on traditional remedies can prolong malaria's persistence, emphasizing the need for consultation with healthcare professionals. Public health campaigns to discourage inappropriate use of traditional medicine and promote clinically validated treatments are crucial. Integrated control strategies, including the use of insecticide-treated bed nets (ITNs), antimalarial drug treatment, and awareness campaigns, are the most effective approach to malaria control. Implementing these strategies simultaneously leads to a faster and more substantial reduction in malaria transmission compared to using any single intervention.

The study suggests that public health policies in malaria-endemic regions should focus on strengthening healthcare infrastructure, regulating traditional medicine, intensifying public health campaigns, and adopting integrated malaria control programs. These include expanding healthcare facilities, ensuring access to antimalarial treatments, and promoting collaboration between traditional practitioners and modern healthcare providers. This holistic approach is more likely to achieve sustainable malaria eradication.

8 Limitations of the study

While this Study provides valuable insights, it has some limitations. The model assumes homogeneous mixing within the population, which may not fully capture the complexities of real-world malaria transmission. Additionally, the study relies on parameter values from existing literature, which may not fully reflect the specific conditions of all malaria-endemic regions. Finally, this study is considered deterministic and governed by a system of integer-order ordinary differential equations. Considering the dynamics of malaria transmission with fractional-order approaches to investigate the memory effect and cost-effectiveness analysis is our future study.

Data availability statement

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

Ethics statement

Written informed consent was obtained from all individuals and, where applicable, from the legal guardians or next of kin of minors, for the publication of any potentially identifiable images or data included in this article.

Author contributions

SJ: Conceptualization, Data curation, Funding acquisition, Investigation, Methodology, Resources, Software, Visualization, Writing – original draft. GD: Formal analysis, Investigation, Methodology, Supervision, Validation, Visualization, Writing – review & editing. CD: Conceptualization, Formal analysis, Project administration, Supervision, Validation, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Conflict of interest

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

Generative AI statement

The author(s) declare that no Gen AI was used in the creation of this manuscript.

Publisher's note

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

Supplementary material

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

Footnotes

1. ^https://www.cdc.gov/dpdx/malaria/index.html

References

1. Pucca MB, de Sousa TN, Cardoso de Melo G, Viana GMR. Challenges for diagnosis, treatment, and elimination of malaria. Front Trop Dis. (2024) 5:1394693. doi: 10.3389/fitd.2024.1394693

Crossref Full Text | Google Scholar

2. World Health Organization. (2023). World malaria report 2023. World Health Organization.

Google Scholar

3. Life Cycle of Culex Mosquitoes-CDC. Available online at: https://www.cdc.gov/mosquitoes/about/life-cycle-of-culex-mosquitoes.html#:~:text=Adult%20female%20Culex%20mosquitoes%20lay,develop%20into%20adult%2C%20flying%20mosquitoes

Google Scholar

4. Urama CE, Manasseh CO, Ukwueze ER, Ogbuabor JE. Choices and determinants of malaria treatment seeking behaviour by rural households in Enugu state, South-East Nigeria. Int J Health Promot Educ. (2021) 59:156–73.

Google Scholar

5. Kumar A, Kumar M, Nilam A. A study on the stability behavior of an epidemic model with ratio- dependent incidence and saturated treatment. Theory in Biosci. (2020) 139:225–34. doi: 10.1007/s12064-020-00314-6

PubMed Abstract | Crossref Full Text | Google Scholar

6. World Health Organization. Global Technical Strategy for Malaria 2016–2030. World Health Organization (2015).

Google Scholar

7. Hsu E. Reflections on the ‘discovery' of the antimalarial qinghao. Br J Clin Pharmacol. (2006) 61:666–70. doi: 10.1111/j.1365-2125.2006.02673.x

PubMed Abstract | Crossref Full Text | Google Scholar

8. Teng WC, Kiat HH, Suwanarusk R, Koh HL. Medicinal Plants and Malaria: Applications, Trends, Prospects. Boca Raton, FL: CRC Press. (2016).

Google Scholar

9. Akilimali A, Bisimwa C, Aborode AT, Biamba C, Sironge L, Balume A, et al. Self-medication and anti-malarial drug resistance in the Democratic Republic of the Congo (DRC): a silent threat. Trop Med Health. (2022) 50:73. doi: 10.1186/s41182-022-00466-9

PubMed Abstract | Crossref Full Text | Google Scholar

10. Noronha M, Pawar V, Prajapati A, Subramanian RB. A literature review on traditional herbal medicines for malaria. South Afric J Botany. (2020) 128:292–303. doi: 10.1016/j.sajb.2019.11.017

Crossref Full Text | Google Scholar

11. Benzie IFF, Watchel-Galor S. Herbal Medicine: an introduction to its history, usage, regulation, current trends, research needs. Herbal Med. (2011) 464:b10787. doi: 10.1201/b10787

PubMed Abstract | Crossref Full Text | Google Scholar

12. Appiah EO, Appiah S, Oti-Boadi E, Oppong-Besse A, Awuah DB, Asiedu PO, et al. Practices of herbal management of malaria among trading mothers in Shai Osudoku District, Accra. PLoS ONE. (2022) 17:e0271669. doi: 10.1371/journal.pone.0271669

PubMed Abstract | Crossref Full Text | Google Scholar

13. Martcheva M. An Introduction to Mathematical Epidemiology. New York: Springer. (2015). p. 9–31.

Google Scholar

14. Ross R. The Prevention of Malaria. London: John Murray. (1911).

Google Scholar

15. Gebremichael SM, Mekonnen TT. Relapse effect on the dynamics of malaria in humans and mosquitoes: a mathematical model analysis. IOSR J Mathemat. (2019) 15:46–59.

Google Scholar

16. Olaniyi S, Falowo OD, Okosun KO, Mukamuri M, Obabiyi OS, Adepoju OA. Effect of saturated treatment on malaria spread with optimal intervention. Alexandria Eng J. (2023) 65:443–59. doi: 10.1016/j.aej.2022.09.024

Crossref Full Text | Google Scholar

17. Woldegerima WA, Ouifki R, Banasiak J. Mathematical analysis of the impact of transmission-blocking drugs on the population dynamics of malaria. Appl Math Comput. (2021) 400:126005. doi: 10.1016/j.amc.2021.126005

Crossref Full Text | Google Scholar

18. Bakary T, Boureima S, Sado T. A mathematical model of malaria transmission in a periodic environment. J. Biol. Dynam. (2018) 12:400–32. doi: 10.1080/17513758.2018.1468935

PubMed Abstract | Crossref Full Text | Google Scholar

19. Elakhe OA, Isere AO, Akerejola RF. Mathematical model of malaria transmission with anti-malarial herbal therapy as control. Studies (2023) 6:1–16. doi: 10.52589/AJMSS-1276JR4U

PubMed Abstract | Crossref Full Text | Google Scholar

20. Haringo AT, Obsu LL, Bushu FK. A mathematical model of malaria transmission with media-awareness and treatment interventions. J Appl Mathemat Comp. (2024) 2024:1–39. doi: 10.1007/s12190-024-02154-9

Crossref Full Text | Google Scholar

21. Collins OC, Duffy KJ. A mathematical model for the dynamics and control of malaria in Nigeria. Infectious Dis Model. (2022) 7:728–41. doi: 10.1016/j.idm.2022.10.005

PubMed Abstract | Crossref Full Text | Google Scholar

22. Abioye AI, Ibrahim MO, Peter OJ, Ogunseye HA. Optimal control on a mathematical model of malaria. Sci Bull, Series A: Appl Math Phy. (2020) 82:177–90.

Google Scholar

23. Keno TD, Dano LB, Makinde OD. Modeling optimal control analysis for malaria transmission with role of climate variability. Computat Mathemat Methods. (2022) 2022:9667396. doi: 10.1155/2022/9667396

Crossref Full Text | Google Scholar

24. Olaniyi S, Okosun KO, Adesanya SO, Lebelo RS. Modelling malaria dynamics with partial immunity and protected travellers: optimal control and cost-effectiveness analysis. J Biol Dyn. (2020) 14:90–115. doi: 10.1080/17513758.2020.1722265

PubMed Abstract | Crossref Full Text | Google Scholar

25. Tchoumi SY, Chukwu CW, Diagne ML, Rwezaura H, Juga ML, Tchuenche JM. Optimal control of a two-group malaria transmission model with vaccination. Netw Model Analy Health Inform Bioinform. (2022) 12:7. doi: 10.1007/s13721-022-00403-0

PubMed Abstract | Crossref Full Text | Google Scholar

26. Graz B, Kitua AY, Malebo HM. To what extent can traditional medicine contribute a complementary or alternative solution to malaria control programmes?. Malaria J. (2011) 10:1–7. doi: 10.1186/1475-2875-10-S1-S6

PubMed Abstract | Crossref Full Text | Google Scholar

27. Prakash Raj M, Venkatesh A, Arun Kumar K, Manivel M. Mathematical modeling of the co-infection dynamics of dengue and malaria using delay differential equations. Adv Theor Simul. (2024) 2400609.

Google Scholar

28. Jin X, Jin S, Gao D. Mathematical analysis of the Ross–Macdonald model with quarantine. Bull Math Biol. (2020) 82:47. doi: 10.1007/s11538-020-00723-0

PubMed Abstract | Crossref Full Text | Google Scholar

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

Crossref Full Text | Google Scholar

30. Jan R, Razak NNA, Boulaaras S, Rehman ZU, Bahramand S. Mathematical analysis of the transmission dynamics of viral infection with effective control policies via fractional derivative. Nonlinear Eng. (2023) 12:20220342. doi: 10.1515/nleng-2022-0342

Crossref Full Text | Google Scholar

31. Deebani W, Jan R, Shah Z, Vrinceanu N, Racheriu M. Modeling the transmission phenomena of water-borne disease with non-singular and non-local kernel. Comput Methods Biomech Biomed Engin. (2023) 26:1294–307. doi: 10.1080/10255842.2022.2114793

PubMed Abstract | Crossref Full Text | Google Scholar

32. Perko L. Differential Equations and Dynamical Systems. Cham: Springer Science and Business Media. (2013).

Google Scholar

33. Tumwiine J, Mugisha JYT, Luboobi LS. A mathematical model for the dynamics of malaria in a human host and mosquito vector with temporary immunity. Appl Math Comput. (2007) 189:1953–65. doi: 10.1016/j.amc.2006.12.084

PubMed Abstract | Crossref Full Text | Google Scholar

34. Alshehri A, Shah Z, Jan R. Mathematical study of the dynamics of lymphatic filariasis infection via fractional-calculus. Eur Phys J Plus. (2023) 138:1–15. doi: 10.1140/epjp/s13360-023-03881-x

PubMed Abstract | Crossref Full Text | Google Scholar

35. Castillo-Chavez C, Feng Z, Huang W. On the computation of R0 and its role on global stability. In: Castillo-Chávez C, van den Dreissche P, Kirschner D, Yakubu A, editor. Mathematical Approaches for Emerging and Reemerging Infectious Diseases: An Introduction. (2002).

Google Scholar

36. Jan R, Boulaaras S, Alyobi S, Jawad M. Transmission dynamics of Hand–Foot–Mouth Disease with partial immunity through non-integer derivative. Int J Biomathem. (2023) 16:2250115. doi: 10.1142/S1793524522501157

Crossref Full Text | Google Scholar

37. Allen LJ, Brauer F, Van den Driessche P, Wu J. Mathematical Epidemiology. Berlin: Springer. (2008).

Google Scholar

38. Jan R, Hinçal E, Hosseini K, Razak NNA, Abdeljawad T, Osman MS. Fractional view analysis of the impact of vaccination on the dynamics of a viral infection. Alexandria Eng J. (2024) 102:36–48. doi: 10.1016/j.aej.2024.05.080

PubMed Abstract | Crossref Full Text | Google Scholar

39. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | Crossref Full Text | Google Scholar

40. Deressa CT, Mussa YO, Duressa GF. Optimal control and sensitivity analysis for transmission dynamics of Coronavirus. Results Phys. (2020) 19:103642. doi: 10.1016/j.rinp.2020.103642

PubMed Abstract | Crossref Full Text | Google Scholar

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

PubMed Abstract | Crossref Full Text | Google Scholar

42. Chitnis N, Hyman JM, Cushing JM. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bull Math Biol. (2008) 70:1272–96. doi: 10.1007/s11538-008-9299-0

PubMed Abstract | Crossref Full Text | Google Scholar

43. Mangongo YT, Bukweli JDK, Kampempe JDB, Mabela RM, Munganga JMW. Stability global sensitivity analysis of the transmission dynamics of malaria with relapse and ignorant infected humans. Physica Scripta. (2022) 97:024002. doi: 10.1088/1402-4896/ac4862

Crossref Full Text | Google Scholar

44. Aldila D, Angelina M. Optimal control problem and backward bifurcation on malaria transmission with vector bias. Heliyon. (2021) 7:4. doi: 10.1016/j.heliyon.2021.e06824

PubMed Abstract | Crossref Full Text | Google Scholar

45. Kotola BS Teklu SW Abebaw Abebaw YF Bifurcation optimal control analysis of HIV/AIDS and COVID-19 co-infection model with numerical simulation. PLoS ONE. (2023) 18:e0284759. doi: 10.1371/journal.pone.0284759

PubMed Abstract | Crossref Full Text | Google Scholar

46. Okosun KO, Makinde OD. A co-infection model of malaria and cholera diseases with optimal control. Math Biosci. (2014) 258:19–32. doi: 10.1016/j.mbs.2014.09.008

PubMed Abstract | Crossref Full Text | Google Scholar

47. Ouattara L, Ouedraogo D, Diop O, Guiro A. Analysis and optimal control of a mathematical model of malaria. Nonlinear Dyn Syst Theory. (2024) 24.

Google Scholar

48. Pontryagin LS. Mathematical Theory of Optimal Processes. London: Routledge. (2018).

Google Scholar

49. Agusto FB. Optimal control and temperature variations of malaria transmission dynamics. Complexity. (2020) 2020:5056432. doi: 10.1155/2020/5056432

Crossref Full Text | Google Scholar

50. Fleming WH, Rishel RW. Deterministic and Stochastic Optimal Control. Cham: Springer Science and Business Media. (2012).

Google Scholar

51. Lenhart S, Workman JT. Optimal Control Applied to Biological Models. Boca Raton, FL: Chapman and Hall/CRC. (2007).

Google Scholar