Fractional-order analysis of temperature- and rainfall-dependent mathematical model for malaria transmission dynamics

Malaria remains a substantial public health challenge and economic burden globally. Currently, malaria has been declared as endemic in 85 countries. In this study, we developed and analyzed a fractional-order mathematical model for malaria transmission dynamics that incorporates variability of temperature and rainfall using Caputo-type AB operators. The existence and uniqueness of the model's solutions were established using the Banach fixed-point theorem. The model system's equilibria (both disease-free and endemic) were identified, and lemmas and theorems were developed to prove their stability. Furthermore, we used different temperature ranges and rainfall data, validating them against existing literature. Numerical simulations using the Toufik-Atangana schemes with various fractional-order alpha values revealed that as the value of alpha approaches 1, the behavior of the fractional-order model converges to that of the classical model. The numerical results are promising and are expected to be valuable for future research related to fractional-order models.

1 Introduction

Malaria, a life-threatening mosquito-borne disease, ranks among the deadliest infectious diseases worldwide [see [1] and references therein]. The World Health Organization (WHO) data for 2022 indicate ~249 million malaria cases, leading to 608,000 deaths in 85 countries [2]. This mosquito-borne infectious disease is the fifth leading cause of death from infectious diseases globally (after respiratory infections, HIV/AIDS, diarrheal diseases, and tuberculosis) and the second leading cause in Africa after HIV/AIDS [3]. Moreover, malaria remains a significant cause of mortality and morbidity in many tropical and subtropical regions, particularly within developing nations such as Ethiopia [see [4] and references therein].

In Ethiopia, with a population exceeding 120 million, more than 60% face the risk of contracting malaria [4]. This vulnerability stems from the fact that nearly 70% of Ethiopian land falls within areas suitable for malaria transmission, with altitude and rainfall serving as key risk factors.

Mathematical models of malaria transmission dynamics are valuable tools for understanding the disease, planning for the future, and implementing effective control measures [5]. The first integer-order model of this kind emerged from the study of Ross [6] and Macdonald [7]. While researchers have expanded these models over the years to incorporate various aspects of malaria transmission and control, research has given less attention to the impact of age structure. Recent studies have highlighted the significance of age structure in both the vector population and human population to understand the impact of climate variability on malaria transmission dynamics [see [8–10] and the references therein]. This is because environmental and climatic factors play a significant role in influencing the vector population's dynamics and the biting rate of mosquitoes on humans.

Temperature and rainfall are key drivers of mosquito population dynamics and subsequent malaria transmission, as reported in [11] and the references therein. These factors positively or negatively impact malaria transmission [12, 13]. Studies using real data from 67 sub-Saharan African cities have shown that the optimal temperature range for mosquito growth, and therefore malaria transmission, falls between 16 and 28°C [14]. This link between climate and malaria burden is further supported by data from South Africa's KwaZulu-Natal province [8]. Their findings reveal a clear trend: Within specific ranges, both increasing mean monthly temperature and rainfall are associated with increased malaria burden. Malaria burden rises with temperatures between 17 and 25°C and rainfall between 32 and 110 mm, while decreasing outside these ranges. Interestingly, Okuneye and Gumel [8] also pinpoints the specific temperature and rainfall combinations where malaria transmission peaks. Their data show that the highest transmission rates occur within the ranges of 21–25°C for temperature and 95–125 mm for rainfall. In addition to impacting transmission, temperature also directly affects the lifespan of adult female mosquitoes, which averages ~21 days [see [15] and the reference therein]. Notably, mosquito survival rapidly declines as temperatures exceed 34°C.

According to current literature and researchers' findings, all these models employ integer-order derivatives in their differential equations. However, fractional calculus has become increasingly prevalent in epidemic modeling [see [16–19] and as well some of the references therein]. This fractional approach has demonstrated significant advantages over integer-order modeling, offering a better fit to real data and possessing numerous other beneficial properties [20]. Furthermore, the memory and inheritance features of fractional calculus make it particularly well-suited for modeling and understanding real-world phenomena [20]. Within this framework, a variety of definitions and operators exist, including the Atangana–Baleanu [21], Caputo–Fabrizio [22], and Caputo derivatives [23], and serving as valuable tools for epidemic disease modeling.

Among these, the Caputo and Caputo–Fabrizio derivatives hold greater significance, with the latter particularly excelling due to its non-singular, non-local core and enhanced ability to reflect disease dynamics [[24] and the references within]. While the Atangana–Baleanu operator has found applications in modeling various real-world problems [24–28], its suitability for disease modeling specifically requires further evaluation.

Motivated by previous studies, we present a novel mathematical model of malaria transmission dynamics. Building upon the study of [8], our model analyzes and extends their framework by incorporating specific new compartments while omitting others. We further employ ABC fractional operators to explore the model's dynamics in a non-integer-order setting.

2 Mathematical preliminaries

This section presents important theorems and definitions of fractional calculus, especially some fundamental ideas about the Atangana–Baleanu fractional derivative operators and other related findings, before applying them to our suggested malaria model.

Definition 1. Let Ω ⊆ R be open and p ∈ [1;∞), so H^p(Ω) can be defined as

$\begin{array}{ll}{H}^p(\Omega )=\left\{\omega \in L^2(\Omega ):D^{\alpha }w\in L^2;forall\left|\alpha \right|\le p\right\}.& (1)\end{array}$

Definition 2. Yadeta et al. [29]. The Caputo derivative of fractional order p with n − 1 < p ≤ n; n ∈ N, for an integrable function g ∈ Cⁿ, can be presented as

$\begin{array}{ll}\begin{array}{c}C\\ a\end{array}{D}_t^{p}g\left(t\right)=\frac{1}{\Gamma \left(n-p\right)}{\displaystyle {\int }_a^t{g}^n\left(\gamma \right){\left(t-\gamma \right)}^{n-p-1}d\gamma }.& (2)\end{array}$

Definition 3. Atangana and Baleanu [21]. Let g ∈ H¹(a, b), a < b, γ ∈ [0, 1], therefore, Atangana–Baleanu–Caputo (ABC) fractional derivative of g(t) with order γ is given by

$\begin{array}{ll}\begin{array}{c}ABC\\ a_1\end{array}{D}_T^{\gamma }\left(g\left(t\right)\right)=\frac{M\left(\gamma \right)}{1-\gamma }{\displaystyle {\int }_{a_1}^t{g}^{\prime }\left(\tau \right)}{E}_{\gamma }\left[-\gamma \frac{{\left(t-\tau \right)}^{\gamma }}{1-\gamma }\right]d\tau , & (3)\end{array}$

where M(γ) is positive and is a normalization function fulfilling M(0) = M(1) = 1, and E_γ is the Mittag-Leffler function

$\begin{array}{ll}{E}_{\gamma }(z)={\displaystyle \sum _{k=0}^{\infty }}\frac{z^k}{\Gamma (\gamma k+1)},& (4)\end{array}$

where γ > 0 and z is complex number.

The Mittag-Leffler function with two parameters appears most frequently and has the following form:

$\begin{array}{ll}{E}_{\alpha ,\beta }(z)={\displaystyle \sum _{k=0}^{\infty }}\frac{z^k}{\Gamma (\alpha k+\beta )}.& (5)\end{array}$

Definition 4 Atangana and Baleanu [21]. Let g ∈ H¹(a, b), a < b, γ ∈ [0, 1], then the Atangana–Baleanu Caputo integral of a function g(t) of order γ is defined by

$\begin{array}{c}\begin{array}{c}ABC\\ a_1\end{array}{I}_T^{\gamma }\left\{g\left(t\right)\right\}=\frac{1-\gamma }{M\left(\gamma \right)}g\left(t\right)+\frac{\gamma }{M\left(\gamma \right)\Gamma \left(\gamma \right)}\\ {\displaystyle {\int }_{a_1}^{t}g\left(\tau \right){\left(t-\tau \right)}^{\gamma -1}}d\tau . & (6)\end{array}$

Definition 5. Atangana and Baleanu [21]. The Laplace transform of the Atangana–Baleanu fractional derivative of order α in Caputo sense is given by

$\begin{array}{ll}L\left\{\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }f\left(t\right)\right\}\left(s\right)=\frac{B\left(\alpha \right)\left(s^{\alpha }F\left(s\right)-s^{\alpha -1}f\left(0\right)\right)}{s^{\alpha }\left(1-\alpha \right)+\alpha },s>0, & (7)\end{array}$

where $L$ is the Laplace transform operator.

Theorem 1. Atangana and Baleanu [21]. For a function g ∈ C[a, b], the inequality shown below holds.

$\begin{array}{ll}\left|{|}_{a_1}^{ABC}{D}_t^{\alpha }g(t)\right||<\frac{1-\alpha }{M(\alpha )}|\left|g(t)\right||,& (8)\end{array}$

where $\Vert g\left(t\right)\Vert =\underset{t\in [a,b]}{\max }\left|g\left(t\right)\right|$.

In addition, ABC derivatives satisfy the Lipschitz requirement; thus, we have

$\begin{array}{ll}\left|{|}_{a_1}^{ABC}{D}_t^{\alpha }{g}_1(t){-}_{a_1}^{ABC}{D}_t^{\alpha }{g}_2(t)\right||<k||g_1(t)-g_2(t)||.& (9)\end{array}$

3 Model formulation

3.1 Presumptions of the model

The human population at time t, denoted by (N_h(t)), comprises two categories: adults and people under the age of seventeen. The latter group is further divided into susceptible children(S_c(t)), exposed children (E_c(t)), and infected children (I_c(t)). Adults are categorized into four compartments: susceptible (S_a(t)), exposed (E_a(t)), infected (I_a(t)), and recovered (W_a(t)). This leads to the following equation:

$\begin{array}{c}{N}_h(t)=S_c(t)+E_c(t)+I_c(t)+S_a(t)+E_a(t)\\ +I_a(t)+W_a(t).& (10)\end{array}$

The total mosquito population at time t, denoted by (N_v(t)), is divided into immature (M_A(t)) and mature mosquitoes (M_v(t)). Furthermore, there are three compartments in mature mosquitoes: susceptible mosquitoes (S_v(t)), exposed mosquitoes (E_v(t)), and infected mosquitoes (I_v(t)). So that:

$\begin{array}{ll}{N}_v(t)=M_A(t)+S_v(t)+E_v(t)+I_v(t).& (11)\end{array}$

In the compartmental model (12), π_c is the recruitment rate for children, λ_i(T)(i = c, a, v) are the infection rate of susceptible children, adults, and vectors, respectively. b_m(T) is the temperature-dependent per capital biting rate of mosquitoes, and β_j(j = c, a, m) is the probability of infection per bite for children, adults, and mosquitoes. μ_h is the natural death rate of humans, and ξ_h is the maturation rate of children to adulthood. Immature mosquitoes (eggs, larvae, and pupae) are lumped into a single compartment (M_A(t)) for computational convenience [see [8] and the references therein].

The temperature-dependent egg deposition rate is represented by α_I(T). Following [14], we assume that the immature mosquito population (encompassing larvae and pupae) is limited by the carrying capacity K_v, where K_v exceeds M_A. This parameter reflects the available nutrients and space, as detailed in [8] and its references. Therefore, ${\alpha }_I(T)(1-\frac{M_A}{K_v})(S_v+E_v+I_v)$ represents the logistic growth rate for the immature mosquitoes.

Using the parameters and their definitions given in Tables 1, 2, we can write the model as a system of first-order ordinary differential equations:

$\begin{array}{ll}\begin{array}{c}\begin{array}{c}\frac{d{S}_c\left(t\right)}{dt}={\pi }_c-\left({\lambda }_c\left(T\right)+{\xi }_h+{\mu }_h\right)S_c\left(t\right),\\ \frac{d{E}_c\left(t\right)}{dt}={\lambda }_c\left(T\right)S_c\left(t\right)-\left({\delta }_c+{\mu }_h\right)E_c\left(t\right),\\ \frac{d{I}_c\left(t\right)}{dt}={\delta }_c{E}_c\left(t\right)-\left({\mu }_h+d_1+{\sigma }_c\right)I_c\left(t\right),\end{array}\\ \frac{d{S}_a\left(t\right)}{dt}={\xi }_h{S}_c+\gamma W_a-\left({\lambda }_a\left(T\right)+{\mu }_h\right)S_a,\\ \frac{d{E}_a\left(t\right)}{dt}={\lambda }_a\left(T\right)S_a-\left({\mu }_h+{\delta }_a\right)E_a,\\ \frac{d{I}_a\left(t\right)}{dt}={\delta }_a{E}_a-\left({\sigma }_a+{\mu }_h+d_2\right)I_a,\\ \frac{d{W}_a\left(t\right)}{dt}={\sigma }_c{I}_c+{\sigma }_a{I}_a-\left({\mu }_h+\gamma \right)W_{a,}\\ \frac{\begin{array}{c}d{M}_A\left(t\right)\end{array}}{dt}={\alpha }_I\left(T\right)\left(1-\frac{M_A}{K_v}\right)\left(S_v+E_v+I_v\right)\\ -\left(\eta \left(T,R\right)+{\mu }_A\left(T\right)\right)M_A,\\ \frac{d{S}_v\left(t\right)}{dt}=\eta \left(T,R\right)M_A-\left({\lambda }_v\left(T\right)+{\mu }_v\left(T\right)\right)S_v,\\ \frac{d{E}_v\left(t\right)}{dt}={\lambda }_v\left(T\right)S_v-\left({\mu }_v\left(T\right)+{\alpha }_v\right)E_v,\\ \frac{d{I}_v\left(t\right)}{dt}={\alpha }_v{E}_v-{\mu }_v\left(T\right)I_V,\end{array}\}& (12)\end{array}$

where the infection rates for children (λ_c(T)), adults (λ_a(T)), and vectors (λ_v (T)) are given, respectively, by

$\begin{array}{r}{\lambda }_c(T)=\frac{b_m(T){\beta }_c{I}_v}{N_h},{\lambda }_a(T)=\frac{b_m(T){\beta }_a{I}_v}{N_h},\\ {\lambda }_v(T)=\frac{b_m(T){\beta }_m(I_a+I_c)}{N_h}.& (13)\end{array}$

Table 1

Table 1. Description of variables of the model (12).

Table 2

Table 2. Description of parameter of the model (12).

4 Model analysis

To account for long-term dependencies and complex memory effects in biological systems, we propose a fractional-order model for malaria. We replace the first-order ordinary derivatives in Equation 12—representing population change rates—with the Atangana–Baleanu fractional derivative of order α (0 < α ≤ 1) to incorporate internal memory effects.

$\begin{array}{cc}\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{c}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{S}_c\left(t\right)={\pi }_c-\left({\lambda }_c\left(T\right)+{\xi }_h+{\mu }_h\right)S_c\left(t\right),\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_c\left(t\right)={\lambda }_c\left(T\right)S_c\left(t\right)-\left({\delta }_c+{\mu }_h\right)E_c\left(t\right),\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_c\left(t\right)={\delta }_c{E}_c\left(t\right)-\left({\mu }_h+d_1+{\sigma }_c\right)I_c\left(t\right),\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{S}_a\left(t\right)={\xi }_h{S}_c+\gamma W_a-\left({\lambda }_a\left(T\right)+{\mu }_h\right)S_a,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_a\left(t\right)={\lambda }_a\left(T\right)S_a-\left({\mu }_h+{\delta }_a\right)E_a,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_a\left(t\right)={\delta }_a{E}_a-\left({\sigma }_a+{\mu }_h+d_2\right)I_a,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{W}_a\left(t\right)={\sigma }_c{I}_c+{\sigma }_a{I}_a-\left({\mu }_h+\gamma \right)W_{a,}\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{M}_A\left(t\right)={\alpha }_I\left(T\right)\left(1-\frac{M_A}{K_v}\right)\left(S_v+E_v+I_v\right)\\ -\left(\eta \left(T,R\right)+{\mu }_A\left(T\right)\right)M_A,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{S}_v\left(t\right)=\eta \left(T,R\right)M_A-\left({\lambda }_v\left(T\right)+{\mu }_v\left(T\right)\right)S_v,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_v\left(t\right)={\lambda }_v{S}_v-\left({\mu }_v\left(T\right)+{\alpha }_v\right)E_v,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_v\left(t\right)={\alpha }_v{E}_v-{\mu }_v\left(T\right)I_V,\end{array}\hfill \end{array}\}\hfill \end{array}& (14)\end{array}$

with initial conditions:

$\begin{array}{r}{S}_c(0)>0,S_a(0)>0,E_c(0)>0,E_a(0)>0,\\ I_c(0)>0,I_a(0)>0,W_a(0)>0,M_A(0)>0,\\ S_v(0)>0,E_v(0)>0,I_v(0)>0,& (15)\end{array}$

where $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }$ is Atangana-Beleanu in Caputo sense of order α, and the infection rates for children (λ_c(T)), adults (λ_a(T)), and vectors (λ_m (T)) are given, respectively, by

$\begin{array}{r}{\lambda }_c(T)=\frac{b_m(T){\beta }_c{I}_v}{N_h},{\lambda }_a(T)=\frac{b_m(T){\beta }_a{I}_v}{N_h},\\ {\lambda }_v(T)=\frac{b_m(T){\beta }_m(I_a+I_c)}{N_h}.& (16)\end{array}$

4.1 Existence and uniqueness of solutions

In this section, we will discuss the existence and uniqueness of solutions ABC fractional operators of the model system (13). To do this, we first re-write the model (13) in the following simple form:

$\begin{array}{ll}\begin{array}{c}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }g\left(t\right)=F\left(t,g\left(t\right)\right),\\ g\left(0\right)=g_0\ge 0.\end{array}\}& (17)\end{array}$

In (17),

$\begin{array}{ll}g\left(t\right)={\left(\begin{array}{c}{S}_c(t,S_a(t),E_c(t),E_a(t),I_c(t),I_a(t),I_c(t),\\ W_a(t),M_A(t),S_v(t),E_v(t),I_v(t)\end{array}\right)}^T\in R_{+}^{11},& (18)\end{array}$

and F is a continuous vector function defined as

$\begin{array}{l}F(t,g(t))=\left(\begin{array}{c}\begin{array}{c}{F}_1(t,g(t))\\ F_2(t,g(t))\\ F_3(t,g(t))\end{array}\\ F_4(t,g(t))\\ F_5(t,g(t))\\ F_6(t,g(t))\\ F_7(t,g(t))\\ F_8(t,g(t))\\ F_9(t,g(t))\\ F_{10}(t,g(t))\\ F_{11}(t,g(t))\end{array}\right)\\ =\left(\begin{array}{c}\begin{array}{c}{\pi }_c-({{\lambda }_c(T)+\xi }_h+{\mu }_h)S_c(t),\\ {\lambda }_c(T)S_c(t)-({\delta }_c+{\mu }_h)E_c(t),\\ {\delta }_c{E}_c(t)-({\mu }_h+d_1+{\sigma }_c)I_c(t),\end{array}\\ {\xi }_h{S}_c+\gamma W_a-({\lambda }_a(T)+{\mu }_h)S_a,\\ {\lambda }_a(T)S_a-({\mu }_h+{\delta }_a)E_a,\\ {\delta }_a{E}_a-({\sigma }_a+{\mu }_h+d_2)I_a,\\ {\sigma }_c{I}_c+{\sigma }_a{I}_a-({\mu }_h+\gamma )W_{a,}\\ {\alpha }_I(T)(1-\frac{M_A}{K_v})(S_v+E_v+I_v)\\ -(\eta (T,R)+{\mu }_A(T))M_A\\ \eta (T,R)M_A-{({\lambda }_v(T)+{\mu }_v(T))S}_v,\\ {\lambda }_v(T)S_v-({\mu }_v(T)+{\alpha }_v)E_v,\\ {\alpha }_v{E}_v-{\mu }_v(T)I_V,\end{array}\right),& (19)\end{array}$

and

$\begin{array}{ll}{g}_0={\left(\begin{array}{c}{S}_c(0),S_a(0),E_c(0),E_a(0),I_c(0),I_a(0),I_c(0),\\ W_a(0),M_A(0),S_v(0),E_v(0),I_v(0)\end{array}\right)}^T& (20)\end{array}$

is the initial vector for state variables. The function F also satisfies the Lipschitz criterion, which is as follows:

$\begin{array}{ll}\left|\right|F(t,g_1(t))-F(t,g_2(t))\left|\right|\le L\left|\right|g_1(t)-g_2(t)\left|\right|,& (21)\end{array}$

where L is a Lipschitz constant.

The following lemma helps us to show the existence and uniqueness of the solution to the fractional model system (14):

Lemma 1. The unique solution of the differential equation with fractional order α is given by

$\begin{array}{ll}\{\begin{array}{c}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }g\left(t\right)=F\left(t,g\left(t\right)\right),\\ g\left(0\right)=g_0\ge 0,\end{array}wherefor\alpha \in (0,1],& (22)\end{array}$

is expressed as

$\begin{array}{r}g(t)=g(0)+\frac{1-\alpha }{B(\alpha )}F(t,g(t))\\ +\frac{\alpha }{B(\alpha )\Gamma (\alpha )}{\displaystyle {\int }_0^t}{F(\tau ,g(\tau ))(t-\tau )}^{\alpha -1}d\tau .& (23)\end{array}$

Theorem 2. Under the condition that $L=(\frac{(1-\alpha )K}{B(\alpha )}+$ $\frac{\alpha K}{B(\alpha )\Gamma (\alpha )}{T}_{max}^{\alpha })<1$, then there exists a unique solution to the fractional-order malaria model (14).

Proof.

Let τ = (0, T), and the operator $F:(\tau ,R^{11})\to (\tau ,R^{11})$, as expressed in lemma (23), is given by

$\begin{array}{r}F(g(t))=g(0)+\frac{1-\alpha }{B(\alpha )}F(t,g(t))\\ +\frac{\alpha }{B(\alpha )\Gamma (\alpha )}{\displaystyle {\int }_0^t}{F(\tau ,g(\tau ))(t-\tau )}^{\alpha -1}d\tau & (24)\end{array}$

$g(t)=F(g(t))$ is expressed as

$\begin{array}{l}g(t)=g(0)+\frac{1-\alpha }{B(\alpha )}F(t,g(t))\\ +\frac{\alpha }{B(\alpha )\Gamma (\alpha )}{\displaystyle {\int }_0^t}{F(\tau ,g(\tau ))(t-\tau )}^{\alpha -1}d\tau & (25)\end{array}$

Furthermore, let:

$\begin{array}{ll}\left|\right|g(t)|{|}_{\tau }={\displaystyle \underset{t\in \tau }{sup}}|g(t)|,g(t)\in C,& (26)\end{array}$

where ||.||_τ represents the supremum norm over τ.

Using Equation 25, we obtain:

$\begin{array}{r}\left|\right|F(g_1(t))-F(g_2(t))\left|\right|=\left|\right|\frac{1-\alpha }{B(\alpha )}(F(g_1(t))-F(g_2(t)))\\ +\frac{\alpha }{B(\alpha )\Gamma (\alpha )}{\displaystyle {\int }_0^t}{(F(\tau ,g_1(\tau ))-F(\tau ,g_2(\tau )))(t-\tau )}^{\alpha -1}d\tau \left|\right|& (27)\end{array}$

From:

$\begin{array}{ll}\left|\right|{\displaystyle {\int }_0^t}D(t,x)u(x)dx\left|\right|\le T\left|\right|D(t,x)\left|{|}_{\tau }\right|\left|u(x)\right|{|}_{\tau }& (28)\end{array}$

With:

u(x) ∈ C(τ, R¹¹) and D(t, τ) ∈ C(τ², R), that is:

$\begin{array}{ll}\left|\right|D(t,x)\left|\right|={\displaystyle \underset{t,x\in \tau }{sup}}\left|D(t,x)\right|& (29)\end{array}$

Now, using Equations 28, 29 and the triangular inequality on Equation 27, we obtain

$\begin{array}{r}\left|\right|F(g_1(t))-F(g_2(t))\left|\right|\le \frac{(1-\alpha )K}{B(\alpha )}\left|\right|g_1(t)-g_2(t)|{|}_{\tau }\\ +\frac{\alpha K}{B(\alpha )\Gamma (\alpha )}{\displaystyle {\int }_0^t}\left|\right|g_1(t)-g_2(t)|{|}_{\tau }{(t-\tau )}^{\alpha -1}d\tau & (30)\end{array}$

$\begin{array}{r}\left|\right|F(g_1(t))-F(g_2(t))\left|\right|\le \frac{(1-\alpha )K}{B(\alpha )}{\displaystyle \underset{t\in \tau }{sup}}|g_1(t)-g_2(t)|\\ +\frac{\alpha K}{B(\alpha )\Gamma (\alpha )}{\displaystyle \underset{t\in \tau }{sup}}|g_1(t)-g_2(t)|{\displaystyle {\int }_0^t}{(t-\tau )}^{\alpha -1}d\tau & (31)\end{array}$

$\begin{array}{r}\left|\right|F(g_1(t))-F(g_2(t))\left|\right|\le (\frac{(1-\alpha )K}{B(\alpha )}+\frac{\alpha K}{B(\alpha )\Gamma (\alpha )}{T}_{max}^{\alpha })\\ \left|\right|g_1(t)-g_2(t)|{|}_{\tau }& (32)\end{array}$

Therefore, we obtain

$\begin{array}{ll}\left|\right|F(g_1(t))-F(g_2(t))\left|\right|\le L\left|\right|g_1(t)-g_2(t)|{|}_{\tau },& (33)\end{array}$

where

$\begin{array}{ll}L=\left(\frac{(1-\alpha )K}{B(\alpha )}+\frac{\alpha K}{B(\alpha )\Gamma (\alpha )}{T}_{max}^{\alpha }\right).& (34)\end{array}$

Hence, if L < 1 on $(\tau ,R_{+}^{11})$, the operator $F$ becomes a contraction. Therefore, the fractional malaria model (14) has a unique solution by the Banach fixed-point theorem.

4.2 Positivity and boundedness of solutions

For the Atangana–Baleanu–Caputo fractional derivative models to be stable, have a steady state of existence, and have biological significance, their solutions have to be positive. Here is a lemma that can be used to prove the positivity of the solution to the proposed model.

Lemma 2. (Generalized Mean Value Theorem see [30]). Supposing that g(t) ∈ C[a, b] and $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }g(t)\in C\left[a,b\right]$ for 0 < α ≤ 1, then

$\begin{array}{c}g(t)=g(k)+{\frac{1}{\Gamma (\alpha )}}_0^{ABC}{D}_t^{\alpha }g(\tau ){(t-k)}^{\alpha },\\ with0\le \tau \le t,\forall t\in \left[a,b\right].& (35)\end{array}$

Remark 1 Suppose that g(t) ∈ C[0, b] and $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }g(t)\in C\left[a,b\right]$ for 0 < α ≤ 1 from Lemma 2 one can deduce that

i. if $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }g(t)\ge 0$, ∀t ∈ (0, b], then the function g(t) is non-decreasing and

ii. if $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }g(t)\le 0$, ∀t ∈ (0, b], then the function g(t) is non-increasing.

$\begin{array}{l}\text{Let }\Omega =\{(S_c(t),S_a(t),E_c(t),E_a(t),I_c(t),I_a(t),\\ {W_a(t),M_A(t),S_v(t),E_v(t),I_v(t))}^T\in R_{+}^{11}\}. & (36)\end{array}$

Theorem 3. The solution of the fractional malaria model (14)

$\begin{array}{ll}x(t)=\left(\begin{array}{c}{S}_c(t),S_a(t),E_c(t),E_a(t),I_c(t),I_a(t),\\ W_a(t),M_A(t),S_v(t),E_v(t),I_v(t)\end{array}\right),& (37)\end{array}$

which starts at, t = t₀ remains in Ω for any t ≥ t₀.

Proof. With reference [31], beginning with the first equation in Equation 14, we want to show that S_c(t) ≥ 0 for all t > t₀.

If, to the contrary, it is not true, then there exists a constant t > t₀ such that

$\begin{array}{ll}\{\begin{array}{c}{S}_c(t)>0,fort\in (t,t_0),\\ S_c(t_{1,})=0,\\ S_c(t_1^{+})<0.\end{array}& (38)\end{array}$

Now, from the first equation in Equation 14 and based on Equation 38, we obtain

$\begin{array}{ll}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{S_c\left(t_1\right)|}_{S_c\left(t_1\right)=0}={\pi }_c>0. & (39)\end{array}$

Consequently, invoking Remark 1 of lemma 2, we obtain $S_c(t_1^{+})\ge 0$, which contracts the assertion $S_c(t_1^{+})<0$. Therefore, we conclude that S_c(t) ≥ 0 for all t > t₀. Similarly, we establish the following inequalities for t > t₀:

$\begin{array}{r}{E}_c(t)\ge 0,I_c(t)\ge 0,S_a(t)\ge 0,E_a(t)\ge 0,\\ I_a(t)\ge 0,W_a(t)\ge 0,M_A(t)\ge 0,S_v(t)\ge 0,\\ E_v(t)\ge 0,I_v(t)\ge 0.\end{array}$

This completes the proofs of Theorem 3.

4.3 Invariant region

Given non-negative initial data, we identify the region where the solutions to the system of Equation 14 are viable. Here, we are primarily concerned with demonstrating that the viable region is located in $R_{+}^{11}$, which is a positive invariant region with regard to the model (14) under the initial condition (15).

Theorem 4. Let the Atangana–Baleanu fractional model (14) has a unique solution (N_h, N_v) for all t ≥ 0. Then, the epidemiologically feasible region of the model is given by $\Omega ={\Omega }_h\times {\Omega }_v\subset R_{+}^{11}$, where

$\begin{array}{r}{\Omega }_h=\{\left(S_c\left(t\right),S_a\left(t\right),E_c\left(t\right),E_a\left(t\right),I_c\left(t\right),I_a\left(t\right),W_A(t)\right)\\ \in R_{+}^7:0\le N_h\le \frac{{\pi }_c}{{\mu }_h}\}& (40)\end{array}$

and

$\begin{array}{r}{\Omega }_v=\{\left(M_A(t),S_v(t),E_v(t),I_v(t)\right)\\ \in R_{+}^4:0\le M_v\le \frac{\eta \left(T,R\right)M_A}{{\mu }_v\left(T\right)}\}.& (41)\end{array}$

Proof. Let N_h and N_v represent the total population of humans and mosquitoes, respectively. By adding all the equations corresponding to the human and mosquito components of the system (14), we get

$\begin{array}{ll}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(N_h\left(t\right)\right)={\pi }_c-{\mu }_h{N}_h-\left(d_1{I}_c+d_2{I}_a\right)\le {\pi }_c-{\mu }_h{N}_h & (42)\end{array}$

and

$\begin{array}{ll}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(M_v\left(t\right)\right)=\eta \left(T,R\right)M_A-{\mu }_v\left(T\right)N_v & (43)\end{array}$

To prove theorem 4, we need to show that system (14) has bounded solutions. Biologically, the least possible value of each state of the model system (14) is zero. Next, we determine the upper bound of the states. To do this,

let

$\begin{array}{ll}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(N_h\left(t\right)\right)\le {\pi }_c-{\mu }_h{N}_h & (44)\end{array}$

and

$\begin{array}{ll}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(N_v\left(t\right)\right)=\eta \left(T,R\right)M_A-{\mu }_v\left(T\right)N_v. & (45)\end{array}$

Then, by applying the Laplace transform on both sides of Equation 44, we obtain

$\begin{array}{ll}L\left[\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(N_h\left(t\right)\right)\right]\left(s\right)\le \frac{{\pi }_c}{s}-{\mu }_{h}L\left[N_h\left(t\right)\right]\left(s\right), & (46)\end{array}$

and using definition 5, we obtain

$\begin{array}{ll}\frac{B(\alpha )(s^{\alpha }L\left[N_h(t)\right](s)-s^{\alpha -1}{N}_h(0))}{s^{\alpha }(1-\alpha )+\alpha }\le \frac{{\pi }_c}{s}+\frac{s^{\alpha -1}{N}_h(0)}{s^{\alpha }(1-\alpha )+\alpha },& (47)\end{array}$

where N_h(0) represents the initial value of the total human population.

This implies that

$\begin{array}{r}L\left\{N_h(t)\right\}(s)\le (\frac{s^{\alpha -1}{N}_h(0)B(\alpha )}{B(\alpha )s^{\alpha }+{\mu }_H(s^{\alpha }(1-\alpha )+\alpha )})\\ +\frac{{\pi }_c}{s}(\frac{s^{\alpha }(1-\alpha )+\alpha }{B(\alpha )s^{\alpha }+{\mu }_h(s^{\alpha }(1-\alpha )+\alpha )})& (48)\end{array}$

Therefore,

$\begin{array}{r}L\left\{N_h(t)\right\}(s)\le \frac{{\pi }_c\alpha }{(B(\alpha )+{\mu }_h(1-\alpha ))}(\frac{s^{\alpha -(\alpha +1)}}{s^{\alpha }+\frac{{\mu }_h\alpha }{B(\alpha )+{\mu }_h(1-\alpha )}})\\ +(\frac{{\pi }_c(1-\alpha )}{B(\alpha )+{\mu }_h(1-\alpha )}+\frac{N_h(0)B(\alpha )}{B(\alpha )+{\mu }_h(1-\alpha )})\\ \frac{s^{\alpha -1}}{s^{\alpha }+\frac{{\mu }_h\alpha }{B(\alpha )+{\mu }_h(1-\alpha )}}.& (49)\end{array}$

Applying the inverse Laplace transform on both sides of (49), we get

$\begin{array}{l}{N}_h(t)\le \frac{{\pi }_c\alpha t^{\alpha }}{(B(\alpha )+{\mu }_h(1-\alpha ))}{E}_{\alpha ,\alpha +1}(-k{t}^{\alpha })\\ +(\frac{{\pi }_c(1-\alpha )}{B(\alpha )+{\mu }_h(1-\alpha )}+\frac{N_h(0)B(\alpha )}{B(\alpha )+{\mu }_h(1-\alpha )})E_{\alpha ,1}(-k{t}^{\alpha }),& (50)\end{array}$

where

$\begin{array}{ll}k=\frac{{\mu }_h\alpha }{B(\alpha )+{\mu }_h(1-\alpha )}.& (51)\end{array}$

Furthermore,

$\begin{array}{ll}{E}_{\alpha ,\alpha +1}(-k{t}^{\alpha })=\frac{1}{-k{t}^{\alpha }}(E_{\alpha ,1}(-k{t}^{\alpha })-1).& (52)\end{array}$

Thus,

$\begin{array}{rr}\begin{array}{l}{N}_h\left(t\right)\le \frac{{\pi }_c}{-{\mu }_h}\left(E_{\alpha ,1}\left(-k{t}^{\alpha }\right)-1\right)+(\frac{{\pi }_c\left(1-\alpha \right)}{B\left(\alpha \right)+{\mu }_h\left(1-\alpha \right)}\\ +\frac{N_h\left(0\right)B\left(\alpha \right)}{B\left(\alpha \right)+{\mu }_h\left(1-\alpha \right)})E_{\alpha ,1}\left(-k{t}^{\alpha }\right)=\frac{{\pi }_c}{{\mu }_h} \end{array}& (53)\end{array}$

since $E_{\alpha ,1}(-k{t}^{\alpha })\to 0$ as t → ∞.

Thus, the epidemiologically feasible region for the human population is

$\begin{array}{ll}{\Omega }_h=\left\{(S_c,E_c,I_c,S_a,E_a,I_a,W_a)\in R_{+}^7:0\le N_h\le \frac{{\pi }_c}{{\mu }_h}\right\}.& (54)\end{array}$

Similarly, it can be shown that the feasible region for the mosquito population is

$\begin{array}{ll}{\Omega }_v=\left\{(M_A,S_M,E_M,I_M)\in R_{+}^3:0\le S_v+E_v+I_v\le \frac{\eta (T,R)M_A}{{\mu }_v(T)}\right\}.& (55)\end{array}$

Thus, the proposed model is mathematically well-posed and epidemiologically sound on the region $\Omega ={\Omega }_h\times {\Omega }_v\subset R_{+}^{11}$.

4.4 Existence of equilibrium points and the basic reproduction number

This section presents the disease-free and endemic equilibrium points of our Atangana–Baleanu fractional-order malaria model (14) to analyze their stability and dynamical behavior.

4.4.1 Disease-free equilibrium point

The steady-state solution of the model system (a form of the ABC fractional model) (14) obtained in the absence of disease is known as the disease-free equilibrium (DFE). The autonomous form of Equation 14 exhibits two disease-free equilibria, as stated in Theorem 6 below, depending on the magnitude of the threshold quantity M, where $M=\frac{{\eta }_I{\alpha }_I}{({\eta }_I+{\mu }_A){\mu }_v}$. For this analysis, we consider the special case of the non-autonomous version of model (14) in which the temperature- and rainfall-dependent parameters are constant, specifically: η_I(T, R) = η_I, α_I(T) = α_I, μ_A(T) = μ_A, and μ_v(T) = μ_v.

Theorem 5. The model system (14) possesses two malaria disease-free equilibrium points: the trivial disease-free equilibrium $(E_1^0)$ and the biologically realistic disease-free equilibrium $(E_2^0)$. These equilibria are defined as follows:

• $E_1^0=(S_c^0,0,0,0,0,0,0,0,0,0,0)$ holds when $M\le 1$.

• $E_2^0=(S_c^0,0,0,0,0,0,M_A^0,S_v^0,0,0,0)$ holds when $M>1$.

where

$\begin{array}{ll}{S}_c^0=\frac{{\pi }_c}{{\xi }_h+{\mu }_h}& (56)\end{array}$

$\begin{array}{ll}{M}_A^0=(1-\frac{1}{M})K_v& (57)\end{array}$

$\begin{array}{ll}{S}_v^0=\frac{{\eta }_I}{{\mu }_v}\left[1-\frac{1}{M}\right]K_v& (58)\end{array}$

Proof: To find the disease-free equilibria of the ABC fractional malaria model (14), we set the right-hand side of each equation in system (14) to zero and also set E_c = I_c = E_a = I_a = E_v = I_v = 0. This means all infected and infectious compartments are empty.

From the equations $\begin{array}{c}ABC\\ 0\end{array}D{M}_A(t)=\begin{array}{c}ABC\\ 0\end{array}D{S}_v(t)=0$, we obtain two possibilities: either M_A = 0 or $M_A=K_v(1-\frac{1}{M})$.

Therefore:

i. $E_1^0$: When M_A = 0 and all infected and infectious compartments are empty (E_c = I_c = E_a = I_a = E_v = I_v = 0), and each equation in Equation 14 is equal to zero, we obtain

$\begin{array}{ll}{S}_c^0=\frac{{\pi }_c}{{\xi }_h+{\mu }_h}& (59)\end{array}$

$\begin{array}{r}{E}_c^0=I_c^0=S_a^0=E_a^0=I_a^0\\ =W_a^0=M_A^0=S_v^0=E_v^0=I_v^0=0& (60)\end{array}$

Therefore, $E_1^0=(S_c^0,0,0,0,0,0,0,0,0,0,0)$.

ii. $E_2^0$: When $M_A=K_v(1-\frac{1}{M})$ and all infected and infectious compartments are empty (E_c = I_c = E_a = I_a = E_v = I_v = 0), and each equation in Equation 14 is equal to zero, we obtain

$\begin{array}{ll}{S}_c^0=\frac{{\pi }_c}{{\xi }_h+{\mu }_h}& (61)\end{array}$

$\begin{array}{ll}{E}_c^0=I_c^0=S_a^0=E_a^0=I_a^0=W_a^0=E_v^0=I_v^0=0& (62)\end{array}$

$\begin{array}{ll}{S}_v^0=\frac{{\eta }_I}{{\mu }_v}\left[1-\frac{1}{M}\right]K_v& (63)\end{array}$

Therefore, $E_2^0=(S_c^0,0,0,0,0,0,M_A^0,S_v^0,0,0,0)$.

This completes the proof of Theorem 5.

4.4.2 The basic reproduction number

Definition 6. Diekmann et al. [32] and Van den Driessche and Watmough [33]. The basic reproductive number $R_0=\rho (F{V}^{-1})$ is the spectral radius (largest eigenvalue) of the next-generation matrix, where F represents the Jacobian of the rates of flows from uninfected to infected classes evaluated at the disease-free equilibrium and V is the Jacobian of the rates of all other flows to and from infected classes evaluated at the disease-free equilibrium.

Theorem 6.

i. If $M\le 1$, then the basic reproduction number, R₀, associated with system (14) is zero.

ii. If $M>1$, then the basic reproduction number, R₀, associated with system (14) is:

$\begin{array}{c}{R}_0\\ =\sqrt{\frac{{(b_m)}^2{({\delta }_c)}^2{\beta }_m{\alpha }_v{\pi }_c\eta }{({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)({\mu }_v+{\alpha }_v){(N_h)}^2{({\mu }_v)}^2({\xi }_h+{\mu }_h)}\left[1-\frac{1}{M}\right]k_v},\\ \text{where}\\ M=\frac{\eta {\alpha }_I}{(\eta +{\mu }_A){\mu }_v}.\end{array}$

Proof. To compute the basic reproduction number of the ABC fractional malaria model (14) from the malaria disease-free equilibrium using the next-generation matrix method [33], we focus on the infected compartments (those representing disease progression dynamics) of the model, leading to the following subsystem:

$\begin{array}{cc}\begin{array}{c}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_c=\frac{b_m\left(T\right){\beta }_c{I}_v}{N_h}{S}_c-\left({\delta }_c+{\mu }_h\right)E_c,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_c={\delta }_c{E}_c-\left({\mu }_h+d_1+{\sigma }_c\right)I_c\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_a=\frac{b_m\left(T\right){\beta }_a{I}_v}{N_h}{S}_a-\left({\mu }_h+{\delta }_a\right)E_a,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_a={\delta }_a{E}_a-\left({\sigma }_a+{\mu }_h+d_2\right)I_a,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_v=\frac{b_m{\beta }_m\left(I_a+I_c\right)}{N_h}{S}_v-\left({\mu }_v+{\alpha }_v\right)E_v,\\ \begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_v\left(t\right)={\alpha }_v{E}_v-{\mu }_v{I}_V.\end{array}\}& (64)\end{array}$

Thus, the Jacobian matrix for the infected sub-population Equation 64 is given by

$\begin{array}{c}J(E_c,I_c,E_a,I_a,E_v,I_v)\\ =\left[\begin{array}{cccccc}-({\delta }_c+{\mu }_h)& 0& 0& 0& 0& \frac{{S_{c}b}_m{\beta }_c}{N_h}\\ {\delta }_c& -({\mu }_h+d_1+{\sigma }_c)& 0& 0& 0& 0\\ 0& 0& -({\mu }_h+{\delta }_a)& 0& 0& \frac{{S_{a}b}_m{\beta }_a}{N_h}\\ 0& 0& {\delta }_a& -({\sigma }_a+{\mu }_h+d_2)& 0& 0\\ 0& \frac{{S_{v}b}_m{\beta }_m}{N_h}& 0& \frac{{S_{v}b}_m{\beta }_m}{N_h}& -({\mu }_v+{\alpha }_v)& 0\\ 0& 0& 0& 0& {\alpha }_v& -{\mu }_v\end{array}\right]& (65)\end{array}$

Now, we look at two cases to determine the basic production number.

Case I: Following (65), the transmission matrix F (of new infection terms) and the transition matrix V (of transition terms) corresponding to the fractional ABC model (14) at trivial disease-free equilibrium point ($E_1^0$) are given, respectively, by

$\begin{array}{ll}F=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& \frac{{S_c^{0}b}_m{\beta }_c}{N_h}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]& (66)\end{array}$

and

$\begin{array}{ll}V=\left[\begin{array}{cccccc}({\delta }_c+{\mu }_h)& 0& 0& 0& 0& 0\\ -{\delta }_c& ({\mu }_h+d_1+{\sigma }_c)& 0& 0& 0& 0\\ 0& 0& ({\mu }_h+{\delta }_a)& 0& 0& 0\\ 0& 0& -{\delta }_a& ({\sigma }_a+{\mu }_h+d_2)& 0& 0\\ 0& 0& 0& 0& ({\mu }_v+{\alpha }_v)& 0\\ 0& 0& 0& 0& {-\alpha }_v& {\mu }_v\end{array}\right],& (67)\end{array}$

where $S_c^0=\frac{{\pi }_c}{{\xi }_h+{\mu }_h}$.

$\begin{array}{c}\Rightarrow V^{-1}\\ =\left[\begin{array}{cccccc}\frac{1}{({\delta }_c+{\mu }_h)}& 0& 0& 0& 0& 0\\ \frac{{\delta }_c}{({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)}& \frac{1}{({\mu }_h+d_1+{\sigma }_c)}& 0& 0& 0& 0\\ 0& 0& \frac{1}{({\mu }_h+{\delta }_a)}& 0& 0& 0\\ 0& 0& \frac{{\beta }_a}{({\mu }_h+{\delta }_a)({\sigma }_a+{\mu }_h+d_2)}& \frac{1}{({\sigma }_a+{\mu }_h+d_2)}& 0& 0\\ 0& 0& 0& 0& \frac{1}{({\mu }_v+{\alpha }_v)}& 0\\ 0& 0& 0& 0& \frac{{\alpha }_v}{({\mu }_v+{\alpha }_v){\mu }_v}& \frac{1}{{\mu }_v}\end{array}\right].& (68)\end{array}$

Thus,

$\begin{array}{ll}F{V}^{-1}=\left[\begin{array}{cccccc}0& 0& 0& 0& \frac{{{{-\alpha }_{v}S}_c^{0}b}_m(T){\beta }_c}{N_h^{*}}& \frac{{{{\mu }_{v}S}_c^{0}b}_m(T){\beta }_c}{N_h^{*}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right].& (69)\end{array}$

The basic reproduction number in this particular case is zero because zero is the only eigenvalue of the matrix FV⁻¹. This suggests that the disease will not spread if an infected person is introduced into the community at trivial disease-free site $(E_1^0)$, that is, trivial disease-free site $(E_1^0)$ is locally asymptotically stable for this particular case.

Case II: Following (65), the transmission matrix F (of new infection terms) and the transition matrix V (of transition terms) associated with the fractional ABC model (14) at realistic disease-free equilibrium point $(E_2^0)$ are given, respectively, by

$\begin{array}{ll}F=\left[\begin{array}{cccccc}0& 0& 0& 0& 0& \frac{{S_c^{0}b}_m{\beta }_c}{N_h^{*}}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& \frac{{S_v^{0}b}_m{\beta }_m}{N_h^{*}}& 0& \frac{{S_v^{0}b}_m{\beta }_m}{N_h^{*}}& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]& (70)\end{array}$

and

$\begin{array}{ll}V=\left[\begin{array}{cccccc}({\delta }_c+{\mu }_h)& 0& 0& 0& 0& 0\\ -{\delta }_c& ({\mu }_h+d_1+{\sigma }_c)& 0& 0& 0& 0\\ 0& 0& ({\mu }_h+{\delta }_a)& 0& 0& 0\\ 0& 0& -{\delta }_a& ({\sigma }_a+{\mu }_h+d_2)& 0& 0\\ 0& 0& 0& 0& ({\mu }_v+{\alpha }_v)& 0\\ 0& 0& 0& 0& {-\alpha }_v& {\mu }_v\end{array}\right].& (71)\end{array}$

$\begin{array}{c}\Rightarrow V^{-1}\\ =\left[\begin{array}{cccccc}\frac{1}{({\delta }_c+{\mu }_h)}& 0& 0& 0& 0& 0\\ \frac{{\delta }_c}{({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)}& \frac{1}{({\mu }_h+d_1+{\sigma }_c)}& 0& 0& 0& 0\\ 0& 0& \frac{1}{({\mu }_h+{\delta }_a)}& 0& 0& 0\\ 0& 0& \frac{{\delta }_a}{({\mu }_h+{\delta }_a)({\sigma }_a+{\mu }_h+d_2)}& \frac{1}{({\sigma }_a+{\mu }_h+d_2)}& 0& 0\\ 0& 0& 0& 0& \frac{1}{({\mu }_v+{\alpha }_v)}& 0\\ 0& 0& 0& 0& \frac{{\alpha }_v}{({\mu }_v+{\alpha }_v){\mu }_v}& \frac{1}{{\mu }_v}\end{array}\right].& (72)\end{array}$

As a result, one can get

$\begin{array}{c}\text{where}\\ F{V}^{-1}=\left[\begin{array}{cccccc}0& 0& 0& 0& \frac{{S_c^{0}b}_m{\beta }_c{\alpha }_v}{({\mu }_v+{\alpha }_v)N_h^{*}{\mu }_v}& \frac{{S_c^{0}b}_m{\beta }_c}{N_h^{*}{\mu }_v}\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ \frac{b_m{\delta }_c{S}_v^0{\beta }_m}{({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)N_h^{*}}& \frac{b_m{S}_v^0{\beta }_m}{({\mu }_h+d_1+{\sigma }_c)N_h^{*}}& \frac{b_m{S}_v^0{\beta }_m{\beta }_m}{({\mu }_h+{\delta }_a)({\sigma }_a+{\mu }_h+d_2)N_h^{*}}& \frac{b_m{S}_v^0{\beta }_m}{({\sigma }_a+{\mu }_h+d_2)N_h^{*}}& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]& (73)\end{array}$

From (73), we get the characteristic polynomial

$\begin{array}{ll}0={\lambda }^4\left({\lambda }^{2-}\frac{{(b_m)}^2{{({\delta }_c)}^{2}S}_c^0{S}_v^0{\beta }_m{\alpha }_v}{({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)({\mu }_v+{\alpha }_v){(N_h^{*})}^2{\mu }_v}\right),& (74)\end{array}$

where

$\begin{array}{ll}{S}_c^0=\frac{{\pi }_c}{{\xi }_h+{\mu }_h},N_h^{*}=\frac{{\pi }_c}{{\mu }_h},and{S}_v^0=\frac{{\eta }_I}{{\mu }_v}\left[1-\frac{(\eta +{\mu }_A){\mu }_v}{\eta {\alpha }_I}\right]K_v& (75)\end{array}$

which implies that

$\begin{array}{c}{R}_0\\ =\sqrt{\frac{{(b_m)}^2{({\delta }_c)}^2{\beta }_m{\alpha }_v{\eta }_I{({\mu }_h)}^2}{({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)({\mu }_v+{\alpha }_v){\pi }_c{({\mu }_v)}^2({\xi }_h+{\mu }_h)}\left[1-\frac{({\eta }_I+{\mu }_A){\mu }_v}{{\eta }_I{\alpha }_I}\right]k_v}.& (76)\end{array}$

Therefore,

$\begin{array}{c}{R}_0\\ =\sqrt{\frac{{(b_m)}^2{({\delta }_c)}^2{\beta }_m{\alpha }_v\eta {({\mu }_h)}^2}{({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)({\mu }_v+{\alpha }_v){\pi }_c{({\mu }_v)}^2({\xi }_h+{\mu }_h)}\left[1-\frac{1}{M}\right]k_v},& (77)\end{array}$

where

$\begin{array}{ll}M=\frac{\eta {\alpha }_I}{(\eta +{\mu }_A){\mu }_v}.& (78)\end{array}$

This completes the proof of the Theorem 6.

Theorem 7. The malaria disease-free equilibrium points $E_1^0$ and $E_2^0$ are locally asymptotically stable if R₀ < 1 and unstable if R₀ > 1.

Proof. The Jacobian matrix $J_{E_1^0}$ is found about the malaria disease-free equilibrium$(E_1^0)$:

$\begin{array}{ll}{J}_{E_1^0}=\left(\begin{array}{ccccccccccc}-A_1& 0& 0& 0& 0& 0& 0& 0& 0& 0& -B\\ 0& -A_2& 0& 0& 0& 0& 0& 0& 0& 0& B\\ 0& {\beta }_c& -A_3& 0& 0& 0& 0& 0& 0& 0& 0\\ {\xi }_h& 0& 0& -A_4& 0& 0& \gamma & 0& 0& 0& 0\\ 0& 0& 0& 0& -A_5& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& {\beta }_a& -A_6& 0& 0& 0& 0& 0\\ 0& 0& {\sigma }_c& 0& 0& {\sigma }_a& -A_7& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& -A_8& A_9& A_9& A_9\\ 0& 0& 0& 0& 0& 0& 0& {\eta }_I& -{\mu }_v& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& -A_{10}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& {\alpha }_v& -{\mu }_v\end{array}\right),& (79)\end{array}$

where

$\begin{array}{r}{A}_1=({\xi }_h+{\mu }_h),A_2=({\delta }_c+{\mu }_h),\\ A_3=({\mu }_h+d_1+{\sigma }_c),A_4={\mu }_h,\\ A_5=({\mu }_h+{\delta }_a),A_6=({\sigma }_a+{\mu }_h+d_2),\\ A_7=({\mu }_h+\gamma ),A_8=-{\alpha }_I(T)\frac{1}{K_v}{S}_v^0-(\eta +{\mu }_A),\\ A_9={\alpha }_I(T)(1-\frac{M_A^0}{K_v}),A_{10}=({\mu }_v+{\alpha }_v).& (80)\end{array}$

The eigenvalues of Equation 79 are−(ξ_h + μ_h)−(δ_c + μ_h), −(μ_h + d₁ + σ_c), −μ_h, −(μ_h + δ_a), −(σ_a + μ_h + d₂), −(μ_h + γ), −(μ_v + α_v), −μ_v, and $\frac{-{\mu }_v-(\eta +{\mu }_A)\pm \sqrt{{{\mu }_v}^2+{(\eta +{\mu }_A)}^2-2{\mu }_v(\eta +{\mu }_A)+4{\alpha }_I(T){\eta }_I}}{2}$ have all negative real part, so for R₀ < 1, the disease-free equilibrium is locally asymptotically stable.

The following Jacobian matrix $J_{E_2^0}$ is found about the biologically realistic disease-free equilibrium $(E_2^0)$:

$\begin{array}{l}{J}_{E_2^0}=\\ \left[\begin{array}{ccccccccccc}-A_1& 0& 0& 0& 0& 0& 0& 0& 0& 0& -\frac{b_m{\beta }_c{S}_c^0}{N_h}\\ 0& -A_2& 0& 0& 0& 0& 0& 0& 0& 0& \frac{b_m{\beta }_c{S}_c^0}{N_h}\\ 0& {\beta }_c& -A_3& 0& 0& 0& 0& 0& 0& 0& 0\\ {\xi }_h& 0& 0& -A_4& 0& 0& \gamma & 0& 0& 0& 0\\ 0& 0& 0& 0& -A_5& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& {\beta }_a& -A_6& 0& 0& 0& 0& 0\\ 0& 0& {\sigma }_c& 0& 0& {\sigma }_a& -A_7& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& -A_8& A_9& A_9& A_9\\ 0& 0& \frac{-b_m{\beta }_m{S}_v^0}{N_h}& 0& 0& \frac{-b_m{\beta }_m{S}_v^0}{N_h}& 0& {\eta }_I& -{\mu }_v& 0& 0\\ 0& 0& \frac{b_m{\beta }_m{S}_v^0}{N_h}& 0& 0& \frac{b_m{\beta }_m{S}_v^0}{N_h}& 0& 0& 0& -A_{10}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& {\alpha }_v& -{\mu }_v\end{array}\right],& (81)\end{array}$

where

$\begin{array}{r}{A}_1=({\xi }_h+{\mu }_h),A_2=({\delta }_c+{\mu }_h),\\ A_3=({\mu }_h+d_1+{\sigma }_c),A_4={\mu }_h,A_5=({\mu }_h+{\delta }_a),\\ A_6=({\sigma }_a+{\mu }_h+d_2),A_7=({\mu }_h+\gamma ),\\ A_8=-{\alpha }_I(T)\frac{1}{K_v}{S}_v^0-(\eta +{\mu }_A),\\ A_9={\alpha }_I(T)(1-\frac{M_A^0}{K_v}),A_{10}=({\mu }_v+{\alpha }_v).& (82)\end{array}$

It follows from Equation 81, the eigenvalues −(ξ_h + μ_h), −μ_h, −(μ_h + γ), −(σ_a + μ_h + d₂), −(μ_h + β_a), and ${\lambda }^2+({\mu }_v+{\alpha }_I\frac{1}{K_v}{S}_v^0+(\eta +{\mu }_A))\lambda +({\mu }_v+{\alpha }_I\frac{1}{K_v}{S}_v^0)=0$ containing negative real parts. The remaining (four) eigenvalues are found in the roots of the equation provided below:

$\begin{array}{ll}{p(\lambda )=\lambda }^4+K_1{\lambda }^3+{K_2\lambda }^2+K_3\lambda {+K}_4=0,& (83)\end{array}$

where the coefficients are

$\begin{array}{ll}\{\begin{array}{c}{K}_1=({\delta }_c+{\mu }_h)+({\mu }_h+d_1+{\sigma }_c)+({\mu }_v+{\alpha }_v)+{\mu }_v,\\ K_2=({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)+({\delta }_c+{\mu }_h){\mu }_v\\ +({\mu }_h+d_1+{\sigma }_c){\mu }_v\\ +({\delta }_c+{\mu }_h)({\mu }_v+{\alpha }_v)+({\mu }_h+d_1+{\sigma }_c)({\mu }_v+{\alpha }_v)\\ +({\mu }_v+{\alpha }_v){\mu }_v,\\ K_3=({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c){\mu }_v\\ +({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)({\mu }_v+{\alpha }_v)\\ +({\delta }_c+{\mu }_h)({\mu }_v+{\alpha }_v){\mu }_v+({\mu }_h+d_1+{\sigma }_c)({\mu }_v+{\alpha }_v){\mu }_v,\\ K_4=({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)({\mu }_v)({\mu }_v+{\alpha }_v)\\ -(\frac{b_m{\beta }_c{S}_c^0}{N_h})({\beta }_c)(\frac{b_m{\beta }_m{S}_v^0}{N_h})({\alpha }_v).\end{array}& (84)\end{array}$

Thus, by Routh–Hurwitz criteria, the four eigenvalues found in the roots of ${\lambda }^4+K_1{\lambda }^3+{K_2\lambda }^2+K_3\lambda {+K}_4=0$ will have a negative real part if they satisfy the Routh–Hurwitz criteria, that is, K_i > 0, for i = 1, ⋯  , 4.

It can be easily seen from the first, second, and third equations in Equation 84 that K₁ > 0, K₂ > 0, and K₃ > 0, respectively.

Furthermore, from the fourth equation in Equation 84, we have

$\begin{array}{c}{K}_4=({\delta }_c+{\mu }_h)({\mu }_h+d_1+{\sigma }_c)({\mu }_v)({\mu }_v+{\alpha }_v)(1-R_0^2),\\ so{K}_4>0,& (85)\end{array}$

if R₀ < 1 and K₄ < 0, if R₀ > 1.

Thus, $|arg{\lambda }_j|>\frac{\propto \pi }{2}$ for all 0 < α ≤ 1. Therefore, $E_2^0$ will be locally asymptotically stable for R₀ < 1 and unstable for R₀ > 1 as in [34].

Lemma. Vargas-De-León [35]. Let f(t) ∈ R⁺ be a continuous and differentiable function. Then, for any time instance t ≥ 0.

$\begin{array}{ll}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(f\left(t\right)-f^{\ast }-f^{\ast }\ln \frac{f\left(t\right)}{f^{\ast }}\right)\le \left(1-\frac{f^{\ast }}{f\left(t\right)}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha } & (86)\end{array}$

and

$\begin{array}{ll}{{\frac{1}{2}}_0^{ABC}D}_t^{\alpha }{f}^2(t)\le \alpha {(t)}_0^{ABC}Df(t),& (87)\end{array}$

where 0 < α < 1.

Note that α = 1, the inequalities in Equations 86, 87 become equalities.

Theorem 8. If R₀ < 1, then the malaria-free equilibria, $E_1^0$ and $E_2^0$, of the proposed model (13) are global asymptotic stability.

Proof. To prove this, we construct a candidate Lyapunov function $L:R_{+}^{11}\to R$ [36, 37] such that

$\begin{array}{c}L(S_c,E_c,I_c,S_a,E_a,I_a,W_a,M_A,S_v,E_v,I_v)\\ =({S_c-S}_c^{*}-S_c^{*}ln\frac{S_c}{S_c^{*}})\\ +\begin{array}{c}{E}_c+I_c+S_a+E_a+I_a+W_a+M_A+S_v+E_v+I_v\end{array},& (88)\end{array}$

where $S_c^0=\frac{{\pi }_c}{{\xi }_h+{\mu }_h}$.

Now applying ABC operator on both sides of Equation 88, we get

$\begin{array}{l}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }L=\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(S_c-S_c^{\ast }-S_c^{\ast }\ln \frac{S_c}{S_c^{\ast }}\right)+\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_c+\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_c\\ +\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{S}_a+\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_a+\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_a+\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{W}_a\\ +\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{M}_A+\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{S}_v+\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_v+\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_v& (89)\end{array}$

Substituting Equation 14 into Equation 86 and after some algebraic simplification, we obtain

$\begin{array}{cc}\begin{array}{cc}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }L\le & \begin{array}{l}\frac{-{\pi }_c}{S_c^0{S}_c}{\left(S_c-S_c^0\right)}^2\\ +\left(\frac{b_m\left(T\right){\delta }_c}{N_h}{I}_v+2{\xi }_h+{\mu }_h\right)\\ S_c^0{\pi }_c-{\mu }_h\left(E_c+I_c+S_a+E_a+I_a+W_{a,}\right)\\ -\left(d_1{I}_c+d_2{I}_a\right)+{\alpha }_I\left(T\right)\left(1-\frac{M_A}{K_v}\right)\left(S_v+E_v+I_v\right)\\ -{\mu }_v\left(T\right)M_A-{\mu }_v\left(T\right)\left(S_v+E_v+I_V\right)+{\xi }_h{S}_c\end{array}\end{array}& (90)\end{array}$

As

$\begin{array}{r}\left(\frac{b_m(T){\delta }_c}{N_h}{I}_v+{\xi }_h+{\mu }_h\right)S_c^0{\pi }_c\\ +{\alpha }_I(T)(S_v+E_v+I_v)+{{\xi }_{h}S}_c\ge 0,& (91)\end{array}$

we conclude that

$\begin{array}{ll}\begin{array}{cc}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }L\le & \begin{array}{l}\frac{-{\pi }_c}{S_c^0{S}_c}{\left(S_c-S_c^0\right)}^2-{\mu }_h\left(E_c+I_c+S_a+E_a+I_a+W_{a,}\right)\\ -\left(d_1{I}_c+d_2{I}_a\right)-{\mu }_v\left(T\right)M_A\\ -{\alpha }_I\left(T\right)\frac{M_A}{K_v}\left(S_v+E_v+I_v\right)-{\mu }_v\left(T\right)\left(S_v+E_v+I_V\right)\end{array}\end{array}& (92)\end{array}$

As a result, $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }L\le 0$ when 0 < α < 1 and $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }L=0$ if and only if $S_c=S_c^0,E_c=E_c^0,{I_c^0=I}_c,S_a^0{=S}_a=0,E_a^0=E_a=0,I_a^0=I_a=0,W_a^0=W_{a,}=0,{S_v=S}_v^0,I_c=I_0^c=0,I_a=I_a^0=0,M_A=M_A^0$. Thus, by LaSalle's invariance principle [38], the malaria-free equilibria, $E_1^0$ and $E_2^0$, are globally asymptotically stable.

4.5 Existence of endemic equilibrium

The endemic equilibrium point of the fractional model (14), denoted by

$E^{*}=(S_c^{*},E_c^{*},I_c^{*},S_a^{*},E_a^{*},I_a^{*},W_a^{*},M_A^{*},S_v^{*},E_v^{*},I_v^{*})$, is defined where d₁ = d₂ = 0, which leads to, $N_h=\frac{{\pi }_c}{{\mu }_h}$. This equilibrium point satisfies the following equation:

$\begin{array}{ll}{\lambda }_a^{*}p({\lambda }_a^{*})={\lambda }_a^{*}(D_2{({\lambda }_a^{*})}^2+D_1({\lambda }_a^{*})+D_0)=0.& (93)\end{array}$

Accordingly, the roots of Equation 93 that are either ${\lambda }_a^{*}=$ 0 correspond to the disease-free equilibrium point or the non-zero roots of

$\begin{array}{ll}p({\lambda }_a^{*})=D_2{({\lambda }_a^{*})}^2+D_1({\lambda }_a^{*})+D_0=0,& (94)\end{array}$

where

$\begin{array}{ll}\{\begin{array}{c}{{D_2=\mu }_{h}b}_m(T){\beta }_m{\pi }_c{\mu }_v(T)({\mu }_v(T)+{\alpha }_v){{\beta }_c\delta }_c{f}_2{f}_3{f}_7{f}_4{f}_5{f}_7\\ +{{{\mu }_h{b}_m(T){\beta }_m{\pi }_c{\mu }_v(T)({\mu }_v(T)+{\alpha }_v)\beta }_c\gamma \sigma }_c{{\delta }_a\delta }_c{f}_2{f}_3{f}_7\\ +{\gamma {\sigma }_a\delta }_a{{{\pi }_c({\mu }_v(T))}^2({\mu }_v(T)+{\alpha }_v)\beta }_c+\\ +{\gamma {\sigma }_a\delta }_a{{\pi }_c({\mu }_v(T))}^2({\mu }_v(T)+{\alpha }_v){\beta }_a{f}_1\\ -(\begin{array}{c}{{{\mu }_h{b}_m(T){\beta }_m{\pi }_c{\mu }_v(T)({\mu }_v(T)+{\alpha }_v)\beta }_c\delta }_c{\gamma {\sigma }_a\delta }_a\\ +f_2{f}_3{f}_7{f}_4{f}_5{f}_7{{{\pi }_c({\mu }_v(T))}^2({\mu }_v(T)+{\alpha }_v)\beta }_c\end{array}).\\ D_1={\delta }_a{\xi }_h{f}_2{f}_3{(f_7)}^2+{\mu }_h{b}_m(T){\beta }_m{\pi }_c{\mu }_v(T)({\mu }_v(T)+{\alpha }_v)\\ {\mu }_h{{{\beta }_c\delta }_{c}f}_2{f}_3{f}_7{f}_4{f}_5{f}_7\\ +{{\mu }_{h}b}_m(T){\beta }_a{\alpha }_v{{\mu }_h{b}_m(T){\beta }_m\eta }_I(T,R)M_A^{*}{{\beta }_c\delta }_c{f}_2{f}_3{f}_7{f}_4{f}_5{f}_7\\ +{{{{\mu }_{h}b}_m(T){\beta }_a{\alpha }_v{{\mu }_h{b}_m(T){\beta }_m\eta }_I(T,R)M_A^{*}\beta }_c\gamma \sigma }_c{{\delta }_a\delta }_c{f}_2{f}_3{f}_7-\\ (\begin{array}{c}{{{{\mu }_{h}b}_m(T){\beta }_a{\alpha }_v{{\mu }_h{b}_m(T){\beta }_m\eta }_I(T,R)M_A^{*}\beta }_c\delta }_c{\gamma {\sigma }_a\delta }_a\\ +{\mu }_h{f}_2{f}_3{f}_7{f}_4{f}_5{f}_7\\ {{{\pi }_c({\mu }_v(T))}^2({\mu }_v(T)+{\alpha }_v)\beta }_c{f}_2{f}_3{f}_7{f}_4{f}_5{f}_7\\ {{\pi }_c({\mu }_v(T))}^2({\mu }_v(T)+{\alpha }_v){\beta }_a{f}_1\end{array}),\\ D_0={{\mu }_{h}b}_m(T){\beta }_a{\alpha }_v{{\mu }_h{b}_m(T){\beta }_m\eta }_I(T,R)M_A^{*}{{\beta }_a\delta }_a{\xi }_h{f}_2{f}_3{(f_7)}^2\\ +f_2{f}_3{f}_4{f}_5{(f_7)}^2(\begin{array}{c}{{{\mu }_{h}b}_m(T){\beta }_a{\alpha }_v{{\mu }_h{b}_m(T){\beta }_m\eta }_I(T,R)M_A^{*}\mu }_h{{\beta }_c\delta }_c\\ -{\mu }_h{{\pi }_c({\mu }_v(T))}^2({\mu }_v(T)+{\alpha }_v){\beta }_a{f}_1\end{array}).\end{array}& (95)\end{array}$

Theorem 9. The endemic equilibrium $E^{*}=(S_c^{*},E_c^{*},I_c^{*},S_a^{*},E_a^{*},I_a^{*},W_a^{*},M_A^{*},S_v^{*},E_v^{*},I_v^{*})$ of the model system (13) is globally asymptotically stable in Ω if R₀ > 1.

Proof. To prove this, we define the Lyapunov function as follows:

$\begin{array}{c}L(S_c,E_c,I_c,S_a,E_a,I_a,W_a,M_A,S_v,E_v,I_v)\\ =\left(\begin{array}{c}\left({S_c-S}_c^{*}-S_c^{*}ln\frac{S_c}{S_c^{*}}\right)+\left({E_c-E}_c^{*}-E_c^{*}ln\frac{E_c}{E_c^{*}}\right)\\ +\left({I_c-I}_c^{*}-I_c^{*}ln\frac{I_c}{I_c^{*}}\right)+\\ \left({S_a-S}_a^{*}-S_a^{*}ln\frac{S_a}{S_a^{*}}\right)+\left({E_a-E}_a^{*}-E_a^{*}ln\frac{E_a}{E_a^{*}}\right)\\ +\left({I_a-I}_a^{*}-I_a^{*}ln\frac{I}{I_a^{*}}\right)\\ +\left({W_a-W}_a^{*}-I_a^{*}ln\frac{W_a}{W_a^{*}}\right)\end{array}\right)\\ +\left(\begin{array}{c}\left({M_A-M}_A^{*}-M_A^{*}ln\frac{M_A}{M_A^{*}}\right)+\left({S_v-S}_v^{*}-S_v^{*}ln\frac{S_v}{S_V^{*}}\right)\\ +\left({E_v-E}_v^{*}-E_v^{*}ln\frac{E_v}{E_V^{*}}\right)\\ +\left(I_v-I_v^{*}-I_v^{*}ln\frac{I_v}{I_V^{*}}\right)\end{array}\right).& (96)\end{array}$

L is continuously differentiable function and also

$\begin{array}{r}L(S_c^{*},E_c^{*},I_c^{*},S_a^{*},E_a^{*},I_a^{*},W_a^{*},M_A^{*},S_v^{*},E_v^{*},I_v^{*})=0,\text{ and}\\ L(S_c^{*},E_c^{*},I_c^{*},S_a^{*},E_a^{*},I_a^{*},W_a^{*},M_A^{*},S_v^{*},E_v^{*},I_v^{*})>0\\ forall(S_c^{*},E_c^{*},I_c^{*},S_a^{*},E_a^{*},I_a^{*},W_a^{*},M_A^{*},S_v^{*},E_v^{*},I_v^{*})\\ \ne (S_c,E_c,I_c,S_a,E_a,I_a,W_a,M_A,S_v,E_v,I_v),& (97)\end{array}$

and applying lemma[35], we obtain

$\begin{array}{c}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(L\right)\\ \le \left\{\begin{array}{c}{A}_1\left(1-\frac{S_c^{\ast }}{S_c}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{S}_c\left(t\right)+A_2\left(1-\frac{E_c^{\ast }}{E_c}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_c\left(t\right)\\ +A_3\left(1-\frac{I_c^{\ast }}{I_c}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_c\left(t\right)+\\ A_4\left(1-\frac{S_a^{\ast }}{S_a}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{S}_a\left(t\right)+A_5\left(1-\frac{E_a^{\ast }}{E_a}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_a\left(t\right)\\ +A_6\left(1-\frac{I_a^{\ast }}{I_a}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_a\left(t\right)\\ +A_7\left(1-\frac{W_a^{\ast }}{W_a}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{W}_a\left(t\right)+A_8\left(1-\frac{M_A^{\ast }}{M_A}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{M}_A\left(t\right)\\ A_9\left(1-\frac{S_v^{\ast }}{S_v}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{S}_v\left(t\right)\\ +A_{10}\left(1-\frac{E_v^{\ast }}{E_v}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{E}_v\left(t\right)+A_{11}\left(1-\frac{I_v^{\ast }}{I_v}\right)\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }{I}_I\left(t\right)\end{array}\right\}.& (98)\end{array}$

Applying Equation 14 from Equation 98, we obtain

$\begin{array}{cc}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(L\right)\le \left\{\begin{array}{c}\left(1-\frac{S_c^{\ast }}{S_c}\right)\left({\pi }_c-\left(\frac{b_m\left(T\right){\beta }_c{I}_v}{N_h}+{\xi }_h+{\mu }_h\right)\right)\\ +\left(1-\frac{E_c^{\ast }}{E_c}\right)\left(\frac{b_m\left(T\right){\beta }_c{I}_v}{N_h}{S}_c\left(t\right)-\left({\beta }_c+{\mu }_h\right)E_c\right)\\ +\left(1-\frac{I_c^{\ast }}{I_c}\right)\left({\beta }_c{E}_c\left(t\right)-\left({\mu }_h+d_1+{\sigma }_c\right)I_c\right)\\ +\left(1-\frac{S_a^{\ast }}{S_a}\right)\left({\xi }_h{S}_c+\gamma W_a-\left(\frac{b_m\left(T\right){\beta }_a{I}_v}{N_h}+{\mu }_h\right)S_a\right)\\ +\left(1-\frac{E_a^{\ast }}{E_a}\right)\left(\frac{b_m\left(T\right){\beta }_a{I}_v}{N_h}{S}_a-\left({\mu }_h+{\beta }_a\right)E_a\right)\\ +\left(1-\frac{I_a^{\ast }}{I_a}\right)\left({\beta }_a{E}_a-\left({\sigma }_a+{\mu }_h+d_2\right)I_a\right)\\ +\left(1-\frac{W_a^{\ast }}{W_a}\right)\left({\sigma }_c{I}_c+{\sigma }_a{I}_a-\left({\mu }_h+\gamma \right)W_{a,}\right)\\ +A_8\left(1-\frac{M_A^{\ast }}{M_A}\right)({\alpha }_I\left(T\right)\left(1-\frac{M_A}{K_v}\right)\\ \left(S_v+E_v+I_v\right)-\left({\eta }_I\left(T,R\right)+{\mu }_A\left(T\right)\right)M_A)\\ +\left(1-\frac{S_v^{\ast }}{S_v}\right)({\eta }_I\left(T,R\right)M_A-(\frac{b_m\left(T\right){\beta }_m\left(I_a+I_c\right)}{N_h}\\ +{\mu }_v\left(T\right))S_v)\\ +\left(1-\frac{E_v^{\ast }}{E_v}\right)\left(\frac{b_m\left(T\right){\beta }_m\left(I_a+I_c\right)}{N_h}{S}_v-\left({\mu }_v\left(T\right)+{\alpha }_v\right)E_v\right)\\ +\left(1-\frac{I_v^{\ast }}{I_v}\right)\left({\alpha }_v{E}_v-{\mu }_v\left(T\right)I_V\right)\end{array}\right\}.& (99)\end{array}$

By further rewriting this inequality, we have

$\begin{array}{cc}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(L\right)\le \left\{\begin{array}{c}\left(1-\frac{S_c^{\ast }}{S_c}\right)({\pi }_c-\left(\frac{b_m\left(T\right){\beta }_c{I}_v}{N_h}+{\xi }_h+{\mu }_h\right)\left(S_c-S_c^{\ast }\right)\\ -\left(\frac{b_m\left(T\right){\beta }_c{I}_v}{N_h}+{\xi }_h+{\mu }_h\right)\left(S_c^{\ast }\right))+\left(1-\frac{E_c^{\ast }}{E_c}\right)\\ (\frac{b_m\left(T\right){\beta }_c{I}_v}{N_h}{S}_c-\left({\beta }_c+{\mu }_h\right)\\ \left(E_c-E_c^{\ast }\right)-\left({\beta }_c+{\mu }_h\right)\left(E_c^{\ast }\right))\\ +\left(1-\frac{I_c^{\ast }}{I_c}\right)({\beta }_c{E}_c-\left({\mu }_h+d_1+{\sigma }_c\right)\left(I_c-I_c^{\ast }\right)\\ -\left({\mu }_h+d_1+{\sigma }_c\right)\left(I_c^{\ast }\right))\\ +\left(1-\frac{S_a^{\ast }}{S_a}\right)({\xi }_h{S}_c+\gamma W_a-\left(\frac{b_m\left(T\right){\beta }_a{I}_v}{N_h}+{\mu }_h\right)\\ \left(S_a-S_a^{\ast }\right)-\left(\frac{b_m\left(T\right){\beta }_a{I}_v}{N_h}+{\mu }_h\right)\left(S_a^{\ast }\right))\\ +\left(1-\frac{E_a^{\ast }}{E_a}\right)(\frac{b_m\left(T\right){\beta }_a{I}_v}{N_h}{S}_a-\left({\mu }_h+{\beta }_a\right)\\ \left(E_a-E_a^{\ast }\right)-\left({\mu }_h+{\beta }_a\right)\left(E_a^{\ast }\right))\\ +\left(1-\frac{I_a^{\ast }}{I_a}\right)({\beta }_a{E}_a-\left({\sigma }_a+{\mu }_h+d_2\right)\left(I_a-I_a^{\ast }\right)\\ -\left({\sigma }_a+{\mu }_h+d_2\right)\left(I_a^{\ast }\right))\\ +\left(1-\frac{W_a^{\ast }}{W_a}\right)({\sigma }_c{I}_c+{\sigma }_a{I}_a-\left({\mu }_h+\gamma \right)\\ \left(W_{a,}-W_a^{\ast }\right)-\left({\mu }_h+\gamma \right)\left(W_a^{\ast }\right))\\ +\left(1-\frac{M_A^{\ast }}{M_A}\right)\left(\begin{array}{c}{\alpha }_I\left(T\right)\left(1-\frac{M_A}{K_v}\right)\left(S_v+E_v+I_v\right)\\ -\left({\eta }_I\left(T,R\right)+{\mu }_A\left(T\right)\right)\left(M_A-M_A^{\ast }\right)\\ -\left({\eta }_I\left(T,R\right)+{\mu }_A\left(T\right)\right)\left(M_A^{\ast }\right)\end{array}\right)\\ +\left(1-\frac{S_v^{\ast }}{S_v}\right)\left(\begin{array}{c}{\eta }_I\left(T,R\right)M_A-(\frac{b_m\left(T\right){\beta }_m\left(I_a+I_c\right)}{N_h}\\ +{\mu }_v\left(T\right))\left(S_v-S_v^{\ast }\right)\\ -\left(\frac{b_m\left(T\right){\beta }_m\left(I_a+I_c\right)}{N_h}+{\mu }_v\left(T\right)\right)\left(S_v^{\ast }\right)\end{array}\right)\\ +\left(1-\frac{E_v^{\ast }}{E_v}\right)(\frac{b_m\left(T\right){\beta }_m\left(I_a+I_c\right)}{N_h}{S}_v\\ -\left({\mu }_v\left(T\right)+{\alpha }_v\right)\left(E_v-E_v^{\ast }\right)-\left({\mu }_v\left(T\right)+{\alpha }_v\right)\left(E_v^{\ast }\right))\\ +\left(1-\frac{I_v^{\ast }}{I_v}\right)({\alpha }_v{E}_v-{\mu }_v\left(T\right)\left(I_V-I_v^{\ast }\right)\\ -{\mu }_v\left(T\right)\left(I_v^{\ast }\right))\end{array}\right\}.& (100)\end{array}$

This inequality can be rewritten as

$\begin{array}{ll}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(L\right)=C_1-C_2, & (101)\end{array}$

where

$\begin{array}{c}{C}_1={\pi }_c(1-\frac{S_c^{*}}{S_c})+(\frac{b_m(T){\beta }_c{I}_v}{N_h}+{\xi }_h+{\mu }_h)\frac{{(S_c^{*})}^2}{S_c}\\ +\frac{b_m(T){\beta }_c{I}_v}{N_h}{S}_c+({\beta }_c+{\mu }_h)\frac{{(E_c^{*})}^2}{E_c}+{\beta }_c{E}_c+({\mu }_h+d_1+{\sigma }_c)\frac{{(I_c^{*})}^2}{I_c}\\ +{\xi }_h{S}_c+\gamma W_a+(\frac{b_m(T){\beta }_a{I}_v}{N_h}+{\mu }_h)\frac{{(S_a^{*})}^2}{S_a}\\ +\frac{b_m(T){\beta }_a{I}_v}{N_h}{S}_a+({\mu }_h+{\beta }_a)\frac{{(E_a^{*})}^2}{E_a}+{\beta }_a{E}_a+({\sigma }_a+{\mu }_h+d_2)\\ \frac{{(I_a^{*})}^2}{I_a}+{\sigma }_c{I}_c+{\sigma }_a{I}_a+({\mu }_h+\gamma )\frac{{(W_a^{*})}^2}{W_a}\\ +{\alpha }_I(T)(1-\frac{M_A}{K_v})(S_v+E_v+I_v)+({\eta }_I(T,R)+{\mu }_A(T))\frac{{(M_A^{*})}^2}{M_A}\\ +{\eta }_I(T,R)M_A+(\frac{b_m(T){\beta }_m(I_a+I_c)}{N_h}+{\mu }_v(T))\frac{{(S_v^{*})}^2}{S_v}\\ +\frac{b_m(T){\beta }_m(I_a+I_c)}{N_h}{S}_v+({\mu }_v(T)+{\alpha }_v)\\ \frac{{(E_v^{*})}^2}{E_v}+{\alpha }_v{E}_v+{\mu }_v(T)\frac{{(I_v^{*})}^2}{I_v}& (102)\end{array}$

and

$\begin{array}{c}{C}_2=({\frac{b_m(T){\beta }_c{I}_v}{N_h}+\xi }_h+{\mu }_h)\frac{S_c^{*}}{S_c}\\ +({\frac{b_m(T){\beta }_c{I}_v}{N_h}+\xi }_h+{\mu }_h)\frac{{(S_c-S_c^{*})}^2}{S_c}+({\frac{b_m(T){\beta }_c{I}_v}{N_h}+\xi }_h+{\mu }_h)\\ S_c^{*}+\frac{b_m(T){\beta }_c{I}_v}{N_h}{S}_c\frac{E_c^{*}}{E_c}+({\beta }_c+{\mu }_h)\frac{{(E_c-E_c^{*})}^2}{E_c}+({\beta }_c+{\mu }_h)E_c^{*}\\ +{\beta }_c{E}_c\frac{I_c^{*}}{I_c}+({\mu }_h+d_1+{\sigma }_c)\frac{{(I_c-I_c^{*})}^2}{I_c}\\ +({\mu }_h+d_1+{\sigma }_c)I_c^{*}+({\xi }_h{S}_c+\gamma W_a)\frac{S_a^{*}}{S_a}\\ +(\frac{b_m(T){\beta }_a{I}_v}{N_h}+{\mu }_h)\frac{{(S_a-S_a^{*})}^2}{S_a}\\ +(\frac{b_m(T){\beta }_a{I}_v}{N_h}+{\mu }_h)S_a^{*}+\frac{b_m(T){\beta }_a{I}_v}{N_h}{S}_a\frac{E_a^{*}}{E_a}+({\mu }_h+{\beta }_a)\\ \frac{{(E_a-E_a^{*})}^2}{E_a}+({\mu }_h+{\beta }_a)(E_a^{*})\\ +{\beta }_a{E}_a\frac{I_a^{*}}{I_a}+({\sigma }_a+{\mu }_h+d_2)\frac{{(I_a-I_a^{*})}^2}{I_a}+({\sigma }_a+{\mu }_h+d_2)I_a^{*}\\ +({\sigma }_c{I}_c+{\sigma }_a{I}_a)\frac{W_a^{*}}{W_a}+({\mu }_h+\gamma )\frac{{(W_{a,}-W_a^{*})}^2}{W_{a,}}\\ +({\mu }_h+\gamma )W_a^{*}+{\alpha }_I(T)(1-\frac{M_A}{K_v})(S_v+E_v+I_v)\frac{M_A^{*}}{M_A}\\ +({\eta }_I(T,R)+{\mu }_A(T))\frac{{(M_A-M_A^{*})}^2}{M_A}+({\eta }_I(T,R)+{\mu }_A(T))M_A^{*}\\ +{\eta }_I(T,R)M_A\frac{S_v^{*}}{S_v}+(\frac{b_m(T){\beta }_m(I_a+I_c)}{N_h}+{\mu }_v(T))\frac{{(S_v-S_v^{*})}^2}{S_v}\\ +(\frac{b_m(T){\beta }_m(I_a+I_c)}{N_h}+{\mu }_v(T))\\ +\frac{b_m(T){\beta }_m(I_a+I_c)}{N_h}{S}_v\frac{E_v^{*}}{E_v}+({\mu }_v(T)+{\alpha }_v)\\ \frac{{(E_v-E_v^{*})}^2}{E_v}+({\mu }_v(T)+{\alpha }_v)E_v^{*}+{\alpha }_v{E}_v\frac{I_v^{*}}{I_v}\\ +{\mu }_v(T)\frac{{(I_V-I_v^{*})}^2}{I_V}+{\mu }_v(T)I_v^{*}.& (103)\end{array}$

As a result, $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }(L)<0$ for C₁ < C₂ and $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }(L)=0$ if and only if $S_c^{*}=S_c,E_c^{*}=E_c,I_c^{*}=I_c,S_a^{*}=S_a,E_a^{*}=E_a,I_a^{*}=I_a,W_a^{*}=W_a,M_A^{*}=M_A,S_v^{*}=S_v,E_v^{*}=E_v,$ and $I_v^{*}=I_v$. The largest closed and bounded invariant set in {$(S_c,E_c,I_c,S_a,E_a,I_a,W_a,M_A,S_v,E_v,I_v)\in R_{+}^{11}$: $\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }(L)=0$ is the singleton {E^*}, where E^* is the endemic equilibrium point. As a result, when R₀ > 1 in the region Ω, the unique equilibrium point E^* is globally asymptotically stable, according to the LaSalle invariance principle [38]. This completes the proof of the theorem.

5 Sensitivity analysis

To ensure model predictions remain reliable despite potential uncertainties in parameter values, sensitivity analysis is widely conducted. It reveals how changes in parameters affect the system's overall dynamics, particularly in relation to the basic reproduction number, offering crucial insights. Sensitivity analysis using the fixed-point estimation method described in reference [39] is applied to the fractional malaria model (14). This method analyzes the effects of local changes in model parameters by calculating the normalized forward sensitivity index of a variable v to a parameter (p), which is defined as

$\begin{array}{ll}{\pi }_p^v=\frac{\partial v}{\partial p}\times \frac{p}{v}.& (104)\end{array}$

We created Table 3 using the following values: b_m = 0.29, δ_c = 0.0003, β_m = 0.022, μ_h = 0.00005, π_c = 540, η = 0.343, d₁ = 0.0002, μ_A = 0.1041, α_I = 1.84, ξ_h = 0.000161, and μ_v = 0.05 per day, K_v = 40000, α_v = 0.5, and σ_c = 0.002 as given in [[8] and the references therein]. Figure 1 illustrates a schematic of the mathematical model for malaria transmission (Equation 12), based on the works of [8, 10].

Table 3

Table 3. Temperature-dependent parameters [see [8, 14, 41, 42] and there references therein].

Figure 1

Figure 1. A schematic of the mathematical model for malaria transmission.

Figure 2 depicts the sensitivity analysis of the basic reproduction number [fractional malaria model (14)] for the 14 parameters in Table 4.

Figure 2

Figure 2. Sensitivity to model parameters (Table 2).

Table 4

Table 4. Sensitivity analysis of the basic reproduction number (R₀) at parameter values given above.

From the analysis of Table 4 and Figure 2, it is evident that each parameter has a positive or negative effect on the basic reproduction number (R₀). Parameters with positive signs, such as b_m, δ_c, β_m, α_v, μ_h, η, k_v, and α_I, increase R₀, while those with negative signs, such as π_c, d₁, μ_A, ξ_h, σ_c, and μ_v, decrease it. The parameter with a higher sensitivity index magnitude is more influential than those with smaller magnitudes, as exhibited in Figure 2. For instance, among the parameters given in the fractional malaria model (14) relative to R₀, the per capita death rate for adult mosquitoes (μ_v), the lifetime number of eggs laid (b_m), the maturation rate of immature mosquitoes (η), vector carrying capacity(k_v), and the rate at which exposed children transition to the infected class (δ_c) are the most sensitive parameters, in that order. Therefore, to eliminate or control malaria disease, it is important to focus on controlling these parameters.

6 Model analysis with climate-dependent (temperature and rainfall) parameter

This section examines the model parameters that affect malaria transmission dynamics focusing on temperature and rainfall. The results presented below were obtained using MATLAB software.

6.1 The mosquito maturation rate

The mosquito maturation rate, denoted by η(T,R), depends on both temperature (T) and rainfall (R). It determines the rate at which immature mosquitoes develop into mature adults, as described by the following equation [see [8, 14, 41] and the references therein].

$\begin{array}{ll}\eta (T,R)=\frac{B(T)P_E(R)P_L(R)P_P(R)P_L(T)}{TEA(T)},& (105)\end{array}$

where

• $B(T)=\frac{-0\mathrm{\text{:}}153T^2+8\mathrm{\text{:}}61T-97.7}{-ln(-0.000828T^2+0.0367T+0.522)}$, where B(T) is the lifetime number of eggs laid,

• $P_E(R)=\frac{4\times 0.93}{2500}R(50-R)$, where P_E(R) is the daily survival probabilities of eggs,

• $P_L(R)=\frac{4\times 0.25}{2500}R(50-R)$, where P_L(R) is the daily survival probabilities of larva,

• $P_P(R)=\frac{4\times 0.75}{2500}R(50-R)$, whereP_P(R) is the daily survival probabilities of pupae,

• $P_L(T)=e^{-(0\mathrm{\text{:}}00554T-0\mathrm{\text{:}}06737)}$, where P_L(T) is the temperature-dependent daily probability of survival of larvae,

• $TEA(T)=\frac{1}{-0\mathrm{\text{:}}00094T^2+0.049T-0.552}$, where TEA(T) is the development time from egg to adult mosquito.

• Thus, $\eta (T,R)=\frac{-0\mathrm{\text{:}}153T^2+8\mathrm{\text{:}}61T-97.7}{-ln(-0.000828T^2+0.0367T+0.522)}\frac{4^{*}.93}{2500}R$ $(50-R)\frac{4^{*}.25}{2500}R(50-R)\frac{4^{*}.75}{2500}R(50-R)$ $e^{-(0\mathrm{\text{:}}00554T-0\mathrm{\text{:}}06737)\frac{1}{-0\mathrm{\text{:}}00094T^2+0.049T-0.552}}$

Based on the study of malaria transmission shown in [8, 14, 40, 41] and Figure 3, we examined the effects of temperature on the infected mosquito population in our proposed integer-order malaria model (14). We simulated infected mosquito populations across four temperature ranges: 17–25°C, 21–25°C, 30–32°C, and 35–39°C. The results are shown in Figures 4A–D, respectively.

Figure 3

Figure 3. Relationship between basic reproduction number vs. temperature.

Figure 4

Figure 4. Model simulations 4 show new mosquito infections for temperatures: (A) (17–25°C), (B) (21–25°C), (C) (30–32°C), and (D) (35–39°C).

Furthermore, we investigated the effects of daily rainfall on mosquito development. Mosquito burden are known to peak at 25°C [[8, 40–42] and Figure 4A], so we used this constant temperature (T = 25°C) for our analysis, as in previous studies [see [42] and the reference therein]. Figure 5 depicts the relationship between mosquito maturation rate (η(T,R)) and rainfall (R) in millimeters for a temperature of 25°C. As shown in Figure 5, aquatic mosquitoes cannot survive daily rainfall exceeding 50 millimeters. It is important to note that aquatic mosquitoes cannot survive daily rainfall exceeding 50 mm, which limits vector population growth. These effects are illustrated in Figures 6A–D.

Figure 5

Figure 5. Mosquito maturation rate (η (T,R)) vs. rainfall (R) in mm for (T = 25).

Figure 6

Figure 6. (A–D) Show mosquito stages (aquatic, susceptible, exposed, and infected) at 25°C over time (days) for different rainfall rates (10, 20, 30, and 40 mm).

The mosquito maturation rate, denoted by η(T,R), depends on both temperature (T) and rainfall (R). This rate determines the speed at which immature mosquitoes develop into mature adults. We examined this rate for various temperature and rainfall values: T = 20, 25, 30, and 35°C and R = 10, 20, 30, and 40 mm. The results are shown in Figure 6.

7 Results and discussion

7.1 Results

This section presents fractional-order (14) malaria models using graphs to understand the behavior of the solution trajectories. To draw the graph of the fractional-order model (13), we use the scheme introduced in [43], that is, from the solution of the differential equation with fractional order α given by

$\begin{array}{ll}\{\begin{array}{c}\begin{array}{c}ABC\\ 0\end{array}{D}_t^{\alpha }\left(g\left(t\right)\right)=F\left(t,g\left(t\right)\right),\\ g\left(0\right)=g_0\ge 0,\end{array}wherefor\alpha \in (0,1],& (106)\end{array}$

expressed as

$\begin{array}{r}g(t)=g(0)+\frac{1-\alpha }{B(\alpha )}F(t,g(t))\\ +\frac{\alpha }{B(\alpha )\Gamma (\alpha )}{\displaystyle {\int }_0^t}{F(\tau ,g(\tau ))(t-\tau )}^{\alpha -1}d\tau ,& (107)\end{array}$

where

$\begin{array}{ll}g(t)={\left(\begin{array}{c}{S}_c(t),S_a(t),E_c(t),E_a(t),I_c(t),I_a(t),\\ W_a(t),M_A(t),S_v(t),E_v(t),I_v(t)\end{array}\right)}^T.& (108)\end{array}$

At t = t_n+1, n = 0, 1, 2, ⋯  , we obtain

$\begin{array}{c}g(t_{n+1})-g(0)=\frac{1-\alpha }{B(\alpha )}F(t_n,g(t_n))+\frac{\alpha }{B(\alpha )\Gamma (\alpha )}\\ {\displaystyle {\int }_0^{t_{n+1}}}{F(\tau ,g(\tau ))(t_{n+1}-\tau )}^{n-1}d\tau \\ =\frac{1-\alpha }{B(\alpha )}F(t_n,g(t_n))+\frac{\alpha }{B(\alpha )\Gamma (\alpha )}\\ \sum _{i=0}^n{\displaystyle {\int }_{t_i}^{t_{i+1}}}{F(\tau ,g(\tau ))(t_{n+1}-\tau )}^{n-1}d\tau .& (109)\end{array}$

With the help of interpolation polynomial, we approximate the function f(τ, g(τ))over [t_i, t_i+1].

$\begin{array}{c}F(\tau ,g(\tau ))\cong p_k(\tau )=\frac{F(t_i,g(t_i))}{h}(\tau -t_{i-1})\\ -\frac{F(t_{i-1},g(t_{i-1}))}{h}(\tau -t_i).& (110)\end{array}$

Using Equation 110, Equation 109 takes the form:

$\begin{array}{c}g\left(t_{n+1}\right)=g\left(0\right)+\frac{1-\alpha }{B\left(\alpha \right)}F\left(t_n,g\left(t_n\right)\right)\\ +\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ ({\displaystyle \sum _{i=0}^n\frac{f\left(t_i,g\left(t_i\right)\right)}{h}{\displaystyle {\int }_{t_i}^{t_{i+1}}\left(\tau -t_{i-1}\right){\left(t_{n+1}-\tau \right)}^{n-1}}d}\tau \\ -\frac{f\left(t_{i-1},g\left(t_{i-1}\right)\right)}{h}{\displaystyle {\int }_{t_i}^{t_{i+1}}\left(\tau -t_i\right){\left(t_{n+1}-\tau \right)}^{n-1}})d\tau .& (111)\end{array}$

Solving the integrals involved in Equation 111, we get the approximate solution as below:

$\begin{array}{c}g\left(t_{n+1}\right)=g\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}F\left(t_n,g\left(t_n\right)\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{l}{h}^{\alpha }F\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }F\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}& (112)\end{array}$

Hence, we have the following recursive formulas for the proposed malaria model (14):

$\begin{array}{c}{S}_c\left(t_{n+1}\right)=S_c\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_1\left(t_n,g\left(t_n\right)\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_1\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_1\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}, & (113)\end{array}$

$\begin{array}{c}{E}_c\left(t_{n+1}\right)=E_c\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_2\left(t_n,g\left(t_n\right)\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_2\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_2\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}, & (114)\end{array}$

$\begin{array}{c}{I}_c\left(t_{n+1}\right)=I_c\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_3\left(t_n,g\left(t_n\right)\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_3\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_3\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}, & (115)\end{array}$

$\begin{array}{c}{S}_a\left(t_{n+1}\right)=S_a\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_4\left(t_n,g\left(t_n\right)\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_4\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_4\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}, & (116)\end{array}$

$\begin{array}{c}{E}_a\left(t_{n+1}\right)=E_a\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_5\left(t_n,g\left(t_n\right)\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_5\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_5\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}, & (117)\end{array}$

$\begin{array}{c}{I}_a\left(t_{n+1}\right)=I_a\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_6\left(t_n,g\left(t_n\right)\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_6\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_6\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}, & (118)\end{array}$

$\begin{array}{c}{W}_a\left(t_{n+1}\right)=W_a\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_7\left(t_n,g\left(t_n\right)\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_7\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_7\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)},& (119)\end{array}$

$\begin{array}{c}{M}_A\left(t_{n+1}\right)=M_A\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_8\left(t_n,g\left(t_n\right)\right)+\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_8\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_8\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}, & (120)\end{array}$

$\begin{array}{c}{S}_v\left(t_{n+1}\right)=S_v\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_9\left(t_n,g\left(t_n\right)\right)\\ +\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_9\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\\ \left(n-i+2+\alpha \right)-{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_9\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}, & (121)\end{array}$

$\begin{array}{c}{E}_v\left(t_{n+1}\right)=E_v\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_{10}\left(t_n,g\left(t_n\right)\right)\\ +\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_{10}\left(t_i,g\left(t_i\right)\right)({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)\\ -{\left(n-i\right)}^{\alpha }\left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_{10}\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ \left({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right)\right)\end{array}\end{array}\right)}, & (122)\end{array}$

$\begin{array}{c}{I}_v\left(t_{n+1}\right)=I_v\left(t_0\right)+\frac{1-\alpha }{B\left(\alpha \right)}{F}_{11}\left(t_n,g\left(t_n\right)\right)\\ +\frac{\alpha }{B\left(\alpha \right)\Gamma \left(\alpha \right)}\\ {\displaystyle \sum _{i=0}^n\left(\begin{array}{c}\begin{array}{c}{h}^{\alpha }{F}_{11}\left(t_i,g\left(t_i\right)\right)\\ ({\left(n+1-i\right)}^{\alpha }\left(n-i+2+\alpha \right)-{\left(n-i\right)}^{\alpha }\\ \left(n-i+2+2\alpha \right))\\ -h^{\alpha }{F}_{11}\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ ({\left(n+1-i\right)}^{\alpha +1}-{\left(n-i\right)}^{\alpha }\left(n-i+1+\alpha \right))\end{array}\end{array}\right)}, & (123)\end{array}$

where

$\begin{array}{c}{F}_j\left(t_n,g\left(t_n\right)\right)\\ =F_j\left(t_n,\left(\begin{array}{c}{S}_c\left(t_n\right),S_a\left(t_n\right),E_c\left(t_n\right),E_a\left(t_n\right),I_c\left(t_n\right),I_a\left(t_n\right),\\ W_a\left(t_n\right),M_A\left(t_n\right),S_v\left(t_n\right),E_v\left(t_n\right),I_v\left(t_n\right)\end{array}\right)\right),\\ j=1,2,\cdots ,11& (124)\end{array}$

$\begin{array}{c}{F}_i\left(t_i,g\left(t_i\right)\right)\\ =F_i\left(t_i,\left(\begin{array}{c}{S}_c\left(t_i\right),S_a\left(t_i\right),E_c\left(t_i\right),E_a\left(t_i\right),I_c\left(t_i\right),I_a\left(t_i\right),\\ W_a\left(t_i\right),M_A\left(t_i\right),S_v\left(t_i\right),E_v\left(t_i\right),I_v\left(t_i\right)\end{array}\right)\right),\\ i=1,2,\cdots ,11& (125)\end{array}$

$\begin{array}{c}{F}_i\left(t_{i-1},g\left(t_{i-1}\right)\right)\\ =F_i\left(t_{i-1},\left(\begin{array}{c}{S}_c\left(t_{i-1}\right),S_a\left(t_{i-1}\right),E_c\left(t_{i-1}\right),E_a\left(t_{i-1}\right),\\ I_c\left(t_{i-1}\right),I_a\left(t_{i-1}\right),W_a\left(t_{i-1}\right),M_A\left(t_{i-1}\right),\\ S_v\left(t_{i-1}\right),E_v\left(t_{i-1}\right),I_v\left(t_{i-1}\right)\end{array}\right)\right),\\ i=1,2,\cdots ,11.& (126)\end{array}$

7.2 Discussion

Fractional-order malaria model analysis considering temperature and rainfall with Caputo operators showed correlation between these factors and mosquito population dynamics. Key factors affecting mosquito dynamics identified through sensitivity analysis were adult mosquito death rate, egg laying rate, maturation rate, vector carrying capacity, and exposed child transition rate, aligning with prior studies [8].

Based on the study of malaria transmission shown in [8, 14, 40–42] and Figure 3, we examined the effects of temperature on the infected mosquito population in our proposed integer-order malaria model 14(α = 1). We varied temperatures across four ranges: 17–25°C, 21–25°C, 30–32°C, and 35–39°C. The results are shown in Figures 4A–D, respectively. Figure 4A confirms that within the 17–25°C range, both infections and malaria burden rise with temperature, with the optimal temperature being 25°C [as expected from [8, 14, 40–42]]. However, Figures 4B–D show varying trends at higher temperatures: a peak burden at higher temperatures [Figure 4B, supported by [8]], a decrease with increasing temperature [Figure 4C, consistent with [8] and the references therein and [40]], and a drastic decrease with increasing temperature (Figure 4D) [see [15] and the references therein and [14, 40]].

As Figure 6A shows, the maximum value for aquatic mosquito vector growth is observed at a rainfall of 40 mm. In contrast, the remaining mosquito stages, including the susceptible mosquito (Figure 6B), exposed mosquito (Figure 6C), and infected mosquito (Figure 6D), all require a maximum rainfall of 30 mm for growth. Notably, infected mosquito populations peak at 30 mm of rainfall, suggesting that this is a favorable condition for malaria transmission compared to other rainfall values considered. Thus, we conclude that the peak for malaria transmission is at temperature (T = 25°C) [8, 40–42] and rainfall (R = 30 mm) [42]. Figure 7 reveals that mosquito maturation rate peaks and reaches its minimum at 20°C (19.18 and 10 mm, respectively) despite variations in rainfall (10–40 mm).

Figure 7

Figure 7. Impact of temperature and daily rainfall on mosquito maturation rate. Values used: T = 20, 25, 30, and 35°C and R = 10, 20, 30, and 40 mm.

Figure 8 shows mosquito infection peaks at 25°C for all fractional orders (α = 0.7, 0.8, 0.9, and 1), similar to the classical malaria model 14 (Figure 4A) and aligned with prior findings [8, 40–42]. In addition, Figures 8A–D (α = 0.7, 0.8, 0.9, and 1) reveal a rise in mosquito infections with increasing temperature, consistent with the classical model 14 (Figure 4A) [8], whereas Figures 9A–D for all fractional orders (α = 0.7, 0.8, 0.9, and 1) show a dramatic reduction in mosquito infection as the temperature increases from 35 to 39°C. This is similar to the graph of the classical malaria model (14) as in Figure 4D, which is consistent with the findings in [15] and the reference therein and [40].

Figure 8

Figure 8. (A–D) Exemplify the total numbers of infected vectors for 17–25°C for different values of fractional order alpha: (A) (α = 0.7), (B) (α = 0.8), (C) (α = 0.9), and (D) (α = 1), respectively.

Figure 9

Figure 9. (A–D) Depict the total numbers of infected vectors for 35–39°C for different values of fractional order alpha: (A) (α = 0.7), (B) (α = 0.8), (C) (α = 0.9), and (D) (α = 1), respectively.

Furthermore, from Figures 8, 9, we can see that as the value of the fractional order alpha approaches 1, the results resemble the graph of the classical malaria model 14. For instance, Figures 8A–D (α = 0.7, 0.8, 0.9, and 1) resemble the classical malaria model 14 (α = 1) in Figure 4A, and Figures 9A–D (α = 0.7, 0.8, 0.9, and 1) resemble the classical malaria model 14 (α = 1) (Figure 4D).

Furthermore, Figure 10 shows that, like Figures 8A–D plotted on the different plane, the burden of mosquito peaks shown to occur at 25°C. As the value of the fractional-order derivative approaches one, it resembles Figure 4A, the classical model of 14 (α = 1). One advantage of ABC fractional derivative operators is that we can obtain distinct solutions.

Figure 10

Figure 10. Infected vector vs. time (days) for different temperatures 17–25°C and different alphas (α = 0.6, 0.7, 0.8, 0.9, and 1).

Figure 11A demonstrates that as the number of susceptible children under seventeen decreases, the value of alpha (α) also decreases. This suggests a linear relationship between the two. Conversely, Figure 11B reveals that decreasing the exposed class of children under seventeen lowers alpha (α). This implies faster transitions from susceptible to exposed with a lower alpha value. Finally, Figure 11C shows a steeper curve. This indicates that model (14) relies heavily on past infection data when determining the rate of change in infected individuals.

Figure 11

Figure 11. (A–C) Show the total numbers of susceptible, exposed, and infected children under 17, for different values of alpha.

Figures 12A–C examine adult susceptibility and disease spread. Figure 12A shows a rise and fall in susceptible adults, reflecting growth (new adults) followed by depletion (infections/immunity). Figure 12B reveals slower transitions from susceptible to exposed adults with lower alpha (α). Conversely, Figure 12C shows infected mosquitoes raise with higher alpha (α), implying the model prioritizes past infections in adult disease spread. The decline in infected adults with higher alpha (α) suggests the model incorporates other factors later.

Figure 12

Figure 12. (A–C) Show the total numbers of susceptible, exposed, and infected adults, for different values of alpha.

Furthermore, for our future study, research directions [44–49] will be used. These areas, particularly those related to non-standard finite difference methods in fractional modeling [44, 49] and delay techniques in epidemic models [45–48], align well with the potential applications of our proposed method.

8 Conclusion

This study analyzes malaria transmission dependence on temperature and rainfall using a fractional-order differential model with Atangana–Baleanu operators (in the Caputo sense). It confirms the model's solution existence, uniqueness, and stability. Sensitivity analysis identified adult mosquito mortality, egg laying, maturation, and exposed child transition rates as key factors, aligning with prior research [8].

Simulations confirmed theoretical results: Peak malaria transmission occurred at 25°C (Figure 4A), consistent with [8, 40–42]. Malaria burden increased with temperature (17–25°C) (Figure 4A) [agreeing with [8]]. Figures 4B–D show a decrease with increasing temperature, also supported by [13–15], and the references therein.

Fractional-order model simulations (alpha approaching 1) resembled the classical model. For instance, Figures 8A–D (α = 0.7, 0.8, 0.9, and 1) resemble the classical malaria model in Figure 4A, and Figures 9A–D (α = 0.7, 0.8, 0.9, and 1) resemble Figure 4D. This study lays the groundwork for future research on infectious diseases using fractional derivatives, particularly ABC operators. Further extension could incorporate real-data non-autonomous parts, requiring additional compartments for mosquito and human populations.

Data availability statement

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

Ethics statement

Written informed consent was obtained from the individual(s), and minor(s)' legal guardian/next of kin, for the publication of any potentially identifiable images or data included in this article.

Author contributions

AG: Writing – original draft, Writing – review & editing, Conceptualization, Data curation, Formal analysis, Investigation, Methodology. CD: Writing – review & editing, Conceptualization, Formal analysis, Investigation, Software, Supervision.

Funding

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

Conflict of interest

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

Publisher's note

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

References

1. Gbadamosi B, Adebimpe O, Ojo MM, Oludoun O, Abiodun O, Adesina I. Modeling the impact of optimal control measures on the dynamics of cholera. Model Earth Syst Environ. (2022) 9:1387–400. doi: 10.1007/s40808-022-01570-9

Crossref Full Text | Google Scholar

2. Malaria (2023). Available at: https://www.who.int/news-room/questions-and-answers/item/malaria (accessed November 30, 2023).

Google Scholar

3. Xing Y, Guo Z, Liu J. Backward bifurcation in a malaria transmission model. J Biol Dyn. (2020) 14:368–88. doi: 10.1080/17513758.2020.1771443

Crossref Full Text | Google Scholar

4. Lakew YY, Fikrie A, Godana SB, Wariyo F, Seyoum W. Magnitude of malaria and associated factors among febrile adults in Siraro District Public Health facilities, West Arsi Zone, Oromia, Ethiopia 2022: a facility-based cross-sectional study. Malar J. (2023) 22:259. doi: 10.1186/s12936-023-04697-x

PubMed Abstract | Crossref Full Text | Google Scholar

5. Ndamuzi E, Gahungu P. Mathematical modeling of malaria transmission dynamics: case of Burundi. J Appl Math Phys. (2021) 09:2447–60. doi: 10.4236/jamp.2021.910156

PubMed Abstract | Crossref Full Text | Google Scholar

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

Google Scholar

7. Macdonald G. Epidemiology and Control of Malaria. London: Oxford University Press (1957).

Google Scholar

8. Okuneye K, Gumel AB. Analysis of a temperature- and rainfall-dependent model for malaria transmission dynamics. Math Biosci. (2017) 287:72–92. doi: 10.1016/j.mbs.2016.03.013

PubMed Abstract | Crossref Full Text | Google Scholar

9. Traoré B, Sangaré B, Traoré S. A mathematical model of malaria transmission with structured vector population and seasonality. J Appl Math. (2017) 2017:1–15. doi: 10.1155/2017/6754097

Crossref Full Text | Google Scholar

10. Ducrot A, Sirima SB, Somé B, Zongo P. A mathematical model for malaria involving differential susceptibility, exposedness and infectivity of human host. J Biol Dyn. (2009) 3:574–98. doi: 10.1080/17513750902829393

Crossref Full Text | Google Scholar

11. Oheneba-Dornyo TV, Amuzu S, Maccagnan A, Taylor T. Estimating the impact of temperature and rainfall on malaria incidence in Ghana from 2012 to 2017. Environ Model Assess. (2022) 27:473–89. doi: 10.1007/s10666-022-09817-6

Crossref Full Text | Google Scholar

12. Parham PE, Michael E. Modeling the effects of weather and climate change on malaria transmission. Environ Health Perspect. (2010) 118:620–6. doi: 10.1289/ehp.0901256

Crossref Full Text | Google Scholar

13. Craig MH, Snow RW, le Sueur D. A climate-based distribution model of malaria transmission in sub-saharan Africa. Parasitol Today. (1999) 15:105–11. doi: 10.1016/s0169-4758(99)01396-4

Crossref Full Text | Google Scholar

14. Agusto FB, Gumel AB, Parham PE. Qualitative assessment of the role of temperature variations on malaria transmission dynamics. J Biol Syst. (2015) 23:1550030. doi: 10.1142/s0218339015500308

Crossref Full Text | Google Scholar

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

Crossref Full Text | Google Scholar

16. Shah K, Sinan M, Abdeljawad T, El-Shorbagy MA, Abdalla B, Abualrub MS, et al. Detailed study of a fractal-fractional transmission dynamical model of viral infectious disease with vaccination. Complexity. (2022) 2022:1–21. doi: 10.1155/2022/7236824

Crossref Full Text | Google Scholar

17. Sidi Ammi MR, Tahiri M, Tilioua M, Zeb A, Khan I, Andualem M. Global analysis of a time fractional order spatio-temporal SIR model. Sci Rep. (2022) 12:5751. doi: 10.1038/s41598-022-08992-6

Crossref Full Text | Google Scholar

18. Mandal M, Jana S, Nandi SK, Kar TK. Modelling and control of a fractional-order epidemic model with fear effect. Energy Ecol Environ. (2020) 5:421–32. doi: 10.1007/s40974-020-00192-0

PubMed Abstract | Crossref Full Text | Google Scholar

19. Windarto, Khan MA, Fatmawati. Parameter estimation and fractional derivatives of dengue transmission model. AIMS Math. (2020) 5:2758–79. doi: 10.3934/math.2020178

Crossref Full Text | Google Scholar

20. Tuan NH, Mohammadi H, Rezapour S. A mathematical model for COVID-19 transmission by using the Caputo fractional derivative. Chaos Solit Fract. (2020) 140:110107. doi: 10.1016/j.chaos.2020.110107

Crossref Full Text | Google Scholar

21. Atangana A, Baleanu D. New fractional derivatives with nonlocal and non-singular kernel: theory and application to heat transfer model. Therm Sci. (2016) 20:763–9. doi: 10.2298/tsci160111018a

Crossref Full Text | Google Scholar

22. Caputo M, Fabrizio M. A new definition of fractional derivative without singular kernel. Prog Fract Differ Appl. (2015) 9:15–23.

Google Scholar

23. Caputo M. Linear models of dissipation whose Q is almost frequency independent–II. Geophys J Int. (1967) 13:529–39. doi: 10.1111/j.1365-246x.1967.tb02303.x

Crossref Full Text | Google Scholar

24. Gizaw AK, Deressa CT. Analysis of age-structured mathematical model of malaria transmission dynamics via classical and ABC fractional operators. Math Probl Eng. (2024) 2024:1–24. doi: 10.1155/2024/3855146

Crossref Full Text | Google Scholar

25. Ullah S, Altaf Khan M, Farooq MO, Alzahrani E. A fractional model for the dynamics of tuberculosis (TB) using Atangana-Baleanu derivative. Discrete Cont Dyn Syst. (2020) 13:937–56. doi: 10.3934/dcdss.2020055

Crossref Full Text | Google Scholar

26. Yavuz M. European option pricing models described by fractional operators with classical and generalized Mittag-Leffler kernels. Numer Methods Partial Differ Equ. (2021) 37:342–62. doi: 10.1002/num.22645

Crossref Full Text | Google Scholar

27. Deressa CT, Duressa GF. Analysis of Atangana–Baleanu fractional-order SEAIR epidemic model with optimal control. Adv Differ Equ. (2021) 2021:174. doi: 10.1186/s13662-021-03334-8

PubMed Abstract | Crossref Full Text | Google Scholar

28. Deressa CT, Etemad S, Rezapour S. On a new four-dimensional model of memristor-based chaotic circuit in the context of nonsingular Atangana–Baleanu–Caputo operators. Adv Differ Eq. (2021) 2021:444. doi: 10.1186/s13662-021-03600-9

Crossref Full Text | Google Scholar

29. Yadeta DM, Gizaw AK, Mussa YO. Approximate analytical solution of one-dimensional beam equations by using time-fractional reduced differential transform method. J Appl Math. (2020) 2020:1–13. doi: 10.1155/2020/7627385

Crossref Full Text | Google Scholar

30. Odibat ZM, Shawagfeh NT. Generalized Taylor's formula. Appl Math Comput. (2007) 186:286–93. doi: 10.1016/j.amc.2006.07.102

Crossref Full Text | Google Scholar

31. Adel W, Elsonbaty A, Aldurayhim A, El-Mesady A. Investigating the dynamics of a novel fractional-order monkeypox epidemic model with optimal control. Alexand Eng J. (2023) 73:519–42. doi: 10.1016/j.aej.2023.04.051

Crossref Full Text | Google Scholar

32. Diekmann O, Heesterbeek JAP, Metz JAJ. On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations. J Math Biol. (1990) 28:365–82. doi: 10.1007/bf00178324

PubMed Abstract | Crossref Full Text | Google Scholar

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

Crossref Full Text | Google Scholar

34. Ahmed E, El-Sayed AMA, El-Saka HAA. On some Routh–Hurwitz conditions for fractional order differential equations and their applications in Lorenz, Rössler, Chua and Chen systems. Phys Lett A. (2006) 358:1–4. doi: 10.1016/j.physleta.2006.04.087

Crossref Full Text | Google Scholar

35. Vargas-De-León C. Volterra-type Lyapunov functions for fractional-order epidemic systems. Commun Nonlinear Sci Numer Simul. (2015) 24:75–85. doi: 10.1016/j.cnsns.2014.12.013

PubMed Abstract | Crossref Full Text | Google Scholar

36. Bebernes JW. The stability of dynamical systems (J. P Lasalle). SIAM Rev. (1979) 21:418–20. doi: 10.1137/1021079

Crossref Full Text | Google Scholar

37. Sahnoune MY, Ez-Zetouni A, Akdim K, Zahid M. Qualitative analysis of a fractional-order two-strain epidemic model with vaccination and general non-monotonic incidence rate. Int J Dyn Control. (2023) 11:1532–43. doi: 10.1007/s40435-022-01083-4

PubMed Abstract | Crossref Full Text | Google Scholar

38. Huo J, Zhao H, Zhu L. The effect of vaccines on backward bifurcation in a fractional order HIV model. Nonlinear Anal. (2015) 26:289–305. doi: 10.1016/j.nonrwa.2015.05.014

Crossref Full Text | Google Scholar

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

Crossref Full Text | Google Scholar

40. Mordecai EA, Paaijmans KP, Johnson LR, Balzer C, Ben-Horin T, de Moor E, et al. Optimal temperature for malaria transmission is dramatically lower than previously predicted. Ecol Lett. (2013) 16:22–30. doi: 10.1111/ele.12015

PubMed Abstract | Crossref Full Text | Google Scholar

41. Abdulaziz YA, Mukhtar JB, Rachid Ouifki M. Assessing the role of climate factors on malaria transmission dynamics in South Sudan. Math Biosci. (2019) 310:13–23. doi: 10.1016/j.mbs.2019.01.002

PubMed Abstract | Crossref Full Text | Google Scholar

42. Yiga V, Nampala H, Tumwiine J. Analysis of the model on the effect of seasonal factors on malaria transmission dynamics. J Appl Math. (2020) 2020:1–19. doi: 10.1155/2020/8885558

Crossref Full Text | Google Scholar

43. Toufik M, Atangana A. New numerical approximation of fractional derivative with non-local and non-singular kernel: Application to chaotic models. Eur Phys J Plus. (2017) 132:444. doi: 10.1140/epjp/i2017-11717-0

Crossref Full Text | Google Scholar

44. Raza A, Rocha E, Fadhal E, Ibrahim RIH, Afkar E, Bilal M. The effect of delay techniques on a lassa fever epidemic model. Complexity. (2024) 2024:2075354. doi: 10.1155/2024/2075354

Crossref Full Text | Google Scholar

45. Raza A, Abdella K. Analysis of the dynamics of anthrax epidemic model with delay. Deleted J. (2024) 6:128. doi: 10.1007/s42452-024-05763-y

PubMed Abstract | Crossref Full Text | Google Scholar

46. Alfwzan WF, Raza A, Martin-Vaquero J, Baleanu D, Rafiq M, Ahmed N, et al. Modeling and transmission dynamics of Zika virus through efficient numerical method. AIP Adv. (2023) 13:8945. doi: 10.1063/5.0168945

Crossref Full Text | Google Scholar

47. Raza A. Mathematical modelling of rotavirus disease through efficient methods. Comp Mater Cont. (2022) 72:4727–40. doi: 10.32604/cmc.2022.027044

Crossref Full Text | Google Scholar

48. Deresse AT, Mussa YO, Gizaw AK. Approximate analytical solution of two-dimensional nonlinear time-fractional damped wave equation in the caputo fractional derivative operator. Math Probl Eng. (2022) 2022:1–28. doi: 10.1155/2022/7004412

Crossref Full Text | Google Scholar

49. Gizaw AK, Mussa YO. Approximate analytical solutions of two-dimensional time fractional Kleingordon Equation. Ethiop J Educ Sci. (2022) 16:1–16. Retrieved from: https://ejhs.ju.edu.et/index.php/ejes/article/view/3956

Google Scholar