Modeling the Impact of Human Awareness and Insecticide Use on Malaria Control: A Fractional-Order Approach

Abstract

In this research work, we developed a fractional-order model for the transmission dynamics of malaria, incorporating two control strategies: health education campaigns and the use of insecticides. The theoretical analysis of the model is presented, including the computation of disease-free equilibrium and basic reproduction number. We analyzed the stability of the proposed model using a well-formulated Lyapunov function. Furthermore, model parameter estimation was carried out using real data from malaria cases reported in Zimbabwe. We found that the fractional-order model provided a better fit to the real data compared to the classical integer-order model. Sensitivity analysis of the basic reproduction number was performed using computed partial rank correlation coefficients to assess the effect of each parameter on malaria transmission. Additionally, we conducted numerical simulations to evaluate the impact of memory effects on the spread of malaria. The simulation results indicated that the order of derivatives significantly influences the dynamics of malaria transmission. Moreover, we simulated the model to assess the effectiveness of the proposed control strategies. Overall, the interventions were found to have the potential to significantly reduce the spread of malaria within the population.

1. Introduction

Malaria is a vector-borne disease caused by Plasmodium parasites, which are transmitted to humans through the bites of infected female Anopheles mosquitoes [1]. The primary mosquito species responsible for malaria transmission are Anopheles gambiae and Anopheles funestus, which play a significant role in the malaria transmission cycle [2]. Malaria was identified as a mosquito-transmitted disease in the early 20th century and remains a widespread and dangerous infectious disease, especially in tropical and subtropical regions [3]. Like other vector-borne diseases such as dengue, Zika virus, and yellow fever, malaria is transmitted when an infected mosquito bites a susceptible individual during a blood meal [4]. Infected mosquitoes acquire the parasite by biting an infected individual, though they do not suffer from the disease themselves [5]. Malaria’s incubation period in humans typically ranges from 7 to 30 days, depending on the specific Plasmodium species [6]. Symptoms include fever, chills, headache, muscle pain, and fatigue; in severe cases, it can lead to anemia, organ failure, and death [7]. While some malaria vaccines have been developed with limited success, no universal vaccine exists, making disease management reliant on antimalarial medications and preventive measures, such as insecticide-treated bed nets, indoor residual spraying, and mosquito repellents [8].

Mathematical models are valuable in epidemiology, providing insights into malaria’s transmission dynamics and supporting public health officials in developing control strategies [9,10,11]. By using mathematical models, researchers can simulate various intervention scenarios, study potential control outcomes, and understand malaria’s behavior in different populations and environments [12,13,14]. Fractional calculus, which applies non-integer-order derivatives, has been widely used to model real-world problems, including heat conduction, control theory, chaotic systems, and biological processes [15,16,17]. The literature suggests that non-integer derivatives yield higher accuracy in modeling biological and physical systems than integer-order derivatives [18]. In disease dynamics, fractional-order derivatives are particularly advantageous for capturing memory effects and hereditary properties within biological systems [19]. Given these advantages, fractional-order derivatives have attracted researchers who want to develop more dynamic and efficient mathematical models [20]. Researchers have applied fractional operators to model various epidemics, incorporating different scenarios and constraints [10,17,21]. Unlike integer-order derivatives, fractional-order derivatives capture biological systems’ memory effects and hereditary properties [22], and research indicates that cell membranes possess fractional-order electrical conductance [23], aligning with fractional-order modeling. Moreover, fractional-order derivatives often provide a better fit to real epidemiological data [24]. The most widely used fractional derivatives include the Caputo, Riemann-Liouville, and Atangana-Baleanu derivatives [15].

Recently, mathematical models using fractional derivatives have been explored to study diseases. For instance, [25] applied the Atangana-Baleanu operator to examine personal protective measures in malaria transmission. Although the model was not validated with real data, it suggested that reducing the fractional-order derivative value from 1 to 0 decreased infections in both humans and mosquitoes. In [26], a fractional-order model for malaria and COVID-19 co-infection was developed, calculating the threshold quantity and using the fixed-point theorem to analyze the model; this showed that preventive measures can reduce disease spread. Another fractional-order malaria model was presented in [10], utilizing the Euler and Adam-Bashforth-Moulton schemes for simulation. The results indicated that fractional-order derivatives significantly influence malaria dynamics. For further information on vector-borne disease models using fractional derivatives, readers can refer to [10,16,21,24,27,28,29] and related references.

In this study, we propose a fractional-order model of malaria transmission to explore the effects of human awareness on insecticide use. Although many fractional derivatives are available, our model uses the Caputo derivative due to its advantages, such as being zero for constant functions (consistent with integer-order derivatives) and allowing for the inclusion of local initial conditions [15]. Recent studies have shown that increasing public awareness of insecticide use is crucial for malaria prevention [29]. Health education campaigns that raise awareness of insecticide application may significantly impact malaria prevention [30]. To the best of our knowledge, few mathematical models for malaria transmission have incorporated real reported malaria case data. This study proposes a fractional-order malaria model incorporating two control strategies—health education campaigns and insecticide use—to assess their impact on reducing malaria spread in the population.

The rest of the article is organized as follows: in Section 2, the mathematical model formulation is presented. In Section 3, the basic properties of the model—the positivity of the variables and the boundedness of the trajectories—are presented, and the computation of the reproduction number and the existence of model equilibrium are presented in Section 4. The results and discussion are presented in Section 5, and our conclusion remarks complete the paper in Section 6.

2. Model Formulation

In this study, the fractional-order model of malaria disease dynamics based on Caputo has been formulated and studied. The proposed model demonstrates the interplay between mosquito and human populations. Throughout the document, the subscripts h and v, respectively, denote the mosquito and human populations. The human population is subdivided into four groups based on infection status, namely the susceptible $S_h\left(t\right)$, exposed $E_h\left(t\right)$, infectious $I_h\left(t\right)$, and recovered $R_h\left(t\right)$ classes. Thus, the total population of humans at time t is denoted by $N_h\left(t\right)$, given as

$$\begin{array}{c}\hfill N_h\left(t\right)=S_h\left(t\right)+E_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)\end{array}$$

(1)

Furthermore, the total population of mosquitoes at time t is denoted by $N_v\left(t\right)$, which is subdivided into three groups, namely, the susceptible $S_v\left(t\right)$, exposed $E_v\left(t\right)$, and infectious $I_v\left(t\right)$ classes. Thus, the total population of mosquitoes at time t is defined as

$$\begin{array}{c}\hfill N_v\left(t\right)=S_v\left(t\right)+E_v\left(t\right)+I_v\left(t\right)\end{array}$$

(2)

All the parameters and variables of the proposed model are considered to be non-negative, and the parameters are described as follows: ${\mathrm{\Lambda }}_h$ and ${\mathrm{\Lambda }}_v$ denote the new recruitment rate of humans and vectors, respectively, and they are all considered to be susceptible whenever coming from either birth or immigrants; ${\mu }_h$ and ${\mu }_v$ denote the natural mortality rate of humans and vectors, respectively; ${\alpha }_h$ and ${\alpha }_v$ denote the transfer rate of humans and vectors, respectively, from the exposed class to the infectious class; ${\theta }_h$ represents the natural recovery rate of exposed individuals; ${\gamma }_h$ denotes the recovery rate of infected humans; we have assumed that both recovered individuals (from exposed and infected classes) recover with temporary immunity. Thus, the parameter ${\rho }_h$ represents the waning rate of immunity for recovered humans; $d_h$ represents the disease-induced death rate of infected humans; ${\beta }_h$ denotes the probability of disease transmission from infected mosquito to susceptible human due to successful contact between infected vectors and susceptible humans; ${\beta }_v$ denotes the probability of disease transmission from infected humans to susceptible mosquitoes due to successful contact between infected persons and susceptible vectors; ${\delta }_v$ denotes the mosquito biting rate on humans. Recently, studies have shown that one of the possible ways of minimizing the spread of malaria disease in the population is based on how much the community is aware of the use of prevention measures, such as the use of mosquito nets and insecticides, that reduce the possibility of human contact with mosquitoes. Therefore, in this study, we consider two control strategies, namely, the use of prevention measures such as mosquito nets, denoted by $\omega$, and the use of insecticides, represented by $ϵ$. Based on the above assumptions, we have assumed the following flowchart and nonlinear differential equations.

From the flowchart diagram in Figure 1, the dynamics of the interaction between humans and mosquitoes is given by the following system of ordinary differential equations:

$$\left\{\begin{array}{ccc}{\displaystyle {}_{t_0}{}^{c}D_t^{\varphi }{S}_h\left(t\right)}\hfill & =\hfill & {\displaystyle {\mathrm{\Lambda }}_h+{\rho }_h{R}_h-(1-\omega ){\delta }_v{\beta }_h{S}_h{I}_v-{\mu }_h{S}_h}\hfill \\ {\displaystyle {}_{t_0}{}^{c}D_t^{\varphi }{E}_h\left(t\right)}\hfill & =\hfill & {\displaystyle (1-\omega ){\delta }_v{\beta }_h{S}_h{I}_v-({\mu }_h+{\alpha }_h+{\theta }_h)E_h,}\hfill \\ {\displaystyle {}_{t_0}{}^{c}D_t^{\varphi }{I}_h\left(t\right)}\hfill & =\hfill & {\displaystyle {\alpha }_h{E}_h-({\mu }_h+d_h+{\gamma }_h)I_h,}\hfill \\ {\displaystyle {}_{t_0}{}^{c}D_t^{\alpha }{R}_h\left(t\right)}\hfill & =\hfill & {\displaystyle {\theta }_h{E}_h+{\gamma }_h{I}_h-({\rho }_h+{\mu }_h)R_h,}\hfill \\ {\displaystyle {}_{t_0}{}^{c}D_t^{\varphi }{S}_v\left(t\right)}\hfill & =\hfill & {\displaystyle {\mathrm{\Lambda }}_v-(1-\omega ){\delta }_v{\beta }_v{I}_h{S}_v-({\mu }_v+ϵ)S_v,}\hfill \\ {\displaystyle {}_{t_0}{}^{c}D_t^{\varphi }{E}_v\left(t\right)}\hfill & =\hfill & {\displaystyle (1-\omega ){\delta }_v{\beta }_v{I}_h{S}_v-({\mu }_v+{\alpha }_v+ϵ)E_v,}\hfill \\ {\displaystyle {}_{t_0}{}^{c}D_t^{\varphi }{I}_v\left(t\right)}\hfill & =\hfill & {\displaystyle {\alpha }_v{E}_v-({\mu }_v+ϵ)I_v.}\hfill \end{array}\right.$$

(3)

Preliminaries on the Caputo Fractional Calculus

We begin by introducing the definition of the Caputo fractional derivative and state-related theorems (see [31,32,33,34]) that we utilized to derive important results in this work.

Definition 1. 

Suppose that $\varphi >0,t>a,\varphi ,a,t\in \mathbb{R}.$ The Caputo fractional derivative is given by

$$\begin{array}{c}\hfill {}_0{}^{c}D_t^{\varphi }f\left(t\right){\displaystyle \frac{1}{\mathrm{\Gamma }(n-\psi )}}{\int }_a^t{\displaystyle \frac{f^n\left(\xi \right)}{{(t-\xi )}^{\psi +1-n}}}d\xi ,n-1<\psi ,n\in \mathbb{N}.\end{array}$$

(4)

Definition 2 

(Caputo derivative of a constant [34]). The fractional derivative for a constant function $f\left(t\right)=c$ is zero, that is,

$$\begin{array}{c}\hfill {}_0{}^{c}D_t^{\varphi }c=0,\varphi \in (0,1)\end{array}$$

(5)

Definition 3. 

The Caputo-fractional derivative in the Liouvile-Caputo sense is defined as

$$\begin{array}{c}\hfill {}_0{}^{c}D_t^{\varphi }f\left(t\right)={\displaystyle \frac{M\left(\varphi \right)}{1-\varphi }}{\int }_0^t{f}^{\prime }\left(\tau \right)exp\left[{\displaystyle \frac{\varphi (t-\tau )}{1-\tau }}\right]d\tau \end{array}$$

(6)

whereby $M\left(\varphi \right)$ is a normalized function, such as $M\left(0\right)=M\left(1\right)=1$.

Definition 4. 

The Caputo-fractional integral of order $q(0<\varphi \le 1)$ of the function $f\left(t\right)$ is defined as

$$\begin{array}{c}\hfill {}_0{}^{c}I_t^{\varphi }f\left(t\right)={\displaystyle \frac{2(1-\varphi )}{(2-\varphi )M\left(\varphi \right)}}f\left(t\right)+{\displaystyle \frac{2\varphi }{(2-\varphi )M\left(\varphi \right)}}{\int }_0^{t}f\left(s\right)ds,t\ge 0\end{array}$$

(7)

whereby $M\left(\psi \right)={\displaystyle \frac{2}{2-\varphi }},0<\varphi \le 1$.

3. Basic Properties

Non-Negativity and Boundness of Model Solutions

Theorem 1. 

For model system (3), there exists a unique solution in $(0,\infty )$; however, the solution is always positive for all values of $t\ge 0$ and remains in ${\mathcal{R}}_{+}^2$.

Proof. 

From model system (3), we first show that ${\mathcal{R}}_{+}^2=${($N_h,N_v$)$\in {\mathcal{R}}_{+}^2:N_h\ge 0,N_v\ge 0\}$ is a positive invariant. Now, we have to demonstrate that each hyperplane bounding the positive orthant and the vector field points to ${\mathcal{R}}_{+}^2.$ Let us consider the following cases:

• Case 1: Let us assume that there exists a $t_{*}>t_0$ such that $N_h\left(t_{*}\right)=0,$ and $N_h\left(t\right)<0$ for $t\in (t_{*},t_1),$ where $t_1$ is sufficiently close to $t_{*},$ if $N_h\left(t_{*}\right)=0$; then, we see that
  ${\displaystyle {}_0{}^{c}D_t^{\varphi }{N}_h\left(t_{*}\right)-{\mathrm{\Lambda }}_h>0}$. This implies that ${\displaystyle {}_0{}^{c}D_t^{\varphi }{N}_h\left(t\right)>0}$ for all $t\in [t_{*},t_1].$
• Case 2: Let us assume that there exists a $t_{*}>t_0$ such that $N_v\left(t_{*}\right)=0,$ and $N_v\left(t\right)<0$ for $t\in (t_{*},t_1),$ where $t_1$ is sufficiently close to $t_{*},$ if $N_v\left(t_{*}\right)=0$; then, we see that
  ${\displaystyle {}_0{}^{c}D_t^{\varphi }{N}_v\left(t_{*}\right)-{\mathrm{\Lambda }}_v>0}$. This implies that ${\displaystyle {}_0{}^{c}D_t^{\varphi }{N}_v\left(t\right)>0}$ for all $t\in [t_{*},t_1].$
  The above discussion shows that the three hyperplanes bound the orthants, meaning the vector field points to ${\mathcal{R}}_{+}^2.$ This shows that all the solutions of model system (3) remain positive for all $t\ge 0$. □

Theorem 2. 

Let $\mathrm{\Phi }\left(t\right)=(N_h\left(t\right),N_v\left(t\right))$ be the unique solution of model system (3) for all $t\ge 0;$ then, the solution $\mathrm{\Phi }\left(t\right)$ is bounded above, that is, $\mathrm{\Phi }\left(t\right)\in \mathrm{\Omega }$, where Ω is the feasible region defined as

$$\mathrm{\Omega }=\left\{\begin{array}{c}\left(\begin{array}{c}{N}_h\left(t\right)\\ N_V\left(t\right)\end{array}\right)\in {\mathbb{R}}_{+}^2\left|\begin{array}{c}{\displaystyle 0\le N_h\left(t\right)\le C_h,}\hfill \\ {\displaystyle 0\le N_v\left(t\right)\le C_v.}\hfill \end{array}\right.\hfill \end{array}\right\}$$

the interior of which is denoted by int$\left(\mathrm{\Omega }\right)$ and given by

$$int\left(\mathrm{\Omega }\right)=\left\{\begin{array}{c}\left(\begin{array}{c}{N}_h\left(t\right)\\ N_v\left(t\right)\end{array}\right)\in {\mathbb{R}}_{+}^2\left|\begin{array}{c}{\displaystyle 0<N_h\left(t\right)<C_h,}\hfill \\ {\displaystyle 0<N_v\left(t\right)<C_v.}\hfill \end{array}\right.\hfill \end{array}\right\}$$

Proof. 

Here, we prove that the solutions of model system (3) are bounded for all $t\ge 0$. Biologically, the lowest possible value of each of the states of model system (3) is zero. Next, we determine the upper bound of the states. Based on this discussion, it is easy to show that the following conditions hold for the biological relevance of a species. $0\le N_h\left(t\right)\le C_h,$ and $0\le N_v\left(t\right)\le C_r.$ From these conditions, we have

$$\begin{array}{c}\hfill {\displaystyle {}_0{}^{c}D_t^{\varphi }{N}_h\le {\mathrm{\Lambda }}_h-{\mu }_h{N}_h\left(t\right)}\end{array}$$

By using the Laplace transformation conditions, we have

$$\begin{array}{c}\hfill S^{\varphi }L\left[N_h\left(t\right)\right]-S^{\varphi -1}{N}_h\left(0\right)\le {\displaystyle \frac{{\mathrm{\Lambda }}_h^{\varphi }}{S}}-{\mu }_h^{\varphi }L\left[N_h\left(t\right)\right]\end{array}$$

By collecting the likely terms, we have

$$\begin{array}{c}\hfill L\left[N_h\left(t\right)\right]\le {\mathrm{\Lambda }}_h^{\varphi }{\displaystyle \frac{S^{-1}}{S^{\varphi }+{\mu }_h^{\varphi }}}+N_h\left(0\right){\displaystyle \frac{S^{\varphi -1}}{S^{\psi }+{\mu }_h^{\psi }}}.\end{array}$$

$$\begin{array}{c}\hfill ={\mathrm{\Lambda }}_h^{\varphi }{\displaystyle \frac{S^{\varphi -(1+\psi )}}{S^{\varphi }+{\mu }_h^{\varphi }}}+N_h\left(0\right){\displaystyle \frac{S^{\varphi -1}}{S^{\varphi }+{\mu }_h^{\varphi }}}\end{array}$$

By using the inverse Laplace transform, we have

$$\begin{array}{c}{N}_h\left(t\right)\le L^{-1}\left\{{\mathrm{\Lambda }}_h^{\varphi }{\displaystyle \frac{S^{\varphi -(1+\varphi )}}{S^{\varphi }+{\mu }_h^{\varphi }}}\right\}-N_h\left(0\right)L^{-1}\left\{{\displaystyle \frac{S^{\varphi -1}}{S^{\varphi }+{\mu }_h^{\varphi }}}\right\}\hfill \\ \le {\mathrm{\Lambda }}_h^{\varphi }{t}^{\varphi }{E}_{\varphi ,\varphi +1}(-{\mu }_h^{\varphi })t^{\varphi }+N_h\left(0\right)E_{\varphi ,1}(-{\mu }_h^{\varphi })t^{\varphi }\hfill \\ \le {\displaystyle \frac{{\mathrm{\Lambda }}_h^{\varphi }}{{\mu }_h^{\varphi }}}{t}^{\varphi }{E}_{\varphi ,\varphi +1}(-{\mu }_h^{\varphi })t^{\varphi }+N_h\left(0\right)E_{\varphi ,1}(-{\mu }_h^{\varphi })t^{\varphi }\hfill \\ \le Max\left\{{\displaystyle \frac{{\mathrm{\Lambda }}_h^{\varphi }}{{\mu }_h^{\varphi }}},N_h\left(0\right)\right\}\left(t^{\varphi }{E}_{\varphi ,\varphi +1}(-{\mu }_h^{\varphi })t^{\varphi }+E_{\varphi ,1}(-{\mu }_h^{\varphi })t^{\varphi }\right)\hfill \\ \hfill ={\displaystyle \frac{C}{\mathrm{\Gamma }\left(1\right)}}=C_h.\hfill \end{array}$$

where $C_h=Max\left\{{\displaystyle \frac{{\mathrm{\Lambda }}_h^{\varphi }}{{\mu }_h^{\varphi }}},N_h\left(0\right)\right\}$. Therefore, $N_h\left(t\right)$ is bounded above.

From the vector population, we have

$$\begin{array}{c}\hfill {\displaystyle {}_0{}^{c}D_t^{\varphi }{N}_v\le {\mathrm{\Lambda }}_v^{\varphi }-{\mu }_v^{\varphi }{N}_v\left(t\right)}\end{array}$$

By using the Laplace transformation conditions, we have

$$\begin{array}{c}\hfill S^{\psi }L\left[N_v\left(t\right)\right]-S^{\varphi -1}{N}_v\left(0\right)\le {\displaystyle \frac{{\mathrm{\Lambda }}_v^{\varphi }}{S}}-{\mu }_v^{\varphi }L\left[N_v\left(t\right)\right]\end{array}$$

By collecting the likely terms, we have

$$\begin{array}{c}\hfill L\left[N_v\left(t\right)\right]\le {\mathrm{\Lambda }}_v^{\varphi }{\displaystyle \frac{S^{-1}}{S^{\varphi }+{\mu }_v^{\varphi }}}+N_v\left(0\right){\displaystyle \frac{S^{\varphi -1}}{S^{\varphi }+{\mu }_v^{\varphi }}}.\end{array}$$

$$\begin{array}{c}\hfill ={\mathrm{\Lambda }}_v^{\varphi }{\displaystyle \frac{S^{\varphi -(1+\varphi )}}{S^{\varphi }+{\mu }_v^{\varphi }}}+N_v\left(0\right){\displaystyle \frac{S^{\varphi -1}}{S^{\varphi }+{\mu }_v^{\varphi }}}\end{array}$$

By using the inverse Laplace transform, we have

$$\begin{array}{c}{N}_v\left(t\right)\le L^{-1}\left\{{\mathrm{\Lambda }}_v^{\varphi }{\displaystyle \frac{S^{\varphi -(1+\varphi )}}{S^{\varphi }+{\mu }_v^{\varphi }}}\right\}-N_v\left(0\right)L^{-1}\left\{{\displaystyle \frac{S^{\varphi -1}}{S^{\varphi }+{\mu }_v^{\varphi }}}\right\}\hfill \\ \le {\mathrm{\Lambda }}_v^{\varphi }{t}^{\varphi }{E}_{\varphi ,\varphi +1}(-{\mu }_v^{\varphi })t^{\varphi }+N_v\left(0\right)E_{\varphi ,1}(-{\mu }_v^{\varphi })t^{\varphi }\hfill \\ \le {\displaystyle \frac{{\mathrm{\Lambda }}_v^{\varphi }}{{\mu }_v^{\varphi }}}{t}^{\varphi }{E}_{\varphi ,\varphi +1}(-{\mu }_v^{\varphi })t^{\varphi }+N_v\left(0\right)E_{\varphi ,1}(-{\mu }_v^{\varphi })t^{\varphi }\hfill \\ \le Max\left\{{\displaystyle \frac{{\mathrm{\Lambda }}_v^{\varphi }}{{\mu }_v^{\varphi }}},N_v\left(0\right)\right\}\left(t^{\varphi }{E}_{\varphi ,\varphi +1}(-{\mu }_v^{\varphi })t^{\varphi }+E_{\varphi ,1}(-{\mu }_v^{\varphi })t^{\varphi }\right)\hfill \\ \hfill ={\displaystyle \frac{C}{\mathrm{\Gamma }\left(1\right)}}=C_v.\hfill \end{array}$$

where $C_v=Max\left\{{\displaystyle \frac{({\mathrm{\Lambda }}_v^{\varphi }}{{\mu }_v^{\varphi }}},N_v\left(0\right)\right\}$. Therefore $N_v\left(t\right)$ is bounded above, and this completes the proof. □

4. Basic Reproduction Number and Existence of Equilibria

In this section, we use the next-generation method, as presented in [35], to compute the threshold quantity ${\mathcal{R}}_0$, which determines the persistence and extinction of disease in the population. It is believed that when the basic reproduction number ${\mathcal{R}}_0>1$, the disease persists in the population, and it dies out when ${\mathcal{R}}_0<1$. The model system (3) always has a disease-free equilibrium, ${\mathcal{E}}^0$, given by

$$\begin{array}{c}\hfill {\mathcal{E}}^0:\left(S_h^0,E_h^0,I_h^0,R_h^0,S_v^0,E_v^0,I_v^0\right)=\left({\displaystyle \frac{{\mathrm{\Lambda }}_h}{{\mu }_h}},0,0,0,{\displaystyle \frac{{\mathrm{\Lambda }}_v}{{\mu }_v+ϵ}},0,0\right).\end{array}$$

By following the next generation matrix approach as used in [35], the non-negative matrix $\mathcal{F}$ that denotes the generation of new infection and the nonsingular matrix $\mathcal{V}$ that denotes the disease transfer among compartments evaluated at ${\mathcal{E}}^0$ are defined as follows:

$$\mathcal{F}=\left(\begin{array}{cccc}0& 0& 0& (1-\omega ){\delta }_v{\beta }_h{S}_h^{*}\\ 0& 0& 0& 0\\ 0& (1-\omega ){\delta }_v{\beta }_v{S}_v^{*}& 0& 0\\ 0& 0& 0& 0\end{array}\right)$$

(8)

$$\mathcal{V}=\left(\begin{array}{cccc}{\alpha }_h+{\mu }_h+{\theta }_h& 0& {\displaystyle 0}& 0\\ -{\alpha }_h& {\mu }_h+{\gamma }_h+d_h& 0& 0\\ 0& 0& {\mu }_v+{\alpha }_v+ϵ& 0\\ 0& 0& -{\alpha }_v& {\mu }_v+ϵ\end{array}\right)$$

(9)

By using (8) and (9), the basic reproduction number ${\mathcal{R}}_0$ of model system (3) is given by

$${\mathcal{R}}_0={\displaystyle \frac{(1-\omega ){\delta }_v}{{\mu }_v+ϵ}}\sqrt{{\displaystyle \frac{{\beta }_v{\beta }_h{\mathrm{\Lambda }}_v{\mathrm{\Lambda }}_h}{({\alpha }_h+{\mu }_h+{\theta }_h)({\mu }_h+{\gamma }_h+d_h)({\mu }_v+ϵ+{\alpha }_v)}}}$$

The basic reproduction number, ${\mathcal{R}}_0$, is defined as the expected number of secondary cases (mosquito or human) produced in a completely susceptible population by one infectious individual (mosquito or human, respectively) during its lifetime of being infectious. The square root here is due to the fact that the generation of secondary cases in malaria diseases requires two transmission processes.

Theorem 3. 

If ${\mathcal{R}}_0<1,$ then the DFE of system (3) is globally asymptotically stable in Ω; otherwise, it is unstable.

Proof. 

By considering only the infected compartments from (3), one can write in the following ways:

$$\begin{array}{c}\hfill {}_{t_0}{}^{c}D_t^{\varphi }x=(F-V)x,\end{array}$$

where $x={(E_h,A_h,I_h,E_v,I_v)}^T$, with F and V defined as follows:

$$\mathcal{F}=\left(\begin{array}{cccc}0& 0& 0& (1-\omega ){\delta }_v{\beta }_h{S}_h^{*}\\ 0& 0& 0& 0& 0\\ 0& (1-\omega ){\delta }_v{\beta }_v{S}_v^{*}& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right),$$

(10)

$$\mathcal{V}=\left(\begin{array}{cccc}{\alpha }_h+{\mu }_h+{\theta }_h& 0& {\displaystyle 0}& 0\\ -{\alpha }_h& {\mu }_h+{\gamma }_h+d_h& 0& 0\\ 0& 0& {\mu }_v+{\alpha }_v+ϵ& 0\\ 0& 0& -{\alpha }_v& {\mu }_v+{ϵ}_v\end{array}\right)$$

(11)

We can observe that by using direct calculation, we can verify that $V^{-1}F$ is a non-negative and irreducible matrix, and $\rho \left(V^{-1}F\right)={\mathcal{R}}_0.$ It follows from the Perron-Frobenius theorem [36] that $V^{-1}F$ has a positive left eigenvector $\mathbf{w}$ associated with ${\mathcal{R}}_0,$ that is,

$$\begin{array}{c}\hfill \mathbf{w}{V}^{-1}F={\mathcal{R}}_0\mathbf{w}.\end{array}$$

Since $\mathbf{w}{V}^{-1}$ is a positive vector, we propose the following Lyapunov function to study the global stability of disease-free equilibrium:

$$\begin{array}{c}\hfill \mathcal{L}\left(t\right)=\mathbf{w}{V}^{-1}x.\end{array}$$

Differentiating $\mathcal{L}$ along solutions of (3) leads to

$$\begin{array}{ccc}\hfill {\displaystyle {}_{t_0}{}^{c}D_t^{\varphi }\mathcal{L}\left(t\right)}& =& \mathbf{w}{V}^{-1}{}_{t_0}{}^{c}D_t^{\varphi }x\le \mathbf{w}{V}^{-1}(F-V)x\hfill \\ & =& ({\mathcal{R}}_0-1)\mathbf{w}x\le 0if{\mathcal{R}}_0\le 1.\hfill \end{array}$$

It can be easily verified that the largest invariant subset of $\mathrm{\Gamma }$, where ${}_{t_0}{}^{c}D_t^{\varphi }\mathcal{L}\left(t\right)=0$, is the singleton $\left\{{\mathcal{E}}^0\right\}$. Therefore, according to LaSalle’s invariance principle [37], ${\mathcal{E}}^0$ is globally asymptotically stable in $\mathrm{\Omega }$ when ${\mathcal{R}}_0\le 1$. □

4.1. Local Stability of the Disease-Free Equilibrium Point

In this section, we analyze the local stability of the disease-free equilibrium (DFE) point using the Jacobian matrix and the Routh-Hurwitz criterion.

Theorem 4 

(Routh-Hurwitz Criterion). Given the polynomial $P\left(\lambda \right)={\lambda }^n+a_1{\lambda }^{n-1}+...+a_{n-1}\lambda +a_n$, where the coefficients $a_i$ are the real constants, $i=1,...,n$ define the n Hurwitz matrices using the coefficients $a_i$ of the characteristic polynomial: $H_1=\left(a_1\right),$

$$H_1=\left(\begin{array}{cc}{a}_1& 1\\ a_3& a_2\end{array}\right)$$

(12)

and

$$H_2=\left(\begin{array}{ccccccc}{a}_1& 1& {\displaystyle 0}& .& .& .& 0\\ a_3& a_2& a_1& .& .& .& 0\\ a_5& a_4& a_3& .& .& .& 0\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ .& .& .& .& .& .& .\\ 0& 0& 0& 0& 0& 0& a_n\end{array}\right)$$

(13)

where $a_j=0$ if $j>n.$ All the roots of the polynomial $p\left(\lambda \right)$ are negative or have a negative real part if and only if the determinants of all Hurwitz matrices are positive: $det{H}_j>0,j=1,2,3,...,n.$ When $n=4,$ the Routh-Hurwitz criterion simplifies to $H_1>0,H_2>0,H_3>0$ and

$$H_4=\left(\begin{array}{cccc}{a}_1& 1& 0& 0\\ a_3& a_2& 1& 0\\ a_5& a_4& a_3& 1\\ a_7& a_6& a_5& a_4\end{array}\right)$$

(14)

which gives the conditions that $n=4:$ $a_1>0,a_3>0,a_4>0$ and $a_1{a}_2{a}_3>a_3^2+a_1^2{a}_4.$

Theorem 5. 

The DFE is locally asymptomatically stable if ${\mathcal{R}}_0<1$, and it is unstable if ${\mathcal{R}}_0>1$.

Proof. 

The Jacobian matrix of the system is given by

$$H_4=\left(\begin{array}{ccccccc}{m}_1& 0& 0& {\rho }_h& 0& 0& m_2\\ m_3& m_4& 0& 0& 0& 0& m_5\\ 0& {\alpha }_h& m_6& 0& 0& 0& 0\\ 0& {\theta }_h& {\gamma }_h& m_7& 0& 0& 0\\ 0& {\alpha }_h& m_8& 0& m_9& 0& 0\\ 0& 0& m_{10}& 0& & m_{12}& 0\\ 0& 0& 0& 0& 0& {\alpha }_v& -({\mu }_v+ϵ)\end{array}\right)$$

(15)

where $m_1=-(1-\omega ){\delta }_v{\beta }_h{I}_v-{\mu }_h,m_2=-(1-\omega ){\delta }_v{\beta }_h{S}_h,m_3=(1-\omega ){\delta }_v{\beta }_h{I}_v,m_4=-({\mu }_h+d_h+{\gamma }_h)$, $m_5=(1-\omega ){\delta }_v{\beta }_h{S}_h,m_6=-({\mu }_h+d_h+{\gamma }_h),m_7=-({\mu }_h+{\rho }_h),m_8=-(1-\omega ){\delta }_v{\beta }_v{S}_v,m_9=-(1-\omega ){\delta }_v{\beta }_v{I}_h-({\mu }_v+ϵ),m_{10}=(1-\omega ){\delta }_v{\beta }_v{S}_v,m_{11}=(1-\omega ){\delta }_v{\beta }_v{I}_h,m_{12}=-({\mu }_v+{\alpha }_v+ϵ)$.

At the DFE, the Jacobian matrix becomes

$$J=\left(\begin{array}{ccccccc}-{\mu }_h& 0& 0& {\rho }_h& 0& 0& -n_1\\ 0& -n_2& 0& 0& 0& 0& n_1\\ 0& {\alpha }_h& -n_3& 0& 0& 0& 0\\ 0& {\theta }_h& {\gamma }_h& -({\mu }_h+{\rho }_h)& 0& 0& 0\\ 0& 0& -n_4& 0& -({\mu }_v+ϵ)& 0& 0\\ 0& 0& n_4& 0& 0& -n_5& 0\\ 0& 0& 0& 0& 0& {\alpha }_v& -({\mu }_v+ϵ)\end{array}\right)$$

(16)

where $n_1={\displaystyle \frac{(1-\omega ){\delta }_v{\beta }_h{\mathrm{\Lambda }}_h}{{\mu }_h}},n_2=({\mu }_h+{\alpha }_h+{\theta }_h),n_3=({\mu }_h+d_h+{\gamma }_h)$ $n_4={\displaystyle \frac{(1-\omega ){\delta }_v{\beta }_v{\mathrm{\Lambda }}_v}{{\mu }_v+ϵ}}$, $n_5=({\mu }_v+{\alpha }_v+ϵ)$.

The eigenvalues of $J_{DFE}-\mathrm{\Lambda }I$ are given by solving

$$det\left(\begin{array}{ccccccc}-{\mu }_h-\lambda & 0& 0& {\rho }_h& 0& 0& -n_1\\ 0& -n_2-\lambda & 0& 0& 0& 0& n_1\\ 0& {\alpha }_h& -n_3-\lambda & 0& 0& 0& 0\\ 0& {\theta }_h& {\gamma }_h& -({\mu }_h+{\rho }_h)-\lambda & 0& 0& 0\\ 0& 0& -n_4& 0& -({\mu }_v+ϵ)-\lambda & 0& 0\\ 0& 0& n_4& 0& 0& -n_5-\lambda & 0\\ 0& 0& 0& 0& 0& {\alpha }_v& -({\mu }_v+ϵ)-\lambda \end{array}\right)$$

(17)

To obtain ${\lambda }_1=-{\mu }_h,{\lambda }_2=-{\mu }_v-ϵ,{\lambda }_3=-{\rho }_h-{\mu }_h$, which are all negative. The other four eigenvalues are obtained by solving the polynomial:

$$\begin{array}{ccc}\hfill {\lambda }^4& +& \left(2({\mu }_h+{\mu }_v+ϵ)+{\alpha }_h+{\alpha }_v+d_h+{\theta }_h+{\gamma }_h\right){\lambda }^3\hfill \\ & & +\left(({\mu }_h+{\alpha }_h+{\theta }_h)({\mu }_h+d_h+{\gamma }_h)+\left(2({\mu }_v+ϵ)+{\alpha }_v\right)+({\mu }_v+{\alpha }_v+ϵ)({\mu }_v+ϵ)\right){\lambda }^2\hfill \\ & & +(({\mu }_h+{\alpha }_h+{\theta }_h)({\mu }_h+d_h+{\gamma }_h)+\left(2({\mu }_v+ϵ)+{\alpha }_v\right)\hfill \\ & & +(2{\mu }_h+{\alpha }_h+d_h+{\theta }_h+{\gamma }_h)({\mu }_v+{\alpha }_v+ϵ)({\mu }_v+ϵ))\lambda \hfill \\ & & +\left(({\mu }_h+{\alpha }_h+{\theta }_h)({\mu }_h+d_h+{\gamma }_h)({\mu }_v+{\alpha }_v+ϵ)({\mu }_v+ϵ)\left(1-{\mathcal{R}}_0^2\right)\right)=0\hfill \end{array}$$

For the eigenvalues to be negative, we use the Routh-Hurwitz criterion; then, we have $2({\mu }_v+{\mu }_h+ϵ)+{\alpha }_h+{\alpha }_v+d_h+{\theta }_h{\gamma }_h>0,$  $\left(({\mu }_h+{\alpha }_h+{\theta }_h)({\mu }_h+d_h+{\gamma }_h)+\left(2({\mu }_v+ϵ)+{\alpha }_v\right)+({\mu }_v+{\alpha }_v+ϵ)({\mu }_v+ϵ)\right)>0,$  $\left(({\mu }_h+{\alpha }_h+{\theta }_h)({\mu }_h+d_h+{\gamma }_h)+\left(2({\mu }_v+ϵ)+{\alpha }_v\right)+(2{\mu }_h+{\alpha }_h+d_h+{\theta }_h+{\gamma }_h)({\mu }_v+{\alpha }_v+ϵ)({\mu }_v+ϵ)\right)>0,$  $\left(({\mu }_h+{\alpha }_h+{\theta }_h)({\mu }_h+d_h+{\gamma }_h)({\mu }_v+{\alpha }_v+ϵ)({\mu }_v+ϵ)\left(1-{\mathcal{R}}_0^2\right)\right)>0$ and $a_1{a}_2{a}_3>a_3^2+a_1^2{a}_4.$ Hence, this implies that ${\mathcal{R}}_0<1$, and the DFE is locally asymptotically stable. □

4.2. Endemic Equilibrium Point

The endemic equilibrium point $(S_h^{*},E_h^{*},I_h^{*},R_h^{*},S_v^{*},E_v^{*},I_v^{*})$ of model system (3) is obtained by simultaneously solving the below system (18) and is implicitly expressed in terms of $I_h^{*}$:

$$\left\{\begin{array}{ccc}{\displaystyle {\mathrm{\Lambda }}_h+{\rho }_h{R}_h^{*}-(1-\omega ){\delta }_v{\beta }_h{S}_h^{*}{I}_v^{*}-{\mu }_h{S}_h^{*}}\hfill & =\hfill & 0,\hfill \\ {\displaystyle (1-\omega ){\delta }_v{\beta }_h{S}_h^{*}{I}_v^{*}-({\mu }_h+{\alpha }_h+{\theta }_h)E_h^{*}}\hfill & =\hfill & 0,\hfill \\ {\displaystyle {\alpha }_h{E}_h^{*}-({\mu }_h+d_h+{\gamma }_h)I_h^{*}}\hfill & =\hfill & 0,\hfill \\ {\displaystyle {\theta }_h{E}_h^{*}+{\gamma }_h{I}_h^{*}-({\rho }_h+{\mu }_h)R_h^{*}}\hfill & =\hfill & 0,\hfill \\ {\displaystyle {\mathrm{\Lambda }}_v-(1-\omega ){\delta }_v{\beta }_v{I}_h^{*}{S}_v^{*}-({\mu }_v+ϵ)S_v^{*}}\hfill & =\hfill & 0,\hfill \\ {\displaystyle (1-\omega ){\delta }_v{\beta }_v{I}_h^{*}{S}_v^{*}-({\mu }_v+{\alpha }_v+ϵ)E_v^{*}}\hfill & =\hfill & 0,\hfill \\ {\displaystyle {\alpha }_v{E}_v^{*}-({\mu }_v+ϵ)I_v^{*}}\hfill & =\hfill & 0.\hfill \end{array}\right.$$

(18)

From the second equation of system (18), we have the following:

$$E_h^{*}={\displaystyle \frac{(1-\omega ){\delta }_v{\beta }_h{S}_h^{*}{I}_v^{*}}{{\mu }_h+{\alpha }_h+\theta -h}}$$

(19)

By substituting (19) into Equation (4) of system (18) and solving for $R_h^{*}$, we obtain

$$R_h^{*}={\displaystyle \frac{(1-\omega ){\delta }_v{\theta }_h{\S }_h^{*}{I}_v^{*}+{\gamma }_h({\mu }_h+{\alpha }_h+{\theta }_h)I_h^{*}}{({\mu }_h+{\alpha }_h+{\theta }_h)({\rho }_h+{\mu }_h)}}$$

(20)

By substituting Equation (20) in the first equation of system (18) and solving for $S_h^{*}$, we obtain the following:

$$S_h^{*}={\displaystyle \frac{\left({\mathrm{\Lambda }}_h({\rho }_h+{\mu }_h)+{\gamma }_h{\rho }_h{I}_h^{*}\right)({\mu }_{h}h+{\alpha }_h+{\theta }_h)}{\left((1-\omega ){\delta }_v{\beta }_h{I}_v^{*}+{\mu }_h\right)({\mu }_h+{\alpha }_{h}h+{\theta }_h)({\rho }_h+{\mu }_h)-(1-\omega ){\delta }_v{\beta }_h{\theta }_h{\rho }_h{I}_v^{*}}}$$

(21)

By substituting (21) with (19), the simplified expression for (19) becomes

$$E_h^{*}={\displaystyle \frac{(1-\omega ){\delta }_v{\beta }_v{I}_v^{*}\left({\mathrm{\Lambda }}_h({\rho }_h+{\mu }_h)+{\gamma }_h{\rho }_h{I}_h^{*}\right)}{\left((1-\omega ){\delta }_v{\beta }_h{I}_v^{*}+{\mu }_h\right)({\mu }_h+{\alpha }_h+{\theta }_h)({\rho }_h+{\mu }_h)-(1-\omega ){\delta }_v{\beta }_h{\theta }_h{\rho }_h{I}_v^{*}}}$$

(22)

Similarly, for the vector population from Equation (6) of (18), we have

$$E_v^{*}={\displaystyle \frac{(1-\omega ){\delta }_v{\beta }_v{S}_v^{*}{I}_h^{*}}{{\mu }_v+{\alpha }_v+ϵ}}$$

(23)

From Equation (5) of system (18), we obtain the following:

$$S_v^{*}={\displaystyle \frac{{\mathrm{\Lambda }}_v}{(1-\omega ){\delta }_v{\beta }_v{I}_h^{*}+{\mu }_v+ϵ}}$$

(24)

By substituting (24) with (23), we obtain the simplified expression of $E_v^{*}$ as

$$E_v^{*}={\displaystyle \frac{(1-\omega ){\mathrm{\Lambda }}_v{\delta }_v{\beta }_v{I}_h^{*}}{({\mu }_v+{\alpha }_v+ϵ)\left((1-\omega ){\delta }_v{\beta }_v{I}_h^{*}+{\mu }_v+ϵ\right)}}$$

(25)

where, from Equation (7) of system (18) and Equation (25), we have

$$I_v^{*}={\displaystyle \frac{(1-\omega ){\mathrm{\Lambda }}_v{\alpha }_v{\delta }_v{\beta }_v{I}_h^{*}}{({\mu }_v+{\alpha }_v+ϵ)({\mu }_v+ϵ)\left((1-\omega ){\delta }_v{\beta }_v{I}_h^{*}+{\mu }_v+ϵ\right)}}$$

(26)

Therefore, Equations (20)–(23), (25), and (26) give the endemic equilibrium point (EEP) in terms of $I_h^{*}$.

5. Results and Discussion

In this section, we present and analyze the numerical solution of the model by utilizing the Euler and Adam-Bashforth-Moulton scheme. The numerical results provide valuable insights into the behavior of the mathematical model system and the effects of its fractional order. Initially, we selected appropriate parameter values along with the initial conditions for the classes: $S_h\left(0\right)=500000,$ $E_h\left(0\right)=1000,$ $I_h\left(0\right)=100,$ $R_h\left(0\right)=20,$ $S_v\left(0\right)=1500,$ $E_v\left(0\right)=10,$ and $I_v\left(0\right)=10.$ Some parameters in this model have been fitted with the real data of the reported cases to align with the realism of the proposed model, as shown in

Table 1. Definition of model parameters and values.

Symbol	Definition	Value	Units	Source
${\beta }_h$	disease transmission from mosquito to human	0.001	${\mathrm{day}}^{-1}$	[9,38]
${\beta }_v$	disease transmission from human to mosquito	0.0001	${\mathrm{day}}^{-1}$	[9,38]
${\mu }_h$	natural mortality rate of human	$1/(60\ast 365)$	${\mathrm{day}}^{-1}$	[9,38]
${\mu }_v$	natural mortality rate of vector	$1/15$	${\mathrm{day}}^{-1}$	[9,38]
${\alpha }_h$	progression rate of human, from incubation to infectious	1/17	${\mathrm{day}}^{-1}$	[9,39]
${\alpha }_v$	progression rate of vector, from incubation to infectious	1/18	${\mathrm{day}}^{-1}$	[9,39]
${\gamma }_h$	progression rate of human, from infectious to recovered class	$0.012$	${\mathrm{day}}^{-1}$	[39]
${\mathrm{\Lambda }}_h$	new recruitment of human	10	${\mathrm{day}}^{-1}$	[9,38]
${\mathrm{\Lambda }}_v$	new recruitment of Aedes mosquito	50	${\mathrm{day}}^{-1}$	[39]
${\theta }_h$	progression rate of exposed human to recovered class	$[0, 1]$	${\mathrm{day}}^{-1}$	fitted
$ϵ$	Rate of use of insecticides	$[0, 1]$	${\mathrm{day}}^{-1}$	fitted
$\omega$	Proportion of human progress to infectious class	$[0, 1]$	${\mathrm{day}}^{-1}$	fitted
${\delta }_v$	Rate of mosquito biting a human	3	${\mathrm{day}}^{-1}$	[39]

.

5.1. Model Parameter Estimations

In this section, we perform model fitting using the data of malaria disease for 12 years, as reported in Zimbabwe, to fit the proposed model (3). Some of the parameters used in the simulations were adopted from the literature, as shown in Table 1, and some parameters were estimated by using the root mean square error (RMSE) in the following formula:

$$\begin{array}{c}\hfill RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{k=1}^{12}{(I\left(k\right)-\widehat{I}\left(k\right))}^2},\end{array}$$

(27)

where n is the number of yearly reported malaria cases for 12 years. We assumed the initial population to be as follows: $S_h\left(0\right)=500000,$ $E_h\left(0\right)=1000,$ $I_h\left(0\right)=100,$ $R_h\left(0\right)=20,$ $S_v\left(0\right)=1500,$ $E_v\left(0\right)=10,$ and $I_v\left(0\right)=10.$ In addition, from model (3), the generated new cases are obtained using the term $(1-\omega ){\delta }_v({\beta }_h{S}_h{I}_v+{\beta }_v{S}_v{I}_h)$, which counts the detected cases.

Figure 2 presents the reported malaria cases in Zimbabwe over a 12-year period, showing unpredictable fluctuations, with reductions and increases occurring inconsistently. Figure 3 illustrates that we simultaneously fitted the malaria case data using both the classical integer-order and fractional-order models. The results demonstrate that the fractional-order model aligns more closely with the reported malaria cases compared to the classical integer-order model.

5.2. Sensitivity Analysis

In this section, we used the partial rank correlation coefficient (PRCC) to perform a global sensitivity analysis of model (3) to identify the most significant parameters that influence the spread of disease in the population. Figure 4 shows the PRCC results of the ${\mathcal{R}}_0$ related to the number of new cases generated in the population reveal seven key influential parameters. These include the new recruitment of humans and vectors (${\mathrm{\Lambda }}_h$ and ${\mathrm{\Lambda }}_v$), the force of infections (${\beta }_h$ and ${\beta }_v$), vector biting rate (${\delta }_v$), and incubation periods (${\alpha }_h$ and ${\alpha }_v$), both of which are positively correlated with ${\mathcal{R}}_0$. On the other hand, parameters such as health education ($\omega$), use of insecticides ($ϵ$), and recovery rate of infected individuals (${\gamma }_h$) have a negative correlation with ${\mathcal{R}}_0$, meaning that whenever these parameters increase, the magnitude of ${\mathcal{R}}_0$ decreases and, hence, the disease dies in the population.

Figure 5 and Figure 6 shows the global sensitivity analysis of model (3) on $I_h$ regarding the key parameters that affect the dynamic of disease in the population. Overall, we noted that an increase in the parameters that represent the rate of insecticides, health education campaigns, and vector mortality rate leads to a decrease in the number of infected humans in the population. In contrast, we observed that an increase in the parameters that represent the rate of infections and vector biting rate leads to an increase in the number of infected humans. In particular, one can observe that as the rate of insecticides ($ϵ$) and vector mortality rate (${\mu }_v$) increase above $0.5$, the number of infected humans ceases in the population.

Figure 7 shows a contour plot of ${\mathcal{R}}_0$ (a) as a function of $ϵ$ (use of insecticides) and ${\theta }_h$ (natural recovery rate of exposed vectors). Overall, one can note that in the presence of insecticides, only the value of ($ϵ$) must be greater than $0.3$ to reduce ${\mathcal{R}}_0$ to less than a unit. The results provide the implication that policymakers must implement a $30\%$ reduction in the use of insecticides to eliminate malaria disease in the population. Figure 7b shows the contour plot of ${\mathcal{R}}_0$ as the function of $ϵ$ (use of insecticides) and $\omega$ (health education campaigns). We observed that when education campaigns and the use of insecticides are implemented simultaneously in the population, the rate of insecticides ($ϵ$) must be $0.05$ to reduce the ${\mathcal{R}}_0$ to less than a unit. This result demonstrates the fact that to minimize the spread of malaria in the community, policymakers must put effort into health education campaigns regarding how people can prevent contact with mosquitoes.

To investigate the role of memory effects on the spread of malaria disease, we numerically performed the simulation of system (3) for ${\mathcal{R}}_0<1$ and ${\mathcal{R}}_0>1$, as presented in Figure 8 and Figure 9, respectively. The order of derivative $\varphi$ was varied within the reasonable range and was set to ($\varphi =0.5,0.7,0.9,1$), as mentioned in the literature, where the fractional order $\varphi =1.0$ is the fractional-order model based on the Caputo derivative and becomes a classical ordinary differential model. Based on the numerical illustrations, one can note that as the order of derivatives, $\varphi$, decreased from 1, the memory effects of the system increased, and the model solution increased quickly, peaking earlier and converging with the unique equilibrium point. In particular, for ${\mathcal{R}}_0<1$, the model solution converges to the disease-free equilibrium after 20 days, and for ${\mathcal{R}}_0>1$, the model solutions converge to a unique endemic equilibrium. Furthermore, we observed that when the fractional order approaches 0, the memory effects become strong, and the model solutions converge to their respective equilibrium point earlier than when the order of the derivatives approaches 1. This result was also observed in previous studies.

5.3. Effects of Insecticide Use on the Disease Dynamics

Figure 10 shows the effect of insecticides on the spread of malaria disease; we simulated system (3) using different values of $ϵ$, and the other parameter values are presented in Table 1. From the numerical results, we observe that the use of insecticides has the potential to reduce the spread of malaria in the population. In particular, one can note that for $ϵ=0$, the disease persists in the population, and when $ϵ\ge 0.9$, the magnitude of ${\mathcal{R}}_0$ is less than a unit. Thus, the disease decreases in the population. Additionally, the results in Figure 11 demonstrate the effects of insecticides on a reduction in new cases of malaria reported in the population; one can observe that implementing insecticides at $90\%$ in the community leads to a high reduction in the reported number of cases of the disease in the population.

5.4. Effects of Health Education Campaigns on the Disease Dynamics

In Figure 12, we simulate system (3) to show the effect of health education campaigns ($\omega$) on malaria disease transmission. The effectiveness of health education campaigns can be enhanced by utilizing various media platforms, such as newspapers, television, and social media, to inform the public on how to protect themselves from mosquitoes that transmit the malaria disease. Based on this assertion, we numerically simulated system (3) at $\omega =0$, $\omega =0.5$, and $\omega =0.9$ to demonstrate the effect of health education on minimizing the spread of disease in the population (Figure 13). From the numerical illustrations, one can note that as $\omega$ increases, the number of infections generated decreases. Additionally, we observed that implementing health education campaigns alone can not reduce the magnitude of transmission potential (${\mathcal{R}}_0$) to less than a unit.

6. Concluding Remarks

In this study, we developed and analyzed a novel model for malaria transmission and assessed the impact of health education campaigns and insecticide use on disease spread within the population. We computed the disease-free equilibrium and derived the reproduction number of the model using this next-generation method. A sensitivity analysis of the basic reproduction number was conducted using partial rank correlation coefficients to examine the relationship between the reproduction number and model parameters. Our results show that parameters with negative indices decrease the reproduction number when increased, while parameters with positive indices cause an increase in the reproduction number. Furthermore, we fitted model system (3) with the real data from malaria cases reported in Zimbabwe over 12 years. We observed that the fractional-order model ($\varphi =0.5$) had a better fit compared to the classical model ($\varphi =1$). To investigate the impact of memory effects, we simulated system (3) at various orders of derivatives: $\varphi :0.5,0.7,0.9,$ and 1. Overall, we found that the order of derivatives significantly influences disease dynamics in the population. Specifically, as the order of the derivative decreases from 1, the system’s behavior converges to the unique equilibrium more rapidly, highlighting the importance of memory effects in disease modeling. Finally, we simulated the model to assess the effects of the proposed interventions (health education campaigns and insecticide use). Our numerical results indicate that implementing insecticides at a 90% convergence rate in the population can lead to disease eradication. Conversely, we observed that implementing health education campaigns alone does not achieve disease eradication. The model formulated in this study could be improved in future work by incorporating delays in certain model parameters to assess their effects on disease spread. Additionally, this model could be improved by considering migration dynamics in both human and vector populations, which may further refine our understanding of malaria transmission and control strategies.