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

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]. These mosquitoes, primarily of the species Anopheles gambiae and Anopheles funestus, play a significant role in the transmission cycle of malaria [2]. Malaria was first identified as a disease transmitted by mosquitoes in the early 20th century, and it remains one of the most widespread and dangerous infectious diseases, especially in tropical and subtropical regions [3]. Similar to other vector-borne diseases, such as dengue, Zika virus, and yellow fever, malaria is transmitted to humans when an infected mosquito bites a susceptible individual during a blood meal [4]. Mosquitoes become infected when they bite an individual carrying the malaria parasite, although the mosquitoes themselves do not suffer from the disease [5]. The incubation period for malaria in humans typically ranges from 7 to 30 days after exposure, depending on the specific Plasmodium species [6]. Common symptoms of malaria include fever, chills, headache, muscle pain, fatigue, and in severe cases, anemia, organ failure, and death [7]. There is currently no universal vaccine for malaria, although some vaccines have been developed with limited success. Malaria management relies on antimalarial medications, such as artemisinin-based combination therapies, and preventive measures, including insecticide-treated bed nets, indoor residual spraying, and the use of mosquito repellents [8].

Mathematical models have proven invaluable in epidemiology, offering critical insights into the transmission dynamics of malaria and aiding public health officials in developing strategies for controlling the spread of the disease [18,23]. Through the use of mathematical models, researchers can simulate various intervention scenarios, study the potential outcomes of control measures, and better understand the disease’s behavior in different populations and environments [29,35].

Fractional calculus has been widely utilized to model real-world problems, including heat conduction, control theory, chaotic systems, and biological processes [13,22,28,13]. Fractional-order calculus is a popular area that focuses on the application of non-integer order derivatives in dynamic systems [23]. Literature suggests that modeling physical and real-world problems with non-integer order derivatives is more accurate than using integer-order derivatives [39]. The primary advantage of utilizing fractional-order derivatives in disease dynamics is their ability to effectively capture the memory effects and hereditary properties of species present in biological systems [39]. Due to its ability to capture memory effects, the presence of non-local operators, and its high accuracy, fractional-order derivative has gained and attracted many researchers in developing more dynamic and efficient mathematical models [41]. Researchers and scientists have utilized fractional operators to formulate mathematical models for various epidemics, accounting for different scenarios and constraints [17,22,23]. In contrast to integer-order derivatives, fractional-order derivatives are capable of capturing the memory effects and hereditary properties present in biological systems [26]. Additionally, it has been concluded that the membranes of cells in living organisms possess fractional-order electrical conductance, which can be categorized within groups of fractional-order derivatives [27]. Furthermore, fractional-order derivatives often provide a better fit to real epidemiological data [24]. The most widely utilized fractional derivatives are the Caputo, Riemann-Liouville, and Atangana-Baleanu derivatives [28]. In recent years, mathematical models incorporating fractional derivatives have been developed and studied [23,36,37]. for instance in [36], the Atangana-Baleanu operator was used to investigate the impact of personal protective measures on the dynamics of malaria transmission. Although the authors did not validate their model with real data, their findings indicated that as the fractional-order derivative decreases from 0 to 1, the number of infections in both humans and mosquitoes declines. A fractional-order model for co-infection of malaria and COVID-19 was proposed and studied in [37]. The authors calculated the threshold quantity and used the fixed-point theorem to analyze the model. Their findings showed that by implementing preventive measures to mitigate the risks associated with malaria and COVID-19, disease transmission within the population can be significantly reduced. Mathematical model of malaria disease incorporating fractional-order derivatives was formulated and analyzed in [23]. The authors used the Euler and Adam-Bashforth-Moulton scheme to simulate the model system using estimated parameters. Overall,their results showed that fractional-order derivatives have more influence on the dynamics of malaria disease in the population. For more details on mathematical model of vector borne disease using fractional-order derivatives we refer the readers to [13,15,17,19,21,23,24] and the references therein. In this study, we formulate a fractional-order model of malaria disease transmission and study the effects of human awareness on the use of insecticides. Although various types of fractional derivatives are available in the literature, our model utilizes the fractional-order calculus based on Caputo derivatives. The decision to use the Caputo derivative was influenced by several advantages; for a constant function, the Caputo derivative equals zero, consistent with the result for integer-order differential equations. Furthermore, the Caputo fractional derivative permits the incorporation of local initial conditions in the model’s formulation [28]. Recently studies have shown that one of the possible way of minimizing the spread of malaria in the population is based on how humans aware on the use of insecticides [21]. Additionally, health education campaigns to human that rise awareness of people on how to use insecticides that kill Aedes mosquitoes may significantly have impact to prevent the spread of malaria in the population [33]. To the best of our knowledge very few of mathematical model for malaria transmission formulated and validated with reported cases malaria. Therefore in this study, we propose a fractional-order model of malaria that incorporates two control strategies namely; health education campaigns, and use of insecticides and assess their effects in preventing the spread of malaria disease in the population.

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

2. Model Formulation

In this paper, the fractional-order model of malaria disease dynamics based on Caputo sense have been formulated and studied. The proposed model demonstrates the interplay between mosquito and human populations. Through out the document, the subscript h and v respectively denote the mosquito and human population. The human population is sub-divided into four groups based on the infection status namely; 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)$ class. Thus, the total population of human at time t is denoted by $N_h\left(t\right)$ given as:

$$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)$$

(1)

Furthermore, the total population of mosquitoes at time t is denoted by $N_v\left(t\right)$ which is sub-divided into three groups namely; the susceptible $S_v\left(t\right)$, exposed $E_v\left(t\right)$, and infectious$I_v\left(t\right)$ class. 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 parameters and variables of the proposed model are are considered to be non-negative, and the parameters are described as follows: ${\Lambda }_h$ and ${\Lambda }_v$ denote the new recruitment rate of human and vector respectively, and they all considered to be susceptible whenever coming from either birth or immigrants; ${\mu }_h$ and ${\mu }_v$ denote the natural mortality rate of human and vector respectively; ${\alpha }_h$ and ${\alpha }_v$ denote transfer rate of human and vector respectively, from exposed to infectious class; ${\theta }_h$ represents the natural recover rate of exposed individuals; ${\gamma }_h$ denotes the recover rate of infected humans; we have assumed that both recovered individuals (from exposed and infected class) recover with temporary immunity. Thus, the parameter $rh{o}_h$ represents the warning rate of immunity for recovered humans; $d_h$ represents the disease induced death rate of infected humans; ${\beta }_h$ denotes probability of disease transmission from infected mosquito to susceptible human due to successful contact between infected vector and susceptible humans; ${\beta }_v$ denotes the probability of disease transmission from infected human to susceptible mosquitoes due to successful contact between infected person and susceptible vector; ${\delta }_v$ denotes the mosquito biting rate on the humans; Recently studies have shown that one of the possible way of minimizing the spread of malaria disease in the population is based on how community is aware on the use of prevention measures such as use of mosquito nets and insecticides that reduce the possibility of human in 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 use of insecticides represented by $ϵ$. Based on the above assumptions, we have assumed the following flow chart and non-linear differential equations:

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

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

(3)

3. Basic Properties

Herein, we study the basic properties of system $\left(\right)$ which are essential in the proof of the stability analysis.

3.1. Positivity Of Solutions

In this section, we investigate the asymptotic behavior of orbits starting in the non-negative cone ${\mathbb{R}}_7^{+}$. Obviously, system $\left(\right)$, which is a $C^{\infty }$ differential system admits a unique maximal solution for any associated Cauchy problem. In what follows we claim the following theorem:

Theorem 1.

Let $\left(t_0=0,X_0\left(t_0\right)=\left(S_h\left(0\right),E_h\left(0\right),I_h\left(0\right),R_h\left(0\right),S_v\left(0\right),E_v\left(0\right),I_v\left(0\right)\right)\right)\in {\mathbb{R}}_{+}^7$ and for $T\in [0,+\infty ]$ , $\left([0,T],X\left(t\right)=\left(S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right),I_v\left(t\right)\right)\right)$ the maximal solution of the Cauchy problem associated to system $\left(\right)$. Then, for all $t\in [0,T[$, $X\left(t\right)\in {\mathbb{R}}_{+}^7$, $\forall t\in [0,T]$.

Proof. 

Let

$$\Delta =\{\tilde{t}\in [0,T[,S_h\left(t\right)>0,E_h\left(t\right)>0,I_h\left(t\right)>0,R_h\left(t\right)>0,S_v\left(t\right)>0,E_v\left(t\right)>0,I_v\left(t\right)>0\forall t\in ]0,\tilde{t}\left[\right\}.$$

By continuity of functions $S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right)$ and $I_v\left(t\right)$, one can see that $\Delta \ne 0$. Let $\tilde{T}=sup\Delta$. next, we have to show that $\tilde{T}=T$. Suppose $\tilde{T}<T$ , then one has that $S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),$

$R_h\left(t\right),S_v\left(t\right),E_v\left(t\right)$ and $I_v\left(t\right)$ are non negative on $[0,\tilde{T}[$. At $\tilde{T}$ at least one of the following conditions is satisfied $S_h\left(\tilde{T}\right)=0,E_h\left(\tilde{T}\right)=0,I_h\left(\tilde{T}\right)=0,R_h\left(\tilde{T}\right)=0$, $S_v\left(\tilde{T}\right)=0,E_v\left(\tilde{T}\right)=0$ and $I_v\left(\tilde{T}\right)=0$. Then, from the last equation of system $\left(\right)$ one has that:

$$\begin{array}{ccc}{}_{t0}{}^{c}D_t^{\varphi }{I}_v\left(t\right)\hfill & =\hfill & {\alpha }_h{E}_h-({\mu }_v+ϵ)I_v.\hfill \end{array}$$

(4)

Integrating Eq. (4) from 0 to $\tilde{T}$ yields

$$I_v\left(\tilde{T}\right)=I_v\left(0\right)+{\int }_0^{\tilde{T}}({\alpha }_h{E}_h-({\mu }_v+ϵ)I_v)dt>0.$$

Similarly, one can show that $S_h\left(\tilde{T}\right)>0,E_h\left(\tilde{T}\right)>0,I_h\left(\tilde{T}\right)>0,R_h\left(\tilde{T}\right)>0$, $S_v\left(\tilde{T}\right)>0,$, and $E_v\left(\tilde{T}\right)>0$. This is the contradiction then, $\tilde{T}=T$ and consequently the maximal solution for the steady states

$\left(S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right),I_v\left(t\right)\right)$ of the Cauchy problem associated to system $\left(\right)$ is positive. This completes the proof. □

3.2. Boundedness Of Trajectories

We have the following result about the boundedness of trajectories of system (3).

Lemma 1.

Each non-negative solution of system (3) is bounded on ${\mathbb{R}}_{+}^7$.

Proof. 

First of all, we separate model system 3 into two parts: the human population and the Aedes mosquito populations. Adding the first five equations of model system (3) and using equation (1), yields:

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

(5)

Applying the Gronwall Lemma in (5) one get the following:

$$\begin{array}{cc}\hfill N_h\left(t\right)& \le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}}+\left(N_h\left(0\right)-{\displaystyle \frac{{\Lambda }_h}{{\mu }_h}}\right)exp-{\mu }_h.\hfill \end{array}$$

With this in mind, one has that $N_h\left(t\right)\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}}$ for all $t\ge 0$ if $N_{hi}\left(0\right)\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}}$. This means that: $S_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}},E_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}},A_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}},I_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}}$ and $R_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}}$. Using the same concept we proof that $N_{vi}\left(t\right)\le {\displaystyle \frac{{\Lambda }_v}{{\mu }_v}}.$ This mean that: $S_v\le {\displaystyle \frac{{\Lambda }_v}{{\mu }_v}},E_v\le {\displaystyle \frac{{\Lambda }_v}{{\mu }_v}}$ and $,I_v\le {\displaystyle \frac{{\Lambda }_v}{{\mu }_v}}$. This completes the proof □

Corollary 1.

The region

$$\Omega ={\Omega }_H\times {\Omega }_V,$$

(6)

where

$${\Omega }_H=\left\{(S_h,E_h,I_h,R_h)\in {\mathbb{R}}_{+}^4,S_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}},E_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}},A_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}},I_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}},R_h\le {\displaystyle \frac{{\Lambda }_h}{{\mu }_h}}\right\},$$

(7)

and

$${\Omega }_V=\left\{(S_v,E_v,I_v)\in {\mathbb{R}}_{+}^3,S_v\le {\displaystyle \frac{v}{{\mu }_h}},E_v\le {\displaystyle \frac{{\Lambda }_v}{{\mu }_h}},I_v\le {\displaystyle \frac{{\Lambda }_v}{{\mu }_h}}\right\},$$

(8)

is invariant and attractive for system (3).

Thus, system (3) is mathematically and ecologically well posed and it is sufficient to consider the dynamics of the flow generated by system (3) in $\Omega$.

4. Basic Reproduction Number and Existence of Equilibria

In this section we use the next generation method as presented in [9] 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 persist in the population, and 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,A_h^0,I_h^0,R_h^0,S_v^0,E_v^0,I_v^0\right)=\left({\displaystyle \frac{{\Lambda }_h}{{\mu }_h}},0,0,0,0,{\displaystyle \frac{{\Lambda }_v}{{\mu }_v+{ϵ}_v}},0,0\right).\end{array}$$

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

$${\mathcal{F}}_i=\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)$$

(9)

$${\mathcal{V}}_i=\left(\begin{array}{cccc}{\alpha }_h+{\mu }_h+{\theta }_h& 0& 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)$$

(10)

Using (9) and (10) the basic reproduction number ${\mathcal{R}}_0$ of the model system (3) is given by:

$${\mathcal{R}}_0={\displaystyle \frac{(1-\omega ){\delta }_v}{{\mu }_v+ϵ}}\sqrt{{\displaystyle \frac{{\beta }_v{\beta }_h{\Lambda }_v{\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 to be 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 as infectious. The square root here is due to the fact that, the generation of secondary cases in malaria diseases require two transmission process.

Theorem 2.

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 {}_{t0}{}^{c}D_t^{\theta }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)$$

(11)

,

$$\mathcal{V}=\left(\begin{array}{cccc}{\alpha }_h+{\mu }_h+{\theta }_h& 0& 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)$$

(12)

We can observe that by direct calculation we can verify that $V^{-1}F$ is a nonnegative and irreducible matrix and $\rho \left(V^{-1}F\right)={\mathcal{R}}_0.$ It follows from the Perron-Frobenius theorem [12] that $V^{-1}F$ has 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 {}_{t0}{}^{c}D_t^{\varphi }\mathcal{L}\left(t\right)& =& \mathbf{w}{V}^{-1}{}_{t0}{}^{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}$$

(13)

It can be easily verified that the largest invariant subset of $\Gamma$ where ${}_{t0}{}^{c}D_t^{\varphi }\mathcal{L}\left(t\right)=0$ is the singleton $\left\{{\mathcal{E}}^0\right\}$. Therefore, by LaSalle’s invariance principle [14], ${\mathcal{E}}^0$ is globally asymptotically stable in $\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 criteria.

Theorem 3.

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)$$

(14)

and

$$H_2=\left(\begin{array}{ccccccc}{a}_1& 1& 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)$$

(15)

where as $a_j=0$ if $j>n.$ All the roots of the polynomial $p\left(\lambda \right)$ are negative or have negative real part if and if the determinants of all Hurwitz matrices are positive:

$det{H}_j>0,j=1,2,3,...,n.$ When $n=4,$ Routh-Hurwitz criteria simplify 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)$$

(16)

which gives the conditions that, for $n=4:$$a_1>0,a_3>0,a_4>0$ and $a-1a_{2}a-3>a_3^2+a_1^2{a}_4.$

Theorem 4

(The DFE is locally asymptomatically stable if ${\mathcal{R}}_0<1$ and 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)$$

(17)

where by $m_1=-(1-\omega ){\delta }_v{\beta }_h{I}_v-{\mu }_h,m_2=-(1-\omega ){\delta }_h{\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 }_h{\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)$$

(18)

Where by $n_1={\displaystyle \frac{(1-\omega ){\delta }_v{\beta }_h{\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{\Lambda }_v}{{\mu }_v+ϵ}},n_5=({\mu }_v+{\alpha }_v+ϵ)$ The eigenvalues of $J_{DFE}-\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)$$

(19)

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+ϵ)-{\displaystyle \frac{{\alpha }_h{\alpha }_v{(1-\omega )}^2{\delta }_v{\beta }_h{\beta }_v{\Lambda }_h{\Lambda }_v}{{\mu }_h({\mu }_v+ϵ)}}\right)=0\hfill \end{array}$$

For the eigenvalues to be negative then we use the Routh-Hurwitz criterion, then we have that: $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+ϵ)-{\displaystyle \frac{{\alpha }_h{\alpha }_v{(1-\omega )}^2{\delta }_v{\beta }_h{\beta }_v{\Lambda }_h{\Lambda }_v}{{\mu }_h({\mu }_v+ϵ)}}\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 the model system (3) is obtained by simultaneously solving the below system (20) and implicitly expressed in terms of $I_h^{*}$:

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

(20)

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

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

(21)

substituting (20) into equation four of the system (19) and solving for $R_h^{*}$ we obtain:

$$R_h^{*}={\displaystyle \frac{(1-\omega ){\delta }_h{\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)}}$$

(22)

Substituting equation (21) in first equation of the system (19) and solving for $S_h^{*}$ we get the following:

$$S_h^{*}={\displaystyle \frac{\left({\Lambda }_h({\rho }_h+{\mu }_h)+{\gamma }_h{\rho }_h{I}_h^{*}\right)({\mu }_{h}h+{\alpha }_h+{\theta }_h)}{\left((1-\omega ){\delta }_h{\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^{*}}}$$

(23)

Substituting (22) in (20) the simplified expression for (20) becomes:

$$E_h^{*}={\displaystyle \frac{(1-\omega ){\delta }_v{\beta }_v{I}_v^{*}\left({\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^{*}}}$$

(24)

Similarly for the vector population from six of (19) we have:

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

(25)

And from equation five of system (19) we get the following:

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

(26)

substituting (25) in (24) we obtain the simplified expression of $E_v^{*}$ as:

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

(27)

Where as from equation seven of system (19) and equation (26) we have:

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

(28)

Therefore, equations (21)–(27) gives 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 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, as $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 real data of reported cases to align with realistic of the proposed model as shown in Table 1.

5.1. Model Parameter Estimations

In this section, we perform the 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 (Table 1) and some parameters were estimated using root-mean square error (RMSE) the following formula:

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

(29)

Where n is the number of yearly reported malaria cases for 12 years. We assumed the initial population 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 the model (3) the generated new cases is obtained using the term $(1-\omega ){\delta }_v({\beta }_h{S}_h{I}_v+{\beta }_v{S}_v{I}_h)$ which count for the detected cases.

Figure 2 shows the reported malaria cases for 12 years in Zimbabwe, and the simulation results in Figure 3 show (a) Model fit at $\varphi =0.5$ to the real data of malaria cases. From simulations one can note that model (3) had good fit to the reported real cases of malaria presented in Table 2. In figure (b) we fitted the malaria cases with the classical integer and fractional-order model simultaneously and we observe that fractional-order model fits well to the reported cases of malaria compared to classical integer model.

5.2. Sensitivity Analysis

In this section, we use 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. The PRCC results of the ${\mathcal{R}}_0$ related to the number of new cases generated in the population reveals seven key influential parameters. These include the new recruitment of human and vectors (${\Lambda }_h$ and ${\Lambda }_v$), force of infections (${\beta }_h$ and ${\beta }_v$), vector biting rate (${\sigma }_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 recover rate of infected individuals (${\gamma }_h$) have negative correlation with ${\mathcal{R}}_0$, means that whenever these parameters increase the magnitude of ${\mathcal{R}}_0$ decreases and hence the disease dies in the population.

Figure 5 shows global sensitivity analysis of the model (3) on $I_h$ to the key parameters that affect the dynamic of disease in the population. Overall, we noted that increase the parameters that represent rate of insecticides, health education campaign, and vector mortality rate lead to decrease on the number of infected humans in the population. In contrast, we observed that increase on the parameters that represent rate of infections and vector biting rate lead increase on the number infected humans. In particular one can observe that as the rate of insecticides ($ϵ$) and vector mortality rate (${\mu }_v$) above $0.5$ the number of infected humans dies 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 recover rate of exposed vectors), (b) as the function of $ϵ$ (use of insecticides) and $\omega$ (health education campaigns). we observed that when the education campaign is not implemented in the population, the use of insecticides ($ϵ$) must be greater than $0.3$ to reduce ${\mathcal{R}}_0$ to be less than unit. On the other hand, when the education campaign is implemented in the population the rate of insecticides ($ϵ$) must be $0.05$ to reduce the ${\mathcal{R}}_0$ to less than unit. This results demonstrates that to minimize the spread of malaria in the community policy markers must put effort on health education campaigns on how people can prevent themselves in contact with mosquitoes. In addition, we perform the numerical simulation of the system (3) to show the influence of parameters on the model compartments. We used the fourth-order Runge–Kutta numerical method to generate the solution and the results are presented in figure (8).

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 9 and Figure 10 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 that when the fractional order $\varphi =1.0$ the fractional-order model based on Caputo derivative becomes a classical ordinary differential model. Based on the numerical illustrations, one can note that as the order of derivatives $\varphi$ reduced from 1 the memory effects of system increases, and the model solution increase quickly, peak earlier and converge to the unique equilibrium point. In particular, for ${\mathcal{R}}_0<1$ the model solution converge to the disease free equilibrium after 20 days, and for ${\mathcal{R}}_0>1$ the model solutions converge to unique endemic equilibrium. Furthermore, we observed that when the fractional order approaches to 0 the memory effects become strong and the model solutions converge to their respective equilibrium point earlier than when order of derivatives approach to 1. This results was also observed in the previous studies.

5.3. Effects of insecticides use on the disease dynamics

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

5.4. Effects of health education campaigns on the disease dynamics

In Figure 13 we simulate the system (3) to show the effect of health education campaign ($\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 assertions, we numerically simulate the system (3) at $\omega =0$, $\omega =0.5$ and $\omega =0.9$ to demonstrate the effect of health education in minimizing the spread of disease in the population. From numerical illustrations one can note that as $\omega$ increases the number of infection generated decreases. Additionally, we observed that implementing health education campaign alone the magnitude of transmission potential (${\mathcal{R}}_0$) can not be reduced to less than 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 the use of insecticides on the spread of disease in the population. We computed the disease-free equilibrium and derived the reproduction number of the model using the next-generation method. A sensitivity analysis of the basic reproduction number was performed using partial rank correlation coefficients to examine the relationship between the reproduction number and model parameters. Overall, we found that parameters with negative indices cause a decrease in the reproduction number when increased, while parameters with positive indices lead to an increase in the reproduction number whenever. Furthermore, we fitted the model system (3) with real data of malaria cases reported in Zimbabwe for 12 years. We observed that the fractional order model ($\varphi =0.5$) had a good fit compared to classical model ($\varphi =1$). To investigated the impact of memory effects, we simulated the system (3) at the order of derivatives $\varphi :0.5,0.7,0.9$ and 1. overall, we observed that as the order of derivatives have an influence on the spread of disease in the population. In particular, one can note that as the order of derivative decrease from 1 the solution behavior of the system converge to the unique equilibrium much faster compared to when the order of derivative approaches 1. Finally, we simulated the model to assess the effect of proposed interventions (health education campaign and use of insecticides). From numerical results, We observed that implementing insecticides at $90\%$ in the population leads to eradication of disease in the community. On the other hand, we noted that implementing health education campaign alone the disease can no be eradicated in the population. However, the formulated model in this study can be improved in future by incorporation of delays in some of the model parameters and assess the its effect on the spread of disease in the population. Furthermore, this model can be improved by considering migration on both human and vectors.