Model formulationA deterministic compartmental model on the transmission dynamics of malaria is been proposed. We sub-divided the total human population \(N_{H} \left( t \right)\), into five (5) distinct classes of susceptible humans \(S_{H}\), Exposed humans to malaria only \(E_{M}\), Infectious humans with malaria only \(I_{M}\), treated humans due to malaria only \(T_{M}\), and Recovered humans \(R\) from Malaria. We further sub-divided the mosquitoes population into three (3) distinct classes of susceptible mosquitoes \(S_{V}\), exposed mosquitoes \(E_{V}\) and infectious mosquitoes \(I_{V}\). The recruitment rate of humans into the susceptible population is at the rate of \(\Lambda_{H}\), so that \(\beta_{M}\) is the effective contact rate with the probability of infection per contact with an infected mosquito with malaria. The rate at which an individual exposed to malaria become infectious is given as \(\theta_{M}\).The rate of recovery of humans from the malaria is \(\alpha_{M}\). The rate at which a recovered individual from malaria becomes susceptible again is denoted by \(\omega_{M}\). The natural death rate of humans is represented by \(\mu_{H}\), while the disease-induced death rate due to malaria is given as \(\delta_{M}\). The rate of compliance with the use of treated bed nets, aimed at reducing the burden of malaria, is denoted by \(\phi\), and the biting rate of malaria-carrying mosquitoes per unit time is represented by \(m\). Additionally, the recruitment rate of the malaria-transmitting Anopheles mosquitoes is given by \(\Lambda_{V}\), making \(\beta_{V}\) the effective contact rate, defined as the probability of infection per bite from an infected human. Exposed mosquitoes then progress to the infectious stage at a rate denoted by \(\theta_{V}\).The natural death rate of the mosquitoes is given as \(\mu_{V}\) and the mortality due to quest for blood meal by the mosquitoes is given as \(\delta_{V}\). In this study, mosquitoes dynamics are explicitly incorporated into the model to accurately capture the transmission cycle of malaria. The mosquito population is divided into three compartments: susceptible, exposed, and infected vectors. The model accounts for key biological processes, including mosquito biting rates and transmission probabilities, which together define the force of infection between humans and vectors. The latency period within mosquitoes known as the extrinsic incubation period is represented by the progression from the exposed to the infected class. Although natural mortality is included, the total mosquito population is assumed to remain constant over the study period for analytical simplicity. The use of fractional-order derivatives, specifically in the Caputo sense, introduces memory effects into the model, enabling it to capture the temporal dependence of infection dynamics. This fractional framework enhances the biological realism of vector–host interactions, providing a more nuanced understanding of malaria transmission and the impact of intervention strategies.Table 1 presents a comprehensive summary of the variables and parameters utilized in the model, highlighting their roles within the analysis. Each variable is described alongside its definition. The model parameters, including constants and coefficients, define the framework’s relationships and dynamics15. This table serves as a crucial reference, enabling readers to understand the model’s structure and application while promoting transparency and reproducibility.Table 1 Description of model parameters.Full size tableModel compartmental explanations1.Susceptible humans: This compartment includes all individuals who are not currently infected with malaria but can contract the disease if bitten by an infected mosquito. They have no immunity and have not been exposed to the parasite yet34. The transition from this class to the exposed class occurs through infectious bites from mosquitoes.2.Exposed humans : These are individuals who have been bitten by an infected mosquito and have harbored the malaria parasite but are not yet infectious themselves35. This class represents the latent period in humans during which the parasite is incubating in the liver. They eventually transition to the infected class once symptoms develop or the parasite reaches the bloodstream.3.Infected humans: This compartment includes individuals who have developed malaria symptoms and are capable of transmitting the disease to mosquitoes. They are the primary source of infection for susceptible mosquitoes36. Without treatment, they may experience severe illness or even death, depending on the assumptions in the model.4.Treated humans: Individuals in this class have received medical treatment for malaria. Treatment reduces their infectivity and possibly shortens the duration of infection. They may still carry the parasite briefly, depending on treatment efficacy, but are progressing toward recovery37. This class is essential for modeling control strategies and evaluating the impact of healthcare interventions.5.Recovered humans: This compartment consists of individuals who have fully recovered from malaria and have developed immunity, at least temporarily. They are no longer infectious and do not contribute to disease transmission34. Depending on the model structure, immunity may be lifelong or may wane over time, allowing re-entry into the susceptible class.6.Susceptible mosquitoes: These are healthy mosquitoes that have not been exposed to the malaria parasite. They become exposed by biting an infected human38. This is the starting point for mosquitoes participation in the transmission cycle.7.Exposed mosquitoes: These mosquitoes have bitten an infected human and carry the parasite in incubation, but cannot yet transmit it. This corresponds to the extrinsic incubation period (EIP) of the parasite inside the mosquito34. After this period, they become infectious.8.Infected mosquitoes: Infected mosquitoes have completed the incubation period and can now transmit malaria to susceptible humans when they bite34. Infected mosquitoes remain infectious for the rest of their lifespan, as mosquitoes do not recover from malaria.Model equations$$\begin{aligned} \frac{{dS_{H} }}{dt} & = \Lambda_{H} - \left( {\lambda_{M} + \mu_{H} } \right)S_{H} + \omega_{M} R \\ \frac{{dE_{M} }}{dt} & = \lambda_{M} S_{H} - \left( {\theta_{M} + \mu_{H} } \right)E_{M} \\ \frac{{dI_{M} }}{dt} & = \theta_{M} E_{M} - \left( {\gamma_{M} + \delta_{M} + \mu_{H} } \right)I_{M} \\ \frac{{dT_{M} }}{dt} & = \gamma_{M} I_{M} - \left( {\alpha_{M} + \delta_{M} + \mu_{H} } \right)T_{M} \\ \frac{dR}{{dt}} & = \alpha_{M} T_{M} - \left( {\omega_{M} + \mu_{H} } \right)R \\ \frac{{dS_{V} }}{dt} & = \Lambda_{V} - \left( {\lambda_{V} + \mu_{V} } \right)S_{V} \\ \frac{{dE_{V} }}{dt} & = \lambda_{V} S_{V} - \left( {\theta_{V} + \mu_{V} } \right)E_{V} \\ \frac{{dI_{V} }}{dt} & = \theta_{V} E_{V} - \left( {\delta_{V} + \mu_{V} } \right)I_{V} \\ \end{aligned}$$(5)where$$\lambda_{M} = \frac{{\left( {1 - \phi } \right)m\beta_{M} I_{V} }}{{N_{H} }},\;\;\lambda_{V} = \frac{{m\beta_{V} I_{M} }}{{N_{H} }}$$Figure 1 presents a schematic diagram for the malaria transmission model, visually illustrating the flow of individuals through various compartments of the disease process15. The diagram typically includes key compartments such as Susceptible, Exposed, Infectious, and Recovered individuals, along with any relevant intermediate stages38. Arrows between these compartments represent transitions between states, driven by processes like infection transmission (through mosquito bites), progression from exposed to infectious, recovery, and possibly death due to malaria.Fig. 1Schematic diagram for the model.Full size imageAssumptions in the model formulation1.The human and mosquito populations are assumed to mix homogeneously, meaning all individuals have an equal likelihood of contact. Every susceptible human has the same probability of being bitten by an infected mosquito, and every mosquito has the same chance of biting an infected human35. This also implies uniform susceptibility across the population to malaria infection, regardless of geographic, behavioral, or socio-economic factors.2.The rate of malaria transmission from mosquitoes to humans (and vice versa) is considered constant throughout the study period. This assumes that environmental, climatic, and seasonal factors which in reality can affect mosquito activity and transmission rates are not modeled39. Additionally, mosquitoes (e.g., Anopheles species) are assumed to be evenly distributed across the study area, ignoring any spatial heterogeneity.3.The total human and mosquito populations are assumed to remain constant during the study period. This means that birth and natural death rates are balanced or have negligible impact. The model does not account for changes due to migration, human movement between regions, or population growth37. Similarly, malaria-induced mortality may be considered if specified, but general demographic transitions are excluded.4.The model is designed exclusively for malaria transmission and does not consider the presence of other infectious diseases or co-infections. All individuals are either malaria-infected or malaria-free40. There is also no inclusion of vaccination in the population; hence, all susceptible individuals can acquire the infection upon exposure.5.Once exposed to the malaria parasite, infected individuals follow a fixed disease progression pathway: incubation, symptomatic infection, treatment (if accessible), and then either recovery or death41. Variations in disease course due to individual immune response, treatment failure, or drug resistance are not explicitly modeled. Treatment is assumed to be immediately available and uniformly effective across the population.6.Recovered individuals may acquire temporary immunity, but they are not permanently protected. In endemic settings with high transmission rates, they can become re-infected upon re-exposure38. This reflects realistic malaria dynamics where immunity wanes over time and reinfection is common.7.Infected mosquitoes are assumed to remain infectious until they die; they do not recover or lose the ability to transmit malaria37. Furthermore, no interventions directly targeting mosquitoes treatment (e.g., genetic modification or vaccines for mosquitoes) are considered in the model.8.The density of the mosquitoes population is considered constant. Seasonal or intervention-driven changes in mosquito numbers such as those resulting from rainfall, insecticide spraying, or habitat control are not explicitly modeled42.9.It is assumed that diagnostic tools and malaria treatments are equally accessible across the entire human population34. Differences in healthcare access due to geographic, infrastructural, or economic factors are not included, even though they may influence treatment effectiveness in reality.Fractional Malaria mathematical modelIn this section, we expand the integer-based malaria model presented in Eq. (5) by incorporating the Caputo fractional derivative operator. This enhanced model offers greater flexibility compared to the classical model in Eq. (5), as the fractional-order framework allows for adjustable outputs, enabling varied responses. The fractional-order malaria model is thus introduced as follows:$$\begin{aligned} & {}^{a}D_{t}^{\beta } S_{H} = \Lambda_{H} - \left( {\lambda_{M} + \mu_{H} } \right)S_{H} + \omega_{M} R \hfill \\ &{}^{a}D_{t}^{\beta } E_{M} = \lambda_{M} S_{H} - K_{1} E_{M} \hfill \\ & {}^{a}D_{t}^{\beta } I_{M} = \theta_{M} E_{M} - K_{2} I_{M} \hfill \\ &{}^{a}D_{t}^{\beta } T_{M} = \gamma_{M} I_{M} - K_{3} T_{M} \hfill \\ & {}^{a}D_{t}^{\beta } R = \alpha_{M} T_{M} - K_{4} R \hfill \\ & {}^{a}D_{t}^{\beta } S_{V} = \Lambda_{V} - \left( {\lambda_{V} + \mu_{V} } \right)S_{V} \hfill \\ & {}^{a}D_{t}^{\beta } E_{V} = \lambda_{V} S_{V} - K_{5} E_{V} \hfill \\ & {}^{a}D_{t}^{\beta } I_{V} = \theta_{V} E_{V} - K_{6} I_{V} \hfill \\ \end{aligned}$$(6)where \(K_{1} = \left( {\theta_{M} + \mu_{H} } \right)\), \(K_{2} = \left( {\gamma_{M} + \delta_{M} + \mu_{H} } \right)\), \(K_{3} = \left( {\alpha_{M} + \delta_{M} + \mu_{H} } \right)\), \(K_{4} = \left( {\omega_{M} + \mu_{H} } \right)\), \(K_{5}\) = \(\left( {\theta_{V} + \mu_{V} } \right)\), \(K_{6} = \left( {\delta_{V} + \mu_{V} } \right)\).Subject to the positive initial conditions$$\begin{aligned} & S_{H} \left( 0 \right) = S_{H0} ,\,\,\,E_{M} \left( 0 \right) = E_{M0} ,\,\,\,I_{M} \left( 0 \right) = I_{M0} ,\,\,T_{M} \left( 0 \right) = T_{M0} ,\,\,\,R\left( 0 \right) = R_{0} ,\,\,\,S_{V} \left( 0 \right) = S_{V0} \\ & E_{V} \left( 0 \right) = E_{V0} ,\,\;I_{V} \left( 0 \right) = I_{V0} \\ \end{aligned}$$(7)Positivity of model solutionConsidering the non-negativity of the initial values$$\lim Sup \cdot N_{H} \left( t \right) \le \frac{{\Lambda_{H} }}{{\mu_{{^{H} }} }},$$Secondly, If \(\lim Sup \cdot N_{H0} \left( t \right) \le \frac{{\Lambda_{H} }}{{\mu_{H} }},\) then our model feasible domain is given by:$$\begin{aligned} & \Omega = \left\{ {\left( {S_{H} ,E_{M} ,I_{M} ,T_{M} ,R} \right) \subset R_{ + }^{6} :S_{H} + E_{M} + I_{M} + T_{M} + R \le \frac{{\Lambda_{H} }}{{\mu_{H} }},} \right\}\;{\text{and}} \\ & \Omega_{V} = \left\{ {\left( {S_{V} ,E_{V} ,I_{V} } \right) \subset R_{ + }^{3} :S_{V} + E_{V} + I_{V} \le \frac{{\Lambda_{V} }}{{\mu_{V} }}} \right\}\;{\text{so that}} \\ & \Omega = \Omega_{H} X\;\Omega_{V} \subset R_{ + }^{8} , \\ \end{aligned}$$hence \(\Omega\) is positively invariant.If \(S_{H0} ,E_{M0} ,I_{M0} ,T_{M0} ,R_{0} ,S_{V0} ,E_{V0} ,I_{V0}\) are non-negative, then the solution of model (6) will be non-negative for t > 043. From (6), picking the first equation, we have that$$\begin{aligned} & {}^{a}D_{t}^{\beta } S_{H} = \Lambda_{H} - \left( {\lambda_{M} + \mu_{H} } \right)S_{H} + \omega_{M} R \\ & {}^{a}D_{t}^{\beta } S_{H} + \left( {\lambda_{M} + \mu_{H} } \right)S_{H} = \Lambda_{H} + \omega_{M} R \\ \end{aligned}$$But \(\Lambda_{H} + \theta_{M} R \ge 0\), then$${}^{a}D_{t}^{\beta } S_{H} + \left( {\lambda_{M} + \mu_{H} } \right)S_{H} \ge 0.$$Applying the Laplace transform we obtained;$$\begin{aligned} & L\left[ {{}^{a}D_{t}^{\beta } S_{H} } \right] + L\left[ {\left( {\lambda_{M} + \mu_{H} } \right)S_{H} } \right] \ge 0, \\ & S^{\beta } S_{H} \left( s \right) - S^{\beta - 1} S\left( 0 \right) + \left( {\lambda_{M} + \mu_{H} } \right)S_{H} \left( s \right) \ge 0 \\ & S_{H} \left( s \right) \ge \frac{{S^{\beta - 1} }}{{\,\left( {S^{\beta } + \left( {\lambda_{M} + \mu_{H} } \right)} \right)}}S_{H0} . \\ & S_{H} \left( t \right) \ge E_{r,1} \left( { - \left( {\lambda_{H} + \mu } \right)t^{\gamma } } \right)S_{H0} . \\ \end{aligned}$$(8)Now since the term on the right hand side of Eq. (8) is positive, we conclude that \(S_{H} \ge 0\) for \(t \ge 0\) In the same way, we also have that \(E_{M} \ge 0,\;I_{M} \ge 0\), \(T_{M} \ge 0,\;R \ge 0\), \(S_{V} \ge 0,\;E_{V} \ge 0\), \(I_{V} \ge 0,\;R \ge 0\) that is are positives, therefore, the solution will remain in \(R_{ + }^{8}\) for all \(t \ge 0\) with positive initial conditions15,43.Boundedness of fractional model solutionThe total population of individuals from our model is given by:$$N_{H} (t) = S_{H} (t) + E_{M} (t) + I_{M} (t) + T_{H} (t) + R_{H} (t).$$So from our fractional model (6), we now obtain$$\begin{aligned} &^{a} D_{t}^{\beta } N_{H} (t) = \;^{a} D_{t}^{\beta } S_{H} (t) + \;^{a} D_{t}^{\beta } E_{M} (t) + \;^{a} D_{t}^{\beta } I_{M} (t) + \;^{a} D_{t}^{\beta } T_{H} (t) + \;^{a} D_{t}^{\beta } R(t), \\ &^{a} D_{t}^{\beta } N_{H} (t) = \;\Lambda_{H} - (\delta_{M} I_{M} + \delta_{M} T_{M} ) - \mu_{H} N_{H} (t) \le \;\Lambda_{H} - \;\mu_{H} N_{H} (t) \\ &^{a} D_{t}^{\beta } N_{H} (t)\; \le \Lambda_{H} - \mu_{H} N_{H} \left( t \right) \\ \end{aligned}$$(9)Taking the Laplace transformation of (9) we have;$$\begin{aligned} & L\left[ {^{a} D_{t}^{\beta } N_{H} (t)} \right] \le L\left[ {\Lambda_{H} - \mu_{H} N_{H} \left( t \right)} \right]\,, \\ & S^{\beta } N_{H} (s) - S^{\beta - 1} N_{H} (0) + \mu_{H} N_{H} \left( s \right) \le \frac{{\Lambda_{H} }}{S}, \\ & N_{H} (s) \le \;\frac{{S^{\beta - 1} }}{{\left( {S^{\beta } + \mu_{H} } \right)}}N_{H} (0) + \frac{{\Lambda_{H} }}{{S\left( {S^{\beta } + \mu_{H} } \right)}} \\ \end{aligned}$$(10)By taking the inverse Laplace transform of Eq. (10) yields;$$N_{H} (t) \le E_{\beta ,1} \left( { - \mu_{H} t^{\beta } } \right)N_{H} (0) + \Lambda_{H} E_{\beta ,\beta + 1} \left( { - \mu_{H} t^{\beta } } \right)$$(11)At \(t \to \infty\), the limit of (12) becomes$$\mathop {\lim }\limits_{t \to \infty } SupN_{H} \left( t \right) = \frac{{\Lambda_{H} }}{{\mu_{H} }}.$$This means that, if \(N_{H0} \le \frac{{\Lambda_{H} }}{{\mu_{H} }}\) then \(N_{H} \left( t \right) \le \frac{{\Lambda_{H} }}{{\mu_{H} }}\) which implies that, \(N_{H} \left( t \right)\) is bounded.We can also show in the same way that \(N_{V} \left( t \right)\) is bounded.We now conclude that, this region \(\Omega = \Omega_{H} \;X\;\Omega_{V}\), is well posed and equally feasible epidemiologically15,16,44.Existence and uniqueness of solutions of the modelIn this study, we establish the existence and uniqueness of solutions for the fractional-order malaria model using fixed-point theory and standard results applicable to Caputo-type differential equations. Specifically, we show that under suitable conditions such as the Lipschitz continuity and boundedness of the model functions—the system admits a unique solution over a given time interval. These conditions are verified by analyzing the right-hand side of each differential equation, ensuring they satisfy the criteria for contraction mappings in an appropriate function space. The importance of these results lies in the mathematical integrity and predictive reliability of the model42. Existence guarantees that a solution to the system actually exists for the given initial conditions, while uniqueness ensures that the solution is singular and well-defined, meaning the system does not yield multiple conflicting outcomes from the same starting point. This is particularly critical in the context of public health modeling, where decision-making based on simulations must rest on mathematically sound foundations.Furthermore, in the context of fractional-order models, the presence of memory and non-local behavior adds complexity to the system’s dynamics34. Establishing existence and uniqueness in this setting provides confidence that the fractional system behaves consistently and that its solutions are not artifacts of numerical or analytical instability. This, in turn, justifies the use of fractional-order methods in understanding the dynamics of malaria transmission and evaluating the effectiveness of control strategies.Let \({\text{T}}\) be a non-negative number that is real, considering \(M = \left[ {0,{\text{T}}} \right].\)The set of all continuous function that is defined on M is represented by \(N_{e}^{0} \left( M \right)\) with norm as;$$\left\| X \right\| = Sup\left\{ {\left| {X\left( t \right)} \right|,\; {\text{t}} \in M} \right\}.$$Considering model (6) with initial conditions presented in (7) which can be denoted as an initial value problem (IVP) in (12).$$\begin{aligned} & {}^{a}D_{t}^{\beta } \left( t \right) = Z\left( {t,X\left( t \right)} \right),\;0 < t < T < \infty , \\ & X\left( 0 \right) = X_{0} . \\ \end{aligned}$$(12)where \(X\left( t \right) = \left( {S_{H} \left( t \right),E_{M} \left( t \right),\,I_{M} \left( t \right),\,T_{M} \left( t \right),\,R\left( t \right),S_{V} \left( t \right),E_{V} \left( t \right),I_{V} \left( t \right)} \right)\) represents the classes and Z is a continuous function given as;$$Z\left( {t,X\left( t \right)} \right) = \left( {\begin{array}{*{20}c} {Z_{1} \left( {t,S_{H} \left( t \right)} \right)} \\ {Z_{2} \left( {t,E_{M} \left( t \right)} \right)} \\ {Z_{3} \left( {t,I_{M} \left( t \right)} \right)} \\ {Z_{4} \left( {t,T_{M} \left( t \right)} \right)} \\ {Z_{5} \left( {t,R\left( t \right)} \right)} \\ \begin{gathered} Z_{6} \left( {t,S_{V} \left( t \right)} \right) \hfill \\ Z_{7} \left( {t,E_{V} \left( t \right)} \right) \hfill \\ Z_{8} \left( {t,I_{V} \left( t \right)} \right) \hfill \\ \end{gathered} \\ \end{array} } \right) = \left( \begin{gathered} \Lambda_{H} - \left( {\frac{{\left( {1 - \phi } \right)m\beta_{M} I_{V} }}{{N_{H} }} + \mu_{H} } \right)S_{H} + \omega_{M} R \hfill \\ \left( {\frac{{\left( {1 - \phi } \right)m\beta_{M} I_{V} }}{{N_{H} }}} \right)S_{H} - \left( {\theta_{M} + \mu_{H} } \right)E_{M} \hfill \\ \theta_{M} E_{M} - \left( {\gamma_{M} + \delta_{M} + \mu_{H} } \right)I_{M} \hfill \\ \gamma_{M} I_{M} - \left( {\alpha_{M} + \delta_{M} + \mu_{H} } \right)T_{M} \hfill \\ \alpha_{M} T_{M} - K_{4} R \hfill \\ \Lambda_{V} - \left( {\frac{{m\beta_{V} I_{M} }}{{N_{H} }} + \mu_{V} } \right)S_{V} \hfill \\ \left( {\frac{{m\beta_{V} I_{M} }}{{N_{H} }}} \right)S_{V} - K_{5} E_{V} \hfill \\ \theta_{V} E_{V} - K_{6} I_{V} \hfill \\ \end{gathered} \right)\,\,\,$$(13)Using Proposition (1), we have that,$$\begin{aligned} & S_{H} \left( t \right) = S_{H0} + I_{t}^{\beta } \left[ {\Lambda_{H} - \left( {\frac{{\left( {1 - \phi } \right)m\beta_{M} I_{V} }}{{N_{H} }} + \mu_{H} } \right)S_{H} + \omega_{M} R} \right], \\ & E_{M} \left( t \right) = E_{M0} + I_{t}^{\beta } \left[ {\left( {\frac{{\left( {1 - \phi } \right)m\beta_{M} I_{V} }}{{N_{H} }}} \right)S_{H} - \left( {\theta_{M} + \mu_{H} } \right)E_{M} } \right], \\ & I_{M} \left( t \right) = I_{M0} + I_{t}^{\beta } \left[ {\theta_{M} E_{M} - \left( {\gamma_{M} + \delta_{M} + \mu_{H} } \right)I_{M} } \right], \\ & T_{M} \left( t \right) = T_{M0} + I_{t}^{\beta } \left[ {\gamma_{M} I_{M} - \left( {\alpha_{M} + \delta_{M} + \mu_{H} } \right)T_{M} } \right], \\ & R\left( t \right) = R_{0} + I_{t}^{\beta } \left[ {\alpha_{M} T_{M} - K_{4} R} \right], \\ & S_{V} \left( t \right) = S_{V0} + I_{t}^{\beta } \left[ {\Lambda_{V} - \left( {\frac{{m\beta_{V} I_{M} }}{{N_{H} }} + \mu_{V} } \right)S_{V} } \right], \\ & E_{V} \left( t \right) = E_{V0} + I_{t}^{\beta } \left[ {\left( {\frac{{m\beta_{V} I_{M} }}{{N_{H} }}} \right)S_{V} - K_{5} E_{V} } \right], \\ & I_{V} \left( t \right) = I_{V0} + I_{t}^{\beta } \left[ {\theta_{V} E_{V} - K_{6} I_{V} } \right]. \\ \end{aligned}$$(14)We obtain the Picard iteration of (15) as follows;$$\begin{aligned} & S_{Hn} \left( t \right) = S_{H0} + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} A_{1} } \left( {\lambda ,S_{{H\left( {n - 1} \right)}} \left( \lambda \right)} \right)d\lambda , \\ & E_{Mn} \left( t \right) = E_{M0} + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} A_{2} } \left( {\lambda ,E_{{M\left( {n - 1} \right)}} \left( \lambda \right)} \right)d\lambda , \\ & I_{Mn} \left( t \right) = I_{M0} + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} A_{3} } \left( {\lambda ,I_{{M\left( {n - 1} \right)}} \left( \lambda \right)} \right)d\lambda , \\ & T_{Mn} \left( t \right) = T_{M0} + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} A_{4} } \left( {\lambda ,T_{{M\left( {n - 1} \right)}} \left( \lambda \right)} \right)d\lambda , \\ & R_{n} \left( t \right) = R_{0} + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} A_{5} } \left( {\lambda ,R_{n - 1} \left( \lambda \right)} \right)d\lambda , \\ & S_{Vn} \left( t \right) = S_{V} + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} A_{6} } \left( {\lambda ,S_{{V\left( {n - 1} \right)}} \left( \lambda \right)} \right)d\lambda , \\ & E_{Vn} \left( t \right) = E_{V} + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} A_{7} } \left( {\lambda ,E_{{V\left( {n - 1} \right)}} \left( \lambda \right)} \right)d\lambda , \\ & I_{Vn} \left( t \right) = I_{V} + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} A_{8} } \left( {\lambda ,I_{{V\left( {n - 1} \right)}} \left( \lambda \right)} \right)d\lambda . \\ \end{aligned}$$(15)Transforming the initial value problem in Eq. (13) yields;$$X\left( t \right) = X\left( 0 \right) + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} A} \left( {\lambda ,X\left( \lambda \right)} \right)d\lambda .$$(16)Lemma 1The Lipchitz condition described from (13) is satisfied by mosquitoes \(Z\left( {t,X\left( t \right)} \right)\) on a set \(\left[ {0,T} \right] \times R_{ + }^{8}\) with the Lipchitz constant given as;$$\varphi = \max \left( \begin{gathered} \left( {\beta_{M}^{*} + \mu_{H} } \right),\left( {\mu_{H} + \theta_{M} } \right),\left( {\mu_{H} + \gamma_{M} + \delta_{M} } \right),\left( {\mu_{H} + \alpha_{M} + \delta_{M} } \right) \hfill \\ \left( {\mu_{H} + \omega_{M} } \right),\left( {\beta_{V}^{*} + \mu_{V} } \right),\left( {\mu_{V} + \theta_{V} } \right),\left( {\mu_{V} + \delta_{V} } \right). \hfill \\ \end{gathered} \right).$$Proof$$\begin{aligned} & \left\| {Z_{1} \left( {t,S_{H} } \right) - Z_{1} \left( {t,S_{H1} } \right)} \right\| \\ & \quad = \left\| {\left( {\Lambda_{H} - \left( {\frac{{\left( {1 - \phi } \right)m\beta_{M} I_{V} }}{{N_{H} }} + \mu_{H} } \right)S_{H} + \omega_{M} R} \right) - \left( {\Lambda_{H} - \frac{{\left( {1 - \phi } \right)m\beta_{M} I_{V} }}{{N_{H} }} + \mu_{H} } \right)S_{H1} - \omega_{M} R} \right\| \\ & \quad = \left\| { - \left( {\frac{{\left( {1 - \phi } \right)m\beta_{M} I_{V} }}{{N_{H} }}\left( {S_{H} - S_{H1} } \right) + \mu_{H} \left( {S_{H} - S_{H1} } \right)} \right)} \right\| \\ & \quad \le \beta_{M}^{*} \left\| {S_{H} - S_{H1} } \right\| + \mu_{H} \left\| {S_{H} - S_{H1} } \right\| \\ & \qquad \left\| {Z_{1} \left( {t,S_{H} } \right) - Z_{1} \left( {t,S_{H1} } \right)} \right\| \le \beta_{M}^{*} + \mu_{H} \left\| {S_{H} - S_{H1} } \right\|. \\ \end{aligned}$$Similarly, we can say that,$$\begin{aligned} &\left\| {Z_{2} \left( {t,E_{M} } \right) - Z_{2} \left( {t,E_{M1} } \right)} \right\| \le \left( {\mu_{H} + \theta_{M} } \right)\left\| {E_{M} - E_{M1} } \right\|, \hfill \\ &\left\| {Z_{3} \left( {t,I_{M} } \right) - Z_{3} \left( {t,I_{M1} } \right)} \right\| \le \left( {\mu_{H} + \gamma_{M} + \delta_{M} } \right)\left\| {I_{M} - I_{M1} } \right\|, \hfill \\ &\left\| {Z_{4} \left( {t,T_{M} } \right) - Z_{4} \left( {t,T_{M1} } \right)} \right\| \le \left( {\mu_{H} + \alpha_{M} + \delta_{M} } \right)\left\| {T_{M} - T_{M1} } \right\|, \hfill \\ &\left\| {Z_{5} \left( {t,R} \right) - Z_{5} \left( {t,R_{1} } \right)} \right\| \le \left( {\mu_{H} + \omega_{M} } \right)\left\| {R - R_{1} } \right\|, \hfill \\ & \left\| {Z_{6} \left( {t,S_{V} } \right) - AZ\left( {t,S_{V1} } \right)} \right\| \le \left( {\beta_{V}^{*} + \mu_{V} } \right)\left\| {S_{V} - S_{V1} } \right\|, \\ & \left\| {Z_{7} \left( {t,E_{V} } \right) - Z_{7} \left( {t,E_{V1} } \right)} \right\| \le \left( {\mu_{V} + \theta_{V} } \right)\left\| {E_{V} - E_{V1} } \right\|, \\ & \left\| {Z_{8} \left( {t,I_{V} } \right) - Z_{8} \left( {t,I_{V1} } \right)} \right\| \le \left( {\mu_{V} + \delta_{V} } \right)\left\| {I_{V} - I_{V1} } \right\|, \\ \end{aligned}$$(17)Therefore, we have$$\begin{aligned} & \left\| {Z\left( {t,X_{1} \left( t \right)} \right) - Z\left( {t,X_{2} \left( t \right)} \right)} \right\| \le \varphi \left\| {X_{1} - X_{2} } \right\|, \\ & \varphi = \max \left( \begin{gathered} \left( {\beta_{M}^{*} + \mu_{H} } \right),\left( {\mu_{H} + \theta_{M} } \right),\left( {\mu_{H} + \gamma_{M} + \delta_{M} } \right),\left( {\mu_{H} + \alpha_{M} + \delta_{M} } \right) \hfill \\ \left( {\mu_{H} + \omega_{M} } \right),\left( {\beta_{V}^{*} + \mu_{V} } \right),\left( {\mu_{V} + \theta_{V} } \right),\left( {\mu_{V} + \delta_{V} } \right). \hfill \\ \end{gathered} \right). \\ \end{aligned}$$(18)Lemma 2The initial value problem (6), (7) in (18) exists and will have a unique solution$$X\left( t \right) \in B_{c}^{0} \left( f \right).$$Using Picard–Lindel of and fixed-point theory, we consider the solution of$${\text{X}}\left( t \right) = S\left( {{\text{X}}\left( t \right)} \right),$$where S is defined as the Picard operator expressed as;$${\text{S}}:{\text{B}}_{c}^{0} \left( {f,R_{ + }^{8} } \right) \to B_{c}^{0} \left( {f,R_{ + }^{8} } \right)$$Therefore$$S\left( {X\left( t \right)} \right) = X\left( 0 \right) + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} Z} \left( {\lambda ,X\left( \lambda \right)} \right)d\lambda .$$which becomes$$\begin{aligned} & \left\| {S\left( {X_{1} \left( t \right)} \right) - S\left( {X_{2} \left( t \right)} \right)} \right\| \\ & \quad = \left\| {\frac{1}{\Gamma \left( \beta \right)}\left[ {\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} Z\left( {\lambda ,X_{1} \left( \lambda \right)} \right) - Z\left( {\lambda ,X_{2} \left( \lambda \right)} \right)d\lambda } } \right]} \right\| \\ & \quad \le \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} } \left\| {Z\left( {\lambda ,X_{1} \left( \lambda \right)} \right) - Z\left( {\lambda ,X_{2} \left( \lambda \right)} \right)d\lambda } \right\|. \\ & \quad \le \frac{\varphi }{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - \lambda } \right)^{\beta - 1} } \left\| {X_{1} - X_{2} } \right\|d\lambda . \\ & \qquad \left\| {S\left( {X_{1} \left( t \right)} \right) - S\left( {X_{2} \left( t \right)} \right)} \right\| \le \frac{\varphi }{{\Gamma \left( {\beta + 1} \right)S}}. \\ \end{aligned}$$When \(\frac{\varphi }{{\Gamma \left( {\beta + 1} \right)}}S \le 1\), then the Picard operator gives a contradiction, so (6), (7) solution is unique.The basic reproduction number and model equilibrium pointsThe basic reproduction number \(\left( {R_{0} } \right)\) is computed using the next-generation matrix method, which is based on the linearization of the system around the disease-free equilibrium. Specifically, the matrix is constructed by identifying new infection terms and transition terms among the infected compartments, and \(R_{0}\) is obtained as the spectral radius (dominant eigenvalue) of the resulting next-generation matrix. Although the model is formulated using fractional-order derivatives, the computation of \(R_{0}\) remains structurally consistent with the classical approach, as it depends on the infection and transition structure of the model rather than the order of the derivative. However, the fractional nature of the model influences how the disease evolves over time. In particular, memory effects embedded in the fractional derivatives may lead to slower convergence to equilibrium or prolonged persistence of infection, even when \(R_{0}\) 1\), giving us:$$\lambda_{M}^{*} = \frac{{ - b + \sqrt {b^{2} - 4ac} }}{2a}$$This is the unique endemic equilibrium solution when \({\mathcal{R}}_{0} > 1\).Sensitivity analysis of the malaria modelTo evaluate the impact of model parameters on malaria transmission, we conducted a local sensitivity analysis of the basic reproduction number, R0, using normalized forward sensitivity indices. This method involves computing the partial derivatives of R0 with respect to each parameter, scaled by the parameter value and R0 itself, to determine their relative influence. The results indicated that parameters such as the mosquito biting rate, transmission probability, and treatment rate had the most significant effect on R0. These findings underscore the importance of enhancing treatment coverage and implementing mosquitoes control strategies as effective measures for reducing malaria transmission (Table 2).Table 2 Sensitivity indices of various parameters.Full size tableThe sensitivity index of \(R_{0}\) with respect to a parameter p is given by:$$\Im_{p}^{{R_{0} }} = \frac{p}{{R_{0} }} \cdot \frac{{\partial R_{0} }}{\partial p}$$Given that:$$R_{0} = \sqrt {\frac{{(1 - \phi )m^{2} \beta_{M} \beta_{V} \theta_{M} \theta_{V} \mu_{H}^{2} }}{{\Lambda_{H}^{2} \left( {\theta_{M} + \mu_{H} } \right)\left( {\gamma_{M} + \delta_{M} + \mu_{H} } \right)\left( {\theta_{V} + \mu_{V} } \right)\left( {\delta_{V} + \mu_{V} } \right)}}}$$Figure 2a presents the sensitivity indices for various parameters in the malaria model, highlighting the relative importance of each parameter on the model’s output15. The sensitivity indices are typically represented in a bar chart format, where each bar corresponds to a different model parameter (e.g., transmission rate, recovery rate, mosquito biting rate, etc.).Fig. 2(a) Sensitivity indices for various parameters. (b) Data fitting for the malaria sub-model.Full size imageThe sensitivity analysis of the basic reproduction number \(R_{0}\) reveals the relative influence of key model parameters on malaria transmission dynamics. From the sensitivity indices displayed in the bar chart, it is evident that certain parameters have a positive influence on \(R_{0}\). Specifically, the biting rate of mosquitoes and the contact rate between susceptible humans and infected mosquitoes exhibit positive sensitivity indices, indicating that an increase in these parameters leads to a corresponding increase in the basic reproduction number. This suggests that any measures aimed at reducing mosquito bites such as the use of insecticide-treated nets, indoor residual spraying, or environmental management would be effective in lowering disease transmission. On the other hand, parameters such as the treatment rate and recovery rate of infected individuals exhibit negative sensitivity indices, meaning that increases in these parameters contribute to a decrease in \(R_{0}\). This underscores the importance of promoting timely and effective treatment of malaria cases, as well as improving access to healthcare services, to reduce the overall disease burden in the population. The sensitivity analysis highlights which parameters are most critical to target for effective malaria control and provides insights for prioritizing intervention strategies.Fractional order model numerical resultsUsing the generalized fractional Adams–Bashforth–Moulton method outlined in Eq. (23), we conducted a numerical solution for the fractional-order malaria model. The parameter values utilized in the simulations are detailed in Table 1, which includes a range of fractional-order values43,44. This varying fractional order \(\beta\) was systematically explored to observe their impact on the model’s behavior, providing insights into the dynamics of the disease under different scenarios.Implementation of fractional Adams–Bashforth methodThis study employs the method introduced by Baskonus et al. and further developed by Diethelm and Freed, as documented in references24,43,44,45,46. Using these approaches, the fractional Adam–Bashforth–Moulton method is applied to derive an approximate solution for the fractional malaria model initially presented in Eq. (6). This numerical technique is particularly suited for fractional differential equations, as it captures the intricate dynamics of the system more effectively than traditional methods. For the numerical simulation of the fractional-order malaria model, we adopted the Adams–Bashforth–Moulton (ABM) method, a widely used predictor–corrector approach for solving fractional differential equations. This method effectively handles the nonlocal properties of Caputo-type fractional derivatives, which are integral to capturing memory effects in disease dynamics. The ABM method was chosen for its balance of computational efficiency and numerical accuracy, making it well-suited for the complexity of our model. Accuracy was ensured by employing sufficiently small step sizes and verifying numerical stability through convergence checks across multiple time-step refinements. Consequently, the fractional malaria model described in Eq. (6) is reformulated and expressed with greater precision and flexibility in Eq. (21), as follows:$$\begin{aligned} & {}^{a}D_{t}^{\beta } B\left( t \right) = N\left( {t,b\left( t \right)} \right),\,\,0 < t < \psi \\ & B^{\left( n \right)} \left( 0 \right) = B_{0}^{{^{\left( n \right)} }} ,\,\,\,n = 1,0,\,\,...,b,b = \left[ \beta \right]. \\ \end{aligned}$$(21)where \(B = (S_{H}^{*} ,E_{M}^{*} ,I_{M}^{*} ,T_{M}^{*} ,R^{*} ,S_{V}^{*} ,E_{V}^{*} ,I_{V}^{*} ) \in R_{ + }^{8}\) and \(Q\left( {t,b\left( t \right)} \right)\) is a real valued-function that is continuous.Equation (21) can be therefore be represented using the concept of fractional integral as follows;$$B\left( t \right) = \sum\limits_{n = 0}^{m - 1} {B_{0}^{\left( n \right)} } \frac{{t^{n} }}{n!} + \frac{1}{\Gamma \left( \beta \right)}\int_{0}^{t} {\left( {t - y} \right)}^{\beta - 1} R\left( {y,b\left( y \right)} \right)dy\,\,$$(22)Applying the Approach described in7, Suppose the step size \(h = \frac{\psi }{N},\,\,N \in N\) with a grid that is uniform on \(\left[ {0,\psi } \right].\) Where \(t_{c} = cr,\,\,c = 0,1,\;.\;.\;.N.\) Therefore, the fractional-order model of malaria model presented in (6) can be approximated as:$$\begin{aligned} S_{{H\left( {k + 1} \right)}} \left( t \right) & = S_{H0} + \frac{{h^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\left\{ {\lambda - \Lambda_{H} - \left( {\left( {1 - \phi } \right)m\beta_{M} I_{V}^{n} } \right)\frac{{S_{H}^{n} }}{{N_{H}^{n} }} - \mu_{H} S_{H}^{n} + \omega_{M} R^{n} } \right\} \\ & \quad + \frac{{g^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\sum\limits_{y = 0}^{k} {dy,k + 1\left\{ {\lambda - \Lambda_{H} - \left( {\left( {1 - \phi } \right)m\beta_{M} I_{Vy} } \right)\frac{{S_{Hy} }}{{N_{Hy} }} - \mu_{H} S_{Hy} + \omega_{M} R_{y} } \right\},} \\ E_{{H\left( {k + 1} \right)}} \left( t \right) & = E_{H0} + \frac{{h^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\left\{ {\left( {\left( {1 - \phi } \right)m\beta_{M} I_{V}^{n} } \right)\frac{{S_{H}^{n} }}{{N_{H}^{n} }} - \left( {\theta_{M} + \mu_{H} } \right)E_{H}^{n} } \right\} \\ & \quad + \frac{{g^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\sum\limits_{y = 0}^{k} {dy,k + 1\left\{ {\left( {\left( {1 - \phi } \right)m\beta_{M} I_{Vy} } \right)\frac{{S_{Hy} }}{{N_{Hy} }} - \left( {\theta_{M} + \mu_{H} } \right)E_{My} } \right\}} , \\ I_{{M\left( {k + 1} \right)}} \left( t \right) & = I_{M0} + \frac{{h^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\left\{ {\theta_{M} E_{M}^{n} - \left( {\gamma_{M} + \delta_{M} + \mu_{H} } \right)I_{M}^{n} } \right\} \\ & \quad + \frac{{g^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\sum\limits_{y = 0}^{k} {dy,k + 1\left\{ {\theta_{M} E_{My} - \left( {\gamma_{M} + \delta_{M} + \mu_{H} } \right)I_{My} } \right\},} \\ T_{{M\left( {k + 1} \right)}} \left( t \right) & = T_{M0} + \frac{{h^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\left\{ {\gamma_{M} I_{M}^{n} - \left( {\alpha_{M} + \delta_{M} + \mu_{H} } \right)T_{M}^{n} } \right\} \\ & \quad + \frac{{g^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\sum\limits_{y = 0}^{k} {dy,k + 1\left\{ {\gamma_{M} I_{My} - \left( {\alpha_{M} + \delta_{M} + \mu_{H} } \right)T_{My} } \right\},} \\ R_{k + 1} \left( t \right) & = T_{H0} + \frac{{h^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\left\{ {\alpha_{M} T_{M}^{n} - K_{4} R^{n} } \right\} \\ & \quad + \frac{{g^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\sum\limits_{y = 0}^{k} {dy,k + 1\left\{ {\alpha_{M} T_{My} - K_{4} R_{y} } \right\}} , \\ S_{{V\left( {k + 1} \right)}} \left( t \right) & = S_{V0} + \frac{{h^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\left\{ {\Lambda_{V} - \left( {m\beta_{V} I_{M}^{n} } \right)\frac{{S_{V}^{n} }}{{N_{H}^{n} }} - \mu_{V} S_{V}^{n} } \right\} \\ & \quad + \frac{{g^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\sum\limits_{y = 0}^{k} {dy,k + 1\left\{ {\Lambda_{V} - \left( {m\beta_{V} I_{My} } \right)\frac{{S_{Vy} }}{{N_{Hy} }} - \mu_{V} S_{Vy} } \right\},} \\ E_{{V\left( {k + 1} \right)}} \left( t \right) & = E_{V0} + \frac{{h^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\left\{ {\left( {m\beta_{V} I_{M} } \right)\frac{{S_{V} }}{{N_{H} }} - K_{5} E_{V} } \right\} \\ & \quad + \frac{{g^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\sum\limits_{y = 0}^{k} {dy,k + 1\left\{ {\left( {m\beta_{V} I_{M} } \right)\frac{{S_{V} }}{{N_{H} }} - K_{5} E_{V} } \right\},} \\ I_{{V\left( {k + 1} \right)}} \left( t \right) & = I_{V0} + \frac{{h^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\left\{ {\theta_{V} E_{V}^{n} - K_{6} I_{V}^{n} } \right\} \\ & \quad + \frac{{g^{\gamma } }}{{\Gamma \left( {\beta + 2} \right)}}\sum\limits_{y = 0}^{k} {dy,k + 1\left\{ {\theta_{V} E_{Vy} - K_{6} I_{Vy} } \right\}.} \\ \end{aligned}$$(23)where$$\begin{aligned} & S_{{H\left( {k + 1} \right)}}^{n} \left( t \right) = S_{H0} + \frac{1}{\Gamma \left( \beta \right)}\sum\limits_{y = 0}^{k} {f_{y,k + 1} } \left\{ {\lambda - \Lambda_{H} - \left( {\left( {1 - \phi } \right)m\beta_{M} I_{Vy} } \right)\frac{{S_{Hy} }}{{N_{Hy} }} - \mu_{H} S_{Hy} + \omega_{M} R_{y} } \right\}, \\ & E_{{M\left( {k + 1} \right)}}^{n} \left( t \right) = E_{M0} + \frac{1}{\Gamma \left( \beta \right)}\sum\limits_{y = 0}^{k} {f_{y,k + 1} } \left\{ {\left( {\left( {1 - \phi } \right)m\beta_{M} I_{Vy} } \right)\frac{{S_{Hy} }}{{N_{Hy} }} - \left( {\theta_{M} + \mu_{H} } \right)E_{My} } \right\}, \\ & I_{{M\left( {k + 1} \right)}}^{n} \left( t \right) = I_{M0} + \frac{1}{\Gamma \left( \beta \right)}\sum\limits_{y = 0}^{k} {f_{y,k + 1} } \left\{ {\theta_{M} E_{My} - \left( {\gamma_{M} + \delta_{M} + \mu_{H} } \right)I_{My} } \right\}, \\ & T_{{M\left( {k + 1} \right)}}^{n} \left( t \right) = T_{M0} + \frac{1}{\Gamma \left( \beta \right)}\sum\limits_{y = 0}^{k} {f_{y,k + 1} } \left\{ {\gamma_{M} I_{My} - \left( {\alpha_{M} + \delta_{M} + \mu_{H} } \right)T_{My} } \right\}, \\ & R_{{\left( {k + 1} \right)}}^{n} \left( t \right) = R_{0} + \frac{1}{\Gamma \left( \beta \right)}\sum\limits_{y = 0}^{k} {f_{y,k + 1} } \left\{ {\alpha_{M} T_{My} - K_{4} R_{y} } \right\}, \\ & S_{{V\left( {k + 1} \right)}}^{n} \left( t \right) = S_{V0} + \frac{1}{\Gamma \left( \beta \right)}\sum\limits_{y = 0}^{k} {f_{y,k + 1} } \left\{ {\Lambda_{V} - \left( {m\beta_{V} I_{My} } \right)\frac{{S_{Vy} }}{{N_{Hy} }} - \mu_{V} S_{Vy} } \right\}, \\ & E_{{V\left( {k + 1} \right)}}^{n} \left( t \right) = E_{V0} + \frac{1}{\Gamma \left( \beta \right)}\sum\limits_{y = 0}^{k} {f_{y,k + 1} } \left\{ {\left( {m\beta_{V} I_{M} } \right)\frac{{S_{V} }}{{N_{H} }} - K_{5} E_{V} } \right\}, \\ & I_{{V\left( {k + 1} \right)}}^{n} \left( t \right) = I_{V0} + \frac{1}{\Gamma \left( \beta \right)}\sum\limits_{y = 0}^{k} {f_{y,k + 1} } \left\{ {\theta_{V} E_{Vy} - K_{6} I_{Vy} } \right\}. \\ \end{aligned}$$(24)From (23) and (24) obtained;$$\begin{aligned} & dy,_{K + 1} = K^{\beta + 1} - \left( {k - \beta } \right)\left( {k + \beta } \right)^{\beta } ,\,\,y = 0 \\ & \left( {k - y + 2} \right)^{\beta + 1} + \left( {k - \beta } \right)^{\beta + 1} - 2\left( {k - y + 1} \right)^{{^{\beta + 1} }} ,\,\,1 \le y \le k \\ & 1,\;\;y = k + 1 \\ \end{aligned}$$and$$f_{y,k + 1} = \frac{{h^{\beta } }}{\beta }\left[ {\left( {k - y + 1} \right)^{\beta } \left( {k - y} \right)^{\beta } } \right],\,\,0 \le y \le k.$$Fitting data for malaria modelsFor malaria, the dataset covers annual records from 2014 to 2021, detailing the number of confirmed malaria cases. These data points are organized and presented in the accompanying table for clarity This information is also summarized in the corresponding table below, offering a comprehensive view of the disease trends over time.Table 3 presents the reported malaria cases in Nigeria from 2010 to 2022, with data sourced from the World Health Organization (WHO). The table provides a year-by-year breakdown of the number of malaria cases reported in the country, offering insight into the trends and patterns of malaria incidence over time47.Table 3 Presents the reported malaria cases from 2010 to 2022 in Nigeria. These figures were sourced from the World Health Organization (WHO) 47. Source: World Health Organization47.Full size tableData fitting methodologyTo ensure the collected data aligns with the mathematical sub-models for malaria and tuberculosis, the fmincon toolbox from MATLAB’s Optimization Toolbox was utilized. This optimization approach facilitates the accurate fitting of model parameters, enhancing the model’s ability to reflect real-world data more precisely. By minimizing discrepancies between the model and the observed data, the fmincon toolbox helps refine the parameter values, ensuring that the simulations accurately represent the underlying dynamics of the diseases. The resulting fitted data plots are presented in the figures below18, showcasing the model’s improved agreement with the observed infection patterns.Figure 2b shows the data fitting results for the malaria sub-model, comparing the model’s simulated outcomes with actual reported malaria case data over a specified period15,47. The figure typically includes a curve or line representing the model’s predictions alongside data points representing the observed malaria cases. A good fit between the simulated data and the observed data indicates that the model accurately captures the dynamics of malaria transmission in the study area.Table 4 lists the parameter values used in the numerical simulations, along with their units. These values, drawn from data or literature, were chosen to realistically represent the system’s behavior15,16. The table provides a clear reference, ensuring transparency and helping other researchers build on or validate the findings. Parameter values were obtained from published literature, WHO reports, and regional malaria data. For parameters lacking direct data, we applied the least squares method to fit the model to reported malaria incidence, ensuring realistic and accurate estimates. The fitting was performed within biologically plausible ranges, and sensitivity analysis was used to validate the influence of key parameters. A full list of parameters, their descriptions, and sources or fitted values is provided in Table 4.Table 4 Parameters values used for Numerical Simulation.Full size tableNumerical simulationNumerical simulations were carried out using MATLAB to analyze the dynamics of malaria and tuberculosis models15,16. The fmincon optimization toolbox helped fit model parameters to real-world data, ensuring accurate representation. Fractional-order terms were used to capture long-term dynamics more effectively than traditional models18. Results from these simulations are illustrated through plots, providing insights into disease trends and supporting public health decision-making.Figure 3a indicates that a rise in the parameters \(\gamma\) and \(\theta\) representing treatment rate of infected humans with malaria and Progression rate of exposed malaria to infected malaria class respectively, which increases the incidence of malaria in human populations, as demonstrated by the basic reproduction number \(R_{0}\) surpassing one (1). Figure 3b shows the effects of different between \(\alpha\) and \(\theta\), indicating the recovery rate of individuals with malaria and Progression rate of exposed malaria to infected malaria class in relation to the population’s basic reproduction number \(R_{0}\). The figure showed that an increase in these parameters leads to a peak in the value of \(R_{0}\) to a numerical value below one. This implies that the burden of malaria among humans will eventually be reduced by limiting rate at people are exposed to malaria and putting policies in place to encourage proper treatment. Figure 4a shows the impact of varying \(\delta_{h}\) (Disease induced death rate of infected malaria) and \(\mu\) (Natural death rate of humans) on the basic reproduction number \(R_{0}\). If appropriate measures are not implemented, these factors \(\delta_{h}\) and \(\mu\) can worsen the prevalence of malaria, as evidenced by their effect on the basic reproduction number, causing it to rise above one (1).Fig. 3Surface plot of the basic reproduction number.Full size imageFig. 4Surface plot and Effect of \(\beta\) on \(S_{H} \left( t \right)\).Full size imageFigure 4b illustrates how changes in the fractional-order derivative \(\beta\) impact the susceptible human population over time. The graph shows that as \(\beta\) increases, the number of susceptible individuals decreases. This decline happens because the treatment strategy reduces the number of people at risk of contracting the disease.Fig. 5a demonstrates that the exposed population decreases rapidly as \(\beta\) rise. This shows that even though the overall exposed population is shrinking due to higher infectiousness, there is still a steady influx of individuals transitioning from the susceptible group to the exposed group. Figure 5b, shows that number of infected human population depicted over time. Initially, the number of infected individuals grows due to the flow of people from the exposed group. However, after a certain point, this number starts to decline, thanks to treatment interventions. This highlights the effectiveness of treatment strategies in ultimately reducing the malaria burden within the population. Figures 3 and 4 provide valuable epidemiological insights into how key parameters influence malaria transmission dynamics. Specifically, the Fig. 3a demonstrates that a simultaneous increase in the treatment rate of infected individuals and the progression rate from exposed to infected classes can drive the basic reproduction number \(R_{0}\) above the critical threshold of one. Epidemiologically, this highlights the unintended consequence of rapid disease progression outpacing treatment efficacy suggesting that without timely diagnosis and efficient treatment coverage, higher infection rates may persist despite increased intervention efforts. In contrast, the Fig. 3b shows that enhancing recovery rates while reducing progression from exposure to infection can successfully suppress \(R_{0}\) below unity. This suggests that public health strategies prioritizing early diagnosis and improved therapeutic effectiveness can significantly curb malaria spread. Furthermore, Fig. 4a indicates that high disease-induced and natural death rates contribute to increased \(R_{0}\), reflecting a population under high mortality stress, where the burden of malaria exacerbates overall vulnerability. Figure 4b, showing the effects of fractional-order changes, reveals that disease memory and population behavior over time significantly influence the susceptible pool. As the fractional-order increases, the model captures more long-term interactions and delayed effects, suggesting that sustained treatment campaigns and long-term community health education are essential for reducing susceptibility and preventing resurgence. These findings emphasize the importance of interpreting parameter interactions not just numerically but in terms of real-world health outcomes, thereby making the model outputs more actionable for public health planning. Figure 6a illustrates the changes in the number of individuals receiving treatment over time. At first, the treated population increases as people transition from the infected group to receiving care. Over time, however, this number gradually decreases as successful medical interventions result in recovery from malaria. Figure 6b shows a rise in the recovered population as treatment rates improve. This highlights the effectiveness of healthcare strategies in combating the disease. However, the recovered population eventually starts to decline since individuals who recover from malaria may become susceptible again, re-entering the cycle of disease transmission. Figure 7 illustrates the cumulative number of new malaria cases in relation to the treatment rate between susceptible and infected individuals. The graph shows that as the treatment rate increases, the cumulative number of new cases begins to decline. This demonstrates that higher treatment rates effectively reduce the overall burden and prevalence of malaria within the population.Fig. 5Effect of \(\beta\) on \(E_{m} \left( t \right)\) and \(I_{m} \left( t \right)\).Full size imageFig. 6Effect of \(\beta\) on \(T_{m} \left( t \right)\) and \(R\left( t \right)\).Full size imageFig. 7Effect of \(\sigma\) on the cumulative new cases of Malaria.Full size imageTo assess the reliability and real-world applicability of the developed model, we compared our simulation results with empirical malaria incidence data from the Nigeria Centre for Disease Control (NCDC) and the World Health Organization (WHO). According to the WHO’s 2023 World Malaria Report, Nigeria remains the country with the highest malaria burden globally, accounting for approximately 27% of global cases and 31% of malaria-related deaths. The NCDC reports a national incidence rate of approximately 301 cases per 1000 population, with significant seasonal and regional variation. Our simulation results demonstrate strong qualitative agreement with these observations. For instance, Figs. 3, 4, 5, 6 and 7 illustrate that increases in treatment rates and reductions in exposure and progression rates can bring the basic reproduction number \(R_{0}\) below unity, indicating potential for malaria elimination. This aligns with empirical trends in Nigeria, where improved treatment coverage has led to measurable declines in incidence in urban areas such as Lagos and Abuja. Conversely, our model also captures scenarios where delayed treatment or high exposure sustains \(R_{0}\) > 1, consistent with persistent transmission in rural and underserved regions. Furthermore, the model reflects reinfection dynamics and treatment effects that mirror reported patterns in Nigeria, where individuals frequently experience multiple infections due to incomplete immunity and ongoing transmission. The inclusion of fractional-order derivatives captures memory effects and non-linear transitions in population behavior, which correspond well to real-world delays in treatment-seeking and gradual changes in susceptibility.