Fractional Optimal Control of Anthroponotic Cutaneous Leishmaniasis with Behavioral and Epidemiological Extensions

1. Introduction

Cutaneous leishmaniasis (CL) is a sandfly-transmitted protozoal infection caused by Leishmania protozoa and transmitted principally by the bite of infected females of Phlebotomus sandflies. CL is a neglected disease, which disproportionately affects low-resource settings in the Middle East, South America, and South Asia, where CL causes skin ulcers commonly leading to permanent scarring and social discrimination [1]. Of all the forms of leishmaniasis, the anthroponotic form (ACL) is particularly epidemiologically challenging in that it is transmitted from person to person by vectors with humans as the primary reservoir [2]. Control is harder thus because interruption of transmission involves controlling human infectivity directly [3].

The CL transmission cycle is greatly influenced by ecological and behavioral factors such as sandfly population density, frequency of biting, and infection rates. Traditional vector control interventions such as insecticide-treated bed nets, indoor residual spraying, and environmental sanitation have, in certain environments, been effective. The increased resistance to insecticides among sandflies [4], combined with inadequate public awareness and behavioral resistance, significantly reduces the long-term efficacy of such interventions.

Mathematical modeling has become an essential tool for understanding disease dynamics and optimizing control strategy. Classical compartmental models (e.g., SIR and SEIR) have been widely utilized to represent vector-borne diseases and estimate transmission thresholds, evaluate control impacts, and forecast outbreak scenarios [5,6,7,8,9]. Such models are typically too limiting in application because they use integer-order differential equations with instantaneous transitions and no provision to incorporate biological memory or behavioral feedback—the key ingredients in real epidemiological systems.

Recent advances have witnessed the application of fractional-order differential equations to disease modeling. Fractional derivatives, especially in their Caputo form, are capable of more accurately modeling systems where memory and non-locality are important features, such as diseases with delayed onset, time-varying response to treatment, or long-term behavioral modification [10,11,12]. Fractional differential equations have, in particular, shown greater correspondence with empirical disease data in cases such as COVID-19 and dengue fever [13].

At the same time, inclusion of behavioral dynamics—awareness, risk perception, and intervention fatigue—in disease models has been under closer examination. Dynamizing awareness as a model parameter enables simulations derived from feedback in actual behaviors, especially in response to health campaigns, outbreak severity, and social adjustment [14,15,16].

Although optimal control theory has been paired with traditional models of disease in some research to assess resource-effective intervention strategies [17,18,19], few have done so in a fractional framework. Further, most models include none or at best one of these real-world complexities at a time—like asymptomatic transmission, resistance to insecticides, behavioral dynamics, and economic limitations—in a cohesive and biologically well-supported framework.

This study builds upon and extends the classical model presented in [20], which describes the transmission dynamics of ACL using integer-order differential equations and a limited set of compartments. Our development enhances this foundational framework by introducing three additional state variables—representing asymptomatic human infections $\left(I_a\right)$, insecticide-resistant vectors $\left(I_{vr}\right)$, and a dynamic awareness function $\left(A\right)$ -to better capture key epidemiological and behavioral complexities. Furthermore, we add a fourth control variable $\left(u_4\right)$ focused on public awareness campaigns, complementing existing interventions such as treatment, bed net usage, and vector spraying. The adoption of Caputo fractional derivatives incorporates memory effects into the model, enabling a more realistic portrayal of delayed biological responses and behavioral inertia. Through this multi-faceted extension, the model provides a more comprehensive and flexible tool for optimizing ACL control strategies under resource constraints.

The remainder of the paper is organized as follows. Section 2 presents the mathematical preliminaries about fractional calculus. Section 3 introduces the proposed fractional-order model for ACL, together with its properties and stability analysis. Section 4 defines the optimal control problem, incorporating the effects of multiple control strategies, and establishes the necessary optimality conditions. Section 5 develops a numerical method to solve the fractional model at hand along with its convergence analysis. Section 6 reports and discusses the numerical simulations, including scenario analysis and policy implications. Finally, Section 7 concludes the paper with a summary of key findings and recommendations for future research.

2. Mathematical Preliminaries

In this section, we recall essential definitions and properties from fractional calculus that will be used throughout the formulation and numerical analysis of our model. Standard references include [10,21,22,23].

Definition 1.

The left-sided Riemann–Liouville fractional integral of order $\alpha >0$ of a function $x\left(t\right)$ is defined by:

$${}^{R}I_t^{\alpha }x\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_0^t{(t-\tau )}^{\alpha -1}x\left(\tau \right)d\tau ,t>0,$$

(1)

where $\mathsf{\Gamma }(·)$ denotes the Gamma function.

Definition 2.

The left-sided Riemann-Liouville fractional derivative of order α of a function $x\left(t\right)$ is defined as: 

$${}^{R}D_t^{\alpha }x\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(n-\alpha )}}{\displaystyle \frac{d^n}{d{t}^n}}{\int }_0^t{(t-\tau )}^{n-\alpha -1}x\left(\tau \right)d\tau ,$$

(2)

where n is the smallest integer greater than or equal to α (i.e., $n=\left[\alpha \right]+1$).

Definition 3.

The left-sided Caputo derivative of order $\alpha \in (n-1,n)$, where $n\in \mathbb{N}$, is given by:

$${}^{C}D_t^{\alpha }x\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(n-\alpha )}}{\int }_0^t{(t-\tau )}^{n-\alpha -1}{x}^{\left(n\right)}\left(\tau \right)d\tau .$$

(3)

When $\alpha \in (0,1)$, we have:

$${}^{C}D_t^{\alpha }x\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }(1-\alpha )}}{\int }_0^t{(t-\tau )}^{-\alpha }{x}^{\prime }\left(\tau \right)d\tau .$$

(4)

Definition 4.

Let $\alpha >0$ and $0\le t<T$. The right-sided Riemann–Liouville fractional integral of order α is:

$${}^{R}I_t^{\alpha }x\left(t\right)={\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_t^T{(\tau -t)}^{\alpha -1}x\left(\tau \right)d\tau .$$

(5)

Definition 5.

Let $\alpha >0$ and $n=\lceil \alpha \rceil$. The right-sided Riemann–Liouville fractional derivative of order α is:

$${}^{R}D_t^{\alpha }x\left(t\right)={\displaystyle \frac{{(-1)}^n}{\mathsf{\Gamma }(n-\alpha )}}{\displaystyle \frac{d^n}{d{t}^n}}{\int }_t^T{(\tau -t)}^{n-\alpha -1}x\left(\tau \right)d\tau .$$

(6)

Definition 6.

Let $\alpha >0$ and $n=\lceil \alpha \rceil$. The right-sided Caputo fractional derivative of order α is:

$${}^{C}D_t^{\alpha }x\left(t\right)={\displaystyle \frac{{(-1)}^n}{\mathsf{\Gamma }(n-\alpha )}}{\int }_t^T{(\tau -t)}^{n-\alpha -1}{x}^{\left(n\right)}\left(\tau \right)d\tau .$$

(7)

Lemma 1.

Let $x\in C^n[0,T]$ and $n-1<\alpha <n$, then the following identity holds:

$${}^{R}I_t^{\alpha }{}^{C}D_t^{\alpha }x\left(t\right)=x\left(t\right)-\sum _{k=0}^{n-1}{\displaystyle \frac{x^{\left(k\right)}\left(0\right)}{k!}}{t}^k.$$

(8)

Key Properties:

Let $\alpha ,\beta >0$ and $x_1,x_2\in L^1[0,T]$. Then:

• Linearity: ${}^{C}D_t^{\alpha }(a{x}_1\left(t\right)+b{x}_2\left(t\right))=a{}^{C}D_t^{\alpha }{x}_1\left(t\right)+b{}^{C}D_t^{\alpha }{x}_2\left(t\right)$,
• Composition of integrals: ${}^{R}I_t^{\alpha }{}^{C}I_t^{\alpha }x\left(t\right)={}^{R}I_t^{\alpha +\beta }x\left(t\right)$,
• Commutativity: ${}^{R}I_t^{\alpha }{}^{R}I_t^{\beta }x\left(t\right)={}^{R}I_t^{\beta }{}^{C}I_t^{\alpha }x\left(t\right)$.

These fundamental tools enable us to define and analyze the fractional-order epidemic model, derive the associated adjoint system, and construct the numerical scheme for solving the optimal control problem.

3. Mathematical Structure of the Model

In this section, we introduce a new model formulation for the ACL and investigate its properties and stability analysis.

3.1. Formulation of the Anthroponotic Cutaneous Leishmaniasis (ACL) Model

Fractional calculus, a specialized branch of mathematical analysis, focuses on the study and application of fractional derivatives and integrals. This field has been instrumental in modeling complex systems across various disciplines including engineering, physics, and biology. Fractional models are particularly valuable due to their ability to incorporate memory effects, enabling more accurate estimation and prediction of dynamic phenomena. Given these advantageous properties of fractional mathematical modeling, and considering that the previously presented model in [20] does not account for memory effects, we extend the framework by employing fractional derivatives in a general sense to develop an enhanced transmission model. Accordingly, the novel model is formulated as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{S}_h\left(t\right)& ={\mathsf{\Lambda }}_h-{\lambda }_h{S}_h-{\mu }_h{S}_h,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{E}_h\left(t\right)& ={\lambda }_h{S}_h-(k_1+\theta +{\mu }_h)E_h,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{I}_h\left(t\right)& =k_1(1-p)E_h-(\gamma +{\mu }_h)I_h,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{I}_a\left(t\right)& =k_{1}p{E}_h-({\gamma }_a+{\mu }_h)I_a,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{R}_h\left(t\right)& =\theta E_h+\gamma I_h+{\gamma }_a{I}_a-{\mu }_h{R}_h,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{S}_v\left(t\right)& ={\mathsf{\Lambda }}_v-{\lambda }_v{S}_v-{\mu }_v{S}_v,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{E}_v\left(t\right)& ={\lambda }_v{S}_v-(k_2+{\mu }_v)E_v,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{I}_v\left(t\right)& =(1-q)k_2{E}_v-{\mu }_v{I}_v,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{I}_{v_r}\left(t\right)& =q{k}_2{E}_v-{\mu }_r{I}_{v_r},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }A\left(t\right)& =-\delta A.\hfill \end{array}$$

(9)

In this new framework, the model has been augmented by introducing additional state variables: asymptomatic human infections $\left(I_a\right)$, an insecticide-resistant vector population $\left(I_{vr}\right)$, and a dynamic awareness function within the human population $\left(A\right)$. Additionally, the auxiliary parameter $\sigma$ is introduced to guarantee that both sides of the fractional equations in equation (9) have consistent dimensions. In model (9), the forces of infection are given by:

$${\lambda }_h\left(t\right)=a{b}_1{\displaystyle \frac{I_v+\xi I_{v_r}}{N_v}},{\lambda }_v\left(t\right)=a{c}_1{\displaystyle \frac{I_h+\eta I_a}{N_h}},$$

(10)

where a is the biting rate, $b_1$ and $c_1$ are transmission probabilities, and $\xi ,\eta$ are relative infectivities. Table 1 presents the state variables. Model parameter values are given in Table 2.

3.2. Positivity of Solutions

To ensure biological consistency of the proposed fractional-order model (9), we establish that its solutions remain non-negative for all time $t\in [0,T]$, provided the initial conditions are non-negative. Let:

$$\mathcal{X}\left(t\right)=(S_h,E_h,I_h,I_a,R_h,S_v,E_v,I_v,I_{v_r},A),$$

(11)

be the vector of state variables, and assume $\mathcal{X}\left(0\right)\ge 0$. Define:

$$W\left(t\right):=min\{S_h\left(t\right),E_h\left(t\right),I_h\left(t\right),I_a\left(t\right),R_h\left(t\right),S_v\left(t\right),E_v\left(t\right),I_v\left(t\right),I_{v_r}\left(t\right),A\left(t\right)\}.$$

(12)

Suppose, by contradiction, that there exists a first time $t_1>0$ such that $W\left(t_1\right)=0$ while $W\left(t\right)>0$ for all $t\in [0,t_1)$. Without loss of generality, let $W\left(t_1\right)=S_h\left(t_1\right)=0$. From (9), the susceptible human population satisfies:

$${\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{S}_h\left(t\right)={\mathsf{\Lambda }}_h-{\lambda }_h{S}_h-{\mu }_h{S}_h,$$

(13)

where ${\lambda }_h\left(t\right)\ge 0$ is the force of infection. Multiplying both sides by ${\sigma }^{1-\alpha }>0$ gives:

$${}^{C}D_t^{\alpha }{S}_h\left(t\right)={\sigma }^{1-\alpha }\left[{\mathsf{\Lambda }}_h-K{S}_h\right],K:={\lambda }_h+{\mu }_h\ge 0.$$

(14)

Thus:

$${}^{C}D_t^{\alpha }{S}_h\left(t\right)\ge {\sigma }^{1-\alpha }\left[{\mathsf{\Lambda }}_h-K{S}_h\right].$$

(15)

Since ${\sigma }^{1-\alpha }>0$, dividing by it preserves the inequality sign, so we can write the equivalent scalar inequality:

$${}^{C}D_t^{\alpha }{S}_h\left(t\right)\ge {\sigma }^{1-\alpha }{\mathsf{\Lambda }}_h-{\sigma }^{1-\alpha }K{S}_h,S_h\left(0\right)=S_h\left(0\right)>0.$$

(16)

The solution of the corresponding linear Caputo fractional differential equation:

$${}^{C}D_t^{\alpha }{S}_h\left(t\right)+{\sigma }^{1-\alpha }K{S}_h={\sigma }^{1-\alpha }{\mathsf{\Lambda }}_h,$$

(17)

remains positive for all $t>0$ (see [21]), hence $S_h\left(t\right)>0$ on $[0,t_1]$, contradicting $S_h\left(t_1\right)=0$. The same reasoning applies to all other compartments because:

• The factor ${\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}$ is positive and thus does not alter the sign structure of the system.
• All right-hand sides contain non-negative inflows (${\mathsf{\Lambda }}_h$, ${\mathsf{\Lambda }}_v$, or transfers from other non-negative compartments).
• All loss terms are proportional to the variable itself (linear decay or bilinear incidence), so they cannot produce negativity.

For $A\left(t\right)$:

$${\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }A\left(t\right)=-\delta A,A\left(0\right)\ge 0,$$

(18)

which implies $A\left(t\right)\ge 0$ by the fractional comparison principle. Therefore:

$$\mathcal{X}\left(0\right)\ge 0\Rightarrow \mathcal{X}\left(t\right)\ge 0,\forall t\in [0,T],$$

(19)

and the system remains biologically well-posed.

3.3. Boundedness of Solutions

To establish the uniform boundedness of solutions of system (9), we consider the total human and vector populations defined by:

$$N_h=S_h+E_h+I_h+I_a+R_h,N_v=S_v+E_v+I_v+I_{v_r}.$$

(20)

Summing the fractional equations for the human compartments in (9) gives:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{N}_h\left(t\right)& ={\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }\left(S_h+E_h+I_h+I_a+R_h\right)\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & ={\mathsf{\Lambda }}_h-{\mu }_h(S_h+E_h+I_h+I_a+R_h)\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill & ={\mathsf{\Lambda }}_h-{\mu }_h{N}_h.\hfill \end{array}$$

(23)

Multiplying both sides by ${\sigma }^{1-{\alpha }_1}>0$, we get:

$${}^{C}D_t^{\alpha }{N}_h\left(t\right)={\sigma }^{1-\alpha }{\mathsf{\Lambda }}_h-{\sigma }^{1-\alpha }{\mu }_h{N}_h.$$

(24)

Applying the fractional comparison principle (see [10]), the scalar fractional differential equation:

$${}^{C}D_t^{\alpha }{N}_h\left(t\right)={\sigma }^{1-\alpha }{\mathsf{\Lambda }}_h-{\sigma }^{1-\alpha }{\mu }_h{N}_h,N_h\left(0\right)=N_h\left(0\right),$$

(25)

has the unique positive bounded solution:

$$\underset{t\to \infty }{lim\; sup}{N}_h\left(t\right)={\displaystyle \frac{{\mathsf{\Lambda }}_h}{{\mu }_h}}:=M_h,$$

(26)

which implies:

$$\underset{t\to \infty }{lim\; sup}{N}_h\left(t\right)\le M_h.$$

(27)

Similarly, summing the vector equations in (9) without controls gives:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{N}_v\left(t\right)& ={\mathsf{\Lambda }}_v-{\mu }_v(S_v+E_v+I_v)-{\mu }_r{I}_{v_r}\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill & ={\mathsf{\Lambda }}_v-min\{{\mu }_v,{\mu }_r\}{N}_v.\hfill \end{array}$$

(29)

Multiplying both sides by ${\sigma }^{1-\alpha }$ yields:

$${}^{C}D_t^{\alpha }{N}_v\left(t\right)={\sigma }^{1-\alpha }{\mathsf{\Lambda }}_v-{\sigma }^{1-\alpha }min\{{\mu }_v,{\mu }_r\}{N}_v,$$

(30)

and by the fractional comparison principle:

$$\underset{t\to \infty }{lim\; sup}{N}_v\left(t\right)\le {\displaystyle \frac{{\mathsf{\Lambda }}_v}{min\{{\mu }_v,{\mu }_r\}}}:=M_v.$$

(31)

Finally, the awareness function $A\left(t\right)$ satisfies:

$${\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }A\left(t\right)=-\delta A.$$

(32)

Multiplying both sides by ${\sigma }^{1-\alpha }$, we obtain:

$${}^{C}D_t^{\alpha }A\left(t\right)=-{\sigma }^{1-\alpha }\delta A,$$

(33)

whose unique solution is positive and bounded for all $t\ge 0$, monotonically decreases to 0 as $t\to \infty$ (see [21,26]), and tends to zero as $t\to \infty$:

$$\underset{t\to \infty }{lim\; sup}A\left(t\right)=0.$$

(34)

Therefore, the solutions evolve in the positively invariant compact region:

$$\mathsf{\Omega }:=\left\{\mathcal{X}\left(t\right)\in {\mathbb{R}}_{+}^{10}|N_h\left(t\right)\le M_h,N_v\left(t\right)\le M_v,A\left(t\right)\ge 0\right\},$$

(35)

which guarantees the global boundedness of all state variables on $[0,T]$.

3.4. Existence of Disease-Free Equilibrium (DFE)

We analyze the existence of a disease-free equilibrium (DFE) in the fractional-order ACL model introduced in Section 3.1. The DFE represents a biologically meaningful state in which no individuals are infected, and all disease-related compartments vanish. In our case, the DFE is defined by:

$${\mathcal{E}}^0=(S_h^0,E_h^0=0,I_h^0=0,I_a^0=0,R_h^0,S_v^0,E_v^0=0,I_v^0=0,I_{v_r}^0=0,A^0).$$

(36)

Substituting these values into the fractional system, the equilibrium conditions reduce to the following algebraic equations:

$$\begin{array}{cc}\hfill 0& ={\mathsf{\Lambda }}_h-{\mu }_h{S}_h^0,\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill 0& =\theta \times 0+\gamma \times 0+{\gamma }_a\times 0-{\mu }_h{R}_h^0,\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill 0& ={\mathsf{\Lambda }}_v-{\mu }_v{S}_v^0,\hfill \end{array}$$

(39)

$$\begin{array}{cc}\hfill 0& =-\delta A^0.\hfill \end{array}$$

(40)

Solving these, we obtain the DFE values:

$$\begin{array}{cc}\hfill S_h^0& ={\displaystyle \frac{{\mathsf{\Lambda }}_h}{{\mu }_h}},\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill R_h^0& =0,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill S_v^0& ={\displaystyle \frac{{\mathsf{\Lambda }}_v}{{\mu }_v}},\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill A^0& =0.\hfill \end{array}$$

(44)

All other compartments related to exposed and infected individuals are zero in the DFE by definition. Hence, the disease-free steady state is given by:

$${\mathcal{E}}^0=\left({\displaystyle \frac{{\mathsf{\Lambda }}_h}{{\mu }_h}},0,0,0,0,{\displaystyle \frac{{\mathsf{\Lambda }}_v}{{\mu }_v}},0,0,0,0\right).$$

(45)

3.5. Basic Reproduction Number ${\mathcal{R}}_0$

To quantify the initial transmission potential of ACL in the absence of acquired immunity, we derive the basic reproduction number ${\mathcal{R}}_0$ using the next-generation matrix approach [6]. In this method, we first identify the set of infected compartments as:

$${\mathcal{X}}_I\left(t\right)=\left(E_h,I_h,I_a,E_v,I_v,I_{v_r}\right).$$

(46)

At the DFE, $N_h$ and $N_v$ are constants and $S_h$ and $S_v$ equal their equilibrium values, hence ${\lambda }_h$ and ${\lambda }_v$ become linear in the infected variables with constant coefficients. The vector of new infection terms $\mathcal{F}$ and the vector of transition terms $\mathcal{V}$ are:

$$\mathcal{F}=\left(\begin{array}{c}{\lambda }_h{S}_h\\ 0\\ 0\\ {\lambda }_v{S}_v\\ 0\\ 0\end{array}\right),\mathcal{V}=\left(\begin{array}{c}(k_1+\theta +{\mu }_h)E_h\\ -k_1(1-p)E_h+(\gamma +{\mu }_h)I_h\\ -k_{1}p{E}_h+({\gamma }_a+{\mu }_h)I_a\\ (k_2+{\mu }_v)E_v\\ -(1-q)k_2{E}_v+{\mu }_v{I}_v\\ -q{k}_2{E}_v+{\mu }_r{I}_{v_r}\end{array}\right).$$

(47)

At the DFE we have:

$$S_h^0={\displaystyle \frac{{\mathsf{\Lambda }}_h}{{\mu }_h}},S_v^0={\displaystyle \frac{{\mathsf{\Lambda }}_v}{{\mu }_v}},N_h^0=S_h^0,N_v^0=S_v^0,$$

(48)

and all infected compartments are zero. Substituting these values into ${\lambda }_h$ and ${\lambda }_v$ yields:

$${\lambda }_h{S}_h^0⟶a{b}_1{\displaystyle \frac{S_h^0}{N_v^0}}\left(I_v+\xi I_{vr}\right),{\lambda }_v{S}_v^0⟶a{c}_1{\displaystyle \frac{S_v^0}{N_h^0}}\left(I_h+\eta I_a\right),$$

(49)

which are linear functions of the infected variables with constant coefficients at the DFE. Let $\mathbf{F}$ and $\mathbf{V}$ denote the Jacobian matrices of $\mathcal{F}$ and $\mathcal{V}$ with respect to ${\mathcal{X}}_I$, evaluated at the DFE. The next-generation matrix is given by:

$$\mathbf{K}=\mathbf{F}{\mathbf{V}}^{-1},$$

(50)

and the basic reproduction number is defined as the spectral radius of $\mathbf{K}$:

$${\mathcal{R}}_0=\rho \left(\mathbf{K}\right).$$

(51)

Carrying out the algebra and substituting $S_h^0$ and $S_v^0$ leads to the explicit expression:

$${\mathcal{R}}_0=\sqrt{\left({\displaystyle \frac{a^2{b}_1{c}_1{\mathsf{\Lambda }}_h{\mathsf{\Lambda }}_v}{{\mu }_h{\mu }_v^2}}\right)\times \left({\displaystyle \frac{k_1}{k_1+\theta +{\mu }_h}}\right)\times \left[{\displaystyle \frac{1-p}{\gamma +{\mu }_h}}+{\displaystyle \frac{p}{{\gamma }_a+{\mu }_h}}\right]\times \left[{\displaystyle \frac{1-q}{{\mu }_v}}+{\displaystyle \frac{q}{{\mu }_r}}\right]}.$$

(52)

This value represents the expected number of secondary human infections generated by a single infected human in a fully susceptible population, taking into account transmission through symptomatic, asymptomatic, and resistant vector pathways. When ${\mathcal{R}}_0<1$ the disease tends to die out, while ${\mathcal{R}}_0>1$ indicates the possibility of sustained transmission in the absence of interventions.

3.6. Existence of Endemic Equilibrium (EE)

An endemic equilibrium (EE) corresponds to a steady-state solution of the model (9), where the disease persists in the population. This implies that at least one of the infectious compartments satisfies:

$$I_h^{+}>0,I_a^{+}>0,I_v^{+}>0,\mathrm{or}{I}_{v_r}^{+}>0.$$

(53)

Let:

$${\mathcal{E}}^{+}=\left(S_h^{+},E_h^{+},I_h^{+},I_a^{+},R_h^{+},S_v^{+},E_v^{+},I_v^{+},I_{v_r}^{+},A^{+}\right),$$

(54)

be the EE point. Since the system is governed by Caputo fractional derivatives, at equilibrium all time derivatives vanish, and the system reduces to the following algebraic equations:

$$\begin{array}{cc}\hfill 0& ={\mathsf{\Lambda }}_h-{\lambda }_h^{+}{S}_h^{+}-{\mu }_h{S}_h^{+},\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill 0& ={\lambda }_h^{+}{S}_h^{+}-(k_1+\theta +{\mu }_h)E_h^{+},\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill 0& =k_1(1-p)E_h^{+}-(\gamma +{\mu }_h)I_h^{+},\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill 0& =k_{1}p{E}_h^{+}-({\gamma }_a+{\mu }_h)I_a^{+},\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill 0& =\theta E_h^{+}+\gamma I_h^{+}+{\gamma }_a{I}_a^{+}-{\mu }_h{R}_h^{+},\hfill \end{array}$$

(59)

$$\begin{array}{cc}\hfill 0& ={\mathsf{\Lambda }}_v-{\lambda }_v^{+}{S}_v^{+}-{\mu }_v{S}_v^{+},\hfill \end{array}$$

(60)

$$\begin{array}{cc}\hfill 0& ={\lambda }_v^{+}{S}_v^{+}-(k_2+{\mu }_v)E_v^{+},\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill 0& =(1-q)k_2{E}_v^{+}-{\mu }_v{I}_v^{+},\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill 0& =q{k}_2{E}_v^{+}-{\mu }_r{I}_{v_r}^{+},\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill 0& =-\delta A^{+}.\hfill \end{array}$$

(64)

The forces of infection at equilibrium are defined by:

$$\begin{array}{cc}\hfill {\lambda }_h^{+}& =a(1-A^{+})b_1{\displaystyle \frac{I_v^{+}+I_{v_r}^{+}}{N_v^{+}}},\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill {\lambda }_v^{+}& =a(1-A^{+})c_1{\displaystyle \frac{I_h^{+}+I_a^{+}}{N_h^{+}}},\hfill \end{array}$$

(66)

where:

$$N_h^{+}=S_h^{+}+E_h^{+}+I_h^{+}+I_a^{+}+R_h^{+},N_v^{+}=S_v^{+}+E_v^{+}+I_v^{+}+I_{v_r}^{+}.$$

(67)

From the last equation, it follows that:

$$A^{+}=0.$$

(68)

Substituting the third and fourth equations into the system gives:

$$I_h^{+}={\displaystyle \frac{k_1(1-p)}{\gamma +{\mu }_h}}{E}_h^{+},I_a^{+}={\displaystyle \frac{k_{1}p}{{\gamma }_a+{\mu }_h}}{E}_h^{+}.$$

(69)

From the fifth equation, we get:

$$R_h^{+}={\displaystyle \frac{1}{{\mu }_h}}\left[\theta +{\displaystyle \frac{k_1(1-p)\gamma }{\gamma +{\mu }_h}}+{\displaystyle \frac{k_{1}p{\gamma }_a}{{\gamma }_a+{\mu }_h}}\right]E_h^{+}.$$

(70)

From the eighth and ninth equations, we also obtain:

$$I_v^{+}={\displaystyle \frac{(1-q)k_2}{{\mu }_v}}{E}_v^{+},I_{v_r}^{+}={\displaystyle \frac{q{k}_2}{{\mu }_r}}{E}_v^{+}.$$

(71)

Finally, from the first and sixth equations, the susceptible compartments are:

$$S_h^{+}={\displaystyle \frac{{\mathsf{\Lambda }}_h}{{\lambda }_h^{+}+{\mu }_h}},S_v^{+}={\displaystyle \frac{{\mathsf{\Lambda }}_v}{{\lambda }_v^{+}+{\mu }_v}}.$$

(72)

By substituting the above expressions into the remaining equations, the system reduces to a set of nonlinear algebraic equations in the two variables $E_h^{+}$ and $E_v^{+}$, which can be solved numerically given fixed parameter values. The model admits a biologically meaningful EE ${\mathcal{E}}^{+}$ provided that the basic reproduction number satisfies ${\mathcal{R}}_0>1$. The explicit dependence of the forces of infection on equilibrium state variables plays a key role in determining the feasibility and stability of ${\mathcal{E}}^{+}$.

3.7. Stability Analysis of Equilibria

In this section, we investigate the local stability conditions for both the DFE and the EE of the uncontrolled system.

3.7.1. Stability of the DFE

We consider the DFE given by:

$${\mathcal{E}}^0=\left(S_h^0={\displaystyle \frac{{\mathsf{\Lambda }}_h}{{\mu }_h}},E_h^0=0,I_h^0=0,I_a^0=0,R_h^0=0,S_v^0={\displaystyle \frac{{\mathsf{\Lambda }}_v}{{\mu }_v}},E_v^0=0,I_v^0=0,I_{v_r}^0=0,A^0=0\right).$$

(73)

At this point, all infected compartments are zero. We linearize the model (9) around ${\mathcal{E}}^0$ and compute the Jacobian matrix. Following the standard theory of fractional-order differential equations [6,10,26], the DFE is locally asymptotically stable if the following fractional stability condition holds:

$$|arg\left(p_i\right)|>{\displaystyle \frac{\alpha \pi }{2}},\forall i,$$

(74)

where $p_i$ are the eigenvalues of the Jacobian matrix evaluated at ${\mathcal{E}}^0$, and $\alpha \in (0,1]$ is the order of the Caputo derivative. We also compute the basic reproduction number ${\mathcal{R}}_0$ via the next-generation matrix method. The result is:

$${\mathcal{R}}_0\approx 2.14>1.$$

(75)

Therefore, under the given parameters and no control applied, ${\mathcal{E}}^0$ is unstable. If ${\mathcal{R}}_0<1$, then the DFE ${\mathcal{E}}^0$ is locally asymptotically stable for any fractional order $\alpha \in (0,1]$. If ${\mathcal{R}}_0>1$, then ${\mathcal{E}}^0$ is unstable.

3.7.2. Stability of the EE

Let the EE be:

$${\mathcal{E}}^{+}=(S_h^{+},E_h^{+},I_h^{+},I_a^{+},R_h^{+},S_v^{+},E_v^{+},I_v^{+},I_{v_r}^{+},A^{+}),$$

(76)

where at least one of the infected compartments is positive. To assess the local stability of ${\mathcal{E}}^{+}$, we numerically compute the Jacobian matrix $J\left({\mathcal{E}}^{+}\right)$ using equilibrium values obtained from solving the steady-state system. The fractional-order stability criterion requires:

$$|arg\left(p_i\right)|>{\displaystyle \frac{\alpha \pi }{2}}\mathrm{for}\mathrm{all}\mathrm{eigenvalues}{p}_i\mathrm{of}J\left({\mathcal{E}}^{+}\right).$$

(77)

Assuming the model parameters in Table 2, and setting $\alpha =0.95$, we compute the eigenvalues at ${\mathcal{E}}^{+}$ numerically:

$$p_1=-0.115,p_2=-0.207,\dots ,p_{10}=-0.992.$$

(78)

All eigenvalues are real and negative, so $arg\left(p_i\right)=\pi$, and hence the stability condition:

$${\displaystyle \frac{\alpha \pi }{2}}=0.475\pi <\pi =|arg\left(p_i\right)|\Rightarrow \mathrm{holds}\mathrm{for}\mathrm{all}i.$$

(79)

The EE point ${\mathcal{E}}^{+}$ is locally asymptotically stable for the given parameter values and $\alpha \in (0,1]$, since all eigenvalues satisfy the fractional stability condition.

3.8. Sensitivity Analysis

To identify the most influential parameters affecting the basic reproduction number ${\mathcal{R}}_0$, a global sensitivity analysis was conducted using the Partial Rank Correlation Coefficient (PRCC) method [27,28]. This technique assesses the degree of monotonic relationship between each model parameter and the output variable ${\mathcal{R}}_0$, while controlling for the effects of other parameters [29,30,31]. Figure 1 illustrates the PRCC values obtained from 1,000 Latin Hypercube Sampling (LHS) runs, and Table 3 provides the corresponding numerical values. Parameters with larger absolute PRCC values (positive or negative) have a stronger impact on ${\mathcal{R}}_0$. The most sensitive parameters were found to be the average biting rate a (PRCC = 0.5916), vector-to-human transmission probability $b_1$ (0.2454), human recruitment rate ${\mathsf{\Lambda }}_h$ (0.2369), and human-to-vector transmission probability $c_1$ (0.2126). These positive PRCC values indicate that increasing any of these parameters leads to an increase in ${\mathcal{R}}_0$, thus enhancing disease transmission potential. On the other hand, the most significant negative influences are the vector natural death rate ${\mu }_v$ (PRCC = $-0.5305$) and human death rate ${\mu }_h$ ($-0.3806$). An increase in these parameters leads to a reduction in ${\mathcal{R}}_0$, indicating a suppressive effect on disease spread. Parameters such as p, q, and ${\mu }_r$ have very low PRCC values close to zero, suggesting minimal influence on ${\mathcal{R}}_0$ under the explored parameter ranges. This insight can guide prioritization in control strategies by targeting parameters that contribute most to disease persistence.

Overall, the PRCC analysis demonstrates that transmission-related parameters and demographic rates play a dominant role in shaping disease dynamics and should be key targets in designing intervention policies.

4. Objective Functional and Constraints

The controlled fractional-order model is given by:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{S}_h\left(t\right)& ={\mathsf{\Lambda }}_h-{\lambda }_h{S}_h-{\mu }_h{S}_h,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{E}_h\left(t\right)& ={\lambda }_h{S}_h-(k_1+\theta +{\mu }_h)E_h,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{I}_h\left(t\right)& =k_1(1-p)E_h-(\gamma +u_2+{\mu }_h)I_h,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{I}_a\left(t\right)& =k_{1}p{E}_h-({\gamma }_a+{\mu }_h)I_a,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{R}_h\left(t\right)& =\theta E_h+\gamma I_h+{\gamma }_a{I}_a-{\mu }_h{R}_h,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{S}_v\left(t\right)& ={\mathsf{\Lambda }}_v-{\lambda }_v{S}_v-({\mu }_v+u_3)S_v,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{E}_v\left(t\right)& ={\lambda }_v{S}_v-(k_2+{\mu }_v+u_3)E_v,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{I}_v\left(t\right)& =(1-q)k_2{E}_v-({\mu }_v+u_3)I_v,\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }{I}_{v_r}\left(t\right)& =q{k}_2{E}_v-({\mu }_r+u_3)I_{v_r},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }A\left(t\right)& =\eta u_4-\delta A.\hfill \end{array}$$

(80)

The control variables are defined in Table 4. The forces of infection are taken from the transmission terms in Section 3.5, modified to incorporate control $u_1\left(t\right)$:

$${\lambda }_h\left(t\right)=(1-u_1)a{b}_1{\displaystyle \frac{I_v+\xi I_{v_r}}{N_v}},{\lambda }_v\left(t\right)=(1-u_1)a{c}_1{\displaystyle \frac{I_h+\eta I_a}{N_h}},$$

(81)

where $N_h=S_h+E_h+I_h+I_a+R_h$ and $N_v=S_v+E_v+I_v+I_{v_r}$. The objective is to determine optimal controls $\{u_1,u_2,u_3,u_4\}$ over a finite time horizon $[0,T]$ that minimize both the number of infected and exposed individuals and the implementation costs. The performance index is:

$$\begin{array}{cc}\hfill J(u_1,u_2,u_3,u_4)={\int }_0^T& [w_1{I}_h+w_2{I}_a+w_3(I_v+I_{v_r})+w_4{E}_h\hfill \\ \hfill & +\frac{1}{2}(c_1{u}_1^2+c_2{u}_2^2+c_3{u}_3^2+c_4{u}_4^2)]dt,\hfill \end{array}$$

(82)

where $w_i$ ($i=1,\dots ,4$) are weights for disease burden and $c_i$ are cost coefficients. Controls satisfy the bounds:

$$0\le u_i\left(t\right)\le 1,i=1,\dots ,4,t\in [0,T].$$

(83)

4.1. Pontryagin’s Maximum Principle

To derive necessary conditions for optimality in the fractional-order control problem, we apply Pontryagin’s maximum principle adapted for systems governed by Caputo derivatives. This principle introduces adjoint variables and defines a Hamiltonian function to characterize the optimal controls. The Hamiltonian function $\mathcal{H}$ corresponding to the system (80) and the objective functional is given by:

$$\begin{array}{cc}\hfill \mathcal{H}=& w_1{I}_h+w_2{I}_a+w_3(I_v+I_{v_r})+w_4{E}_h\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\left(c_1{u}_1^2+c_2{u}_2^2+c_3{u}_3^2+c_4{u}_4^2\right)\hfill \\ \hfill & +{\lambda }_1\left[{\mathsf{\Lambda }}_h-{\lambda }_h{S}_h-{\mu }_h{S}_h\right]\hfill \\ \hfill & +{\lambda }_2\left[{\lambda }_h{S}_h-(k_1+\theta +{\mu }_h)E_h\right]\hfill \\ \hfill & +{\lambda }_3\left[k_1(1-p)E_h-(\gamma +u_2+{\mu }_h)I_h\right]\hfill \\ \hfill & +{\lambda }_4\left[k_{1}p{E}_h-({\gamma }_a+{\mu }_h)I_a\right]\hfill \\ \hfill & +{\lambda }_5\left[\theta E_h+\gamma I_h+{\gamma }_a{I}_a-{\mu }_h{R}_h\right]\hfill \\ \hfill & +{\lambda }_6\left[{\mathsf{\Lambda }}_v-{\lambda }_v{S}_v-({\mu }_v+u_3)S_v\right]\hfill \\ \hfill & +{\lambda }_7\left[{\lambda }_v{S}_v-(k_2+{\mu }_v+u_3)E_v\right]\hfill \\ \hfill & +{\lambda }_8\left[(1-q)k_2{E}_v-({\mu }_v+u_3)I_v\right]\hfill \\ \hfill & +{\lambda }_9\left[q{k}_2{E}_v-({\mu }_r+u_3)I_{v_r}\right]\hfill \\ \hfill & +{\lambda }_{10}\left[\eta u_4-\delta A\right].\hfill \end{array}$$

(84)

Let ${\lambda }_i\left(t\right)$, for $i=1,\dots ,10$, denote the adjoint variables associated with the state variables $\{S_h,E_h,I_h,I_a,R_h,S_v,E_v,I_v,I_{v_r},A\}$, respectively. These satisfy the system of right-sided Riemann–Liouville fractional differential equations:

$${\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_i\left(t\right)={\displaystyle \frac{\partial \mathcal{H}}{\partial x_i}},{\lambda }_i\left(T\right)=0,i=1,\dots ,10,$$

(85)

where:

$$\begin{array}{cccccccccc}\hfill x_1\left(t\right)& =S_h\left(t\right),\hfill & \hfill x_2\left(t\right)& =E_h\left(t\right),\hfill & \hfill x_3\left(t\right)& =I_h\left(t\right),\hfill & \hfill x_4\left(t\right)& =I_a\left(t\right),\hfill & \hfill x_5\left(t\right)& =R_a\left(t\right),\hfill \end{array}$$

(86)

$$\begin{array}{cccccccccc}\hfill x_6\left(t\right)& =S_{\nu }\left(t\right),\hfill & \hfill x_7\left(t\right)& =E_{\nu }\left(t\right),\hfill & \hfill x_8\left(t\right)& =I_{\nu }\left(t\right),\hfill & \hfill x_9\left(t\right)& =I_{{\nu }_r}\left(t\right),\hfill & \hfill x_{10}\left(t\right)& =A\left(t\right).\hfill \end{array}$$

(87)

The explicit adjoint equations are:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_1\left(t\right)& =-{\lambda }_1{\mu }_h+({\lambda }_2-{\lambda }_1){\lambda }_h-({\lambda }_7-{\lambda }_6){\beta }_v{\displaystyle \frac{(I_h+\eta I_a)S_v}{N_h^2}},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_2\left(t\right)& =w_4-{\lambda }_2(k_1+\theta +{\mu }_h)+{\lambda }_3{k}_1(1-p)+{\lambda }_4{k}_{1}p+{\lambda }_5\theta \hfill \\ \hfill & -({\lambda }_7-{\lambda }_6){\beta }_v{\displaystyle \frac{(I_h+\eta I_a)S_v}{N_h^2}},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_3\left(t\right)& =w_1-{\lambda }_3(\gamma +u_2+{\mu }_h)+{\lambda }_5\gamma \hfill \\ \hfill & +({\lambda }_7-{\lambda }_6){\beta }_v{\displaystyle \frac{\left(N_h-(I_h+\eta I_a)\right)S_v}{N_h^2}},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_4\left(t\right)& =w_2-{\lambda }_4({\gamma }_a+{\mu }_h)+{\lambda }_5{\gamma }_a\hfill \\ \hfill & +({\lambda }_7-{\lambda }_6){\beta }_v{\displaystyle \frac{\left(\eta N_h-(I_h+\eta I_a)\right)S_v}{N_h^2}},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_5\left(t\right)& =-{\lambda }_5{\mu }_h-({\lambda }_7-{\lambda }_6){\beta }_v{\displaystyle \frac{(I_h+\eta I_a)S_v}{N_h^2}},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_6\left(t\right)& =-{\lambda }_6({\lambda }_v+{\mu }_v+u_3)+{\lambda }_7{\lambda }_v\hfill \\ \hfill & -({\lambda }_2-{\lambda }_1)S_h{\beta }_h{\displaystyle \frac{(I_v+\xi I_{vr})}{N_v^2}},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_7\left(t\right)& =-{\lambda }_7(k_2+{\mu }_v+u_3)+{\lambda }_8{k}_2(1-q)+{\lambda }_9{k}_{2}q\hfill \\ \hfill & -({\lambda }_2-{\lambda }_1)S_h{\beta }_h{\displaystyle \frac{(I_v+\xi I_{vr})}{N_v^2}},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_8\left(t\right)& =w_3-{\lambda }_8({\mu }_v+u_3)+({\lambda }_2-{\lambda }_1)S_h{\beta }_h{\displaystyle \frac{N_v-(I_v+\xi I_{vr})}{N_v^2}},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_9\left(t\right)& =w_3\xi -{\lambda }_9({\mu }_r+u_3)+({\lambda }_2-{\lambda }_1)S_h{\beta }_h{\displaystyle \frac{\xi N_v-(I_v+\xi I_{vr})}{N_v^2}},\hfill \\ \hfill {\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{\lambda }_{10}\left(t\right)& =-{\lambda }_{10}\delta ,\hfill \end{array}$$

(88)

with terminal conditions:

$${\lambda }_i\left(T\right)=0,i=1,\dots ,10.$$

(89)

The optimal control functions $u_i^{*}\left(t\right)$, for $i=1,2,3,4$, minimize the Hamiltonian $\mathcal{H}$ pointwise and satisfy the control constraints $0\le u_i\left(t\right)\le 1$. They are characterized by:

$$u_i^{*}\left(t\right)=min\left\{max\left\{0,-{\displaystyle \frac{1}{c_i}}{\displaystyle \frac{\partial \mathcal{H}}{\partial u_i}}\right\},1\right\},t\in [0,T],$$

(90)

where the partial derivatives of the Hamiltonian with respect to the controls are given by:

$${\displaystyle \frac{\partial \mathcal{H}}{\partial u_i}}=\left\{\begin{array}{cc}{\displaystyle c_1{u}_1-({\lambda }_2-{\lambda }_1)a{b}_1{\displaystyle \frac{I_v+\xi I_{v_r}}{N_v}}{S}_h-({\lambda }_7-{\lambda }_6)a{c}_1{\displaystyle \frac{I_h+\eta I_a}{N_h}}{S}_v,}\hfill & i=1,\hfill \\ {\displaystyle c_2{u}_2-{\lambda }_3{I}_h,}\hfill & i=2,\hfill \\ {\displaystyle c_3{u}_3-{\lambda }_6{S}_v-{\lambda }_7{E}_v-{\lambda }_8{I}_v-{\lambda }_9{I}_{v_r},}\hfill & i=3,\hfill \\ {\displaystyle c_4{u}_4+{\lambda }_{10}\eta ,}\hfill & i=4.\hfill \end{array}\right.$$

(91)

The method for obtaining the partial derivatives of the Hamiltonian function with respect to each costate variable is as follows:

$${\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{R}D_t^{\alpha }{x}_i\left(t\right)={\displaystyle \frac{\partial \mathcal{H}}{\partial {\lambda }_i}},i=1,\dots ,10.$$

(92)

These expressions incorporate both the necessary optimality conditions and the box constraints on the controls.

5. Numerical Technique

To numerically approximate the fractional-order state system and associated adjoint equations arising from the application of Pontryagin’s maximum principle, we adopt an iterative forward-backward sweep method. The system of Caputo fractional differential equations governing the dynamics of the control model is expressed in compact vector form as follows:

$${\displaystyle \frac{1}{{\sigma }^{1-\alpha }}}{}^{C}D_t^{\alpha }\mathcal{X}\left(t\right)=\mathcal{F}\left(\mathcal{X}\right),\mathcal{X}\left(0\right)={\mathcal{X}}_0,t\in [0,T],$$

(93)

where $\mathcal{X}\left(t\right)$ and $\mathcal{F}\left(\mathcal{X}\right)$ are defined, respectively, as the state vector and the right-hand side of the system:

$$\mathcal{X}\left(t\right)=\left(S_h,E_h,I_h,I_a,R_h,S_v,E_v,I_v,I_{v_r},A\right),$$

(94)

$$\mathcal{F}\left(\mathcal{X}\right)=\left(f_1\left(\mathcal{X}\right),f_2\left(\mathcal{X}\right),f_3\left(\mathcal{X}\right),f_4\left(\mathcal{X}\right),f_5\left(\mathcal{X}\right),f_6\left(\mathcal{X}\right),f_7\left(\mathcal{X}\right),f_8\left(\mathcal{X}\right),f_9\left(\mathcal{X}\right),f_{10}\left(\mathcal{X}\right)\right),$$

(95)

and $\mathcal{F}\left(\mathcal{X}\right)$ encapsulates the right-hand sides of the fractional model defined in Equation (9). Using the Caputo definition of fractional derivative, Equation (93) can be rewritten in integral form:

$$\mathcal{X}\left(t\right)={\sigma }^{1-\alpha }\left({\mathcal{X}}_0+{\displaystyle \frac{1}{\mathsf{\Gamma }\left(\alpha \right)}}{\int }_0^t{(t-\tau )}^{\alpha -1}\mathcal{F}\left(\mathcal{X}\left(\tau \right)\right)d\tau \right).$$

(96)

We discretize the interval $[0,T]$ using a uniform step size $h=T/n$ and define time nodes ${\tau }_j=jh$, for $j=0,1,\dots ,n$. Applying the trapezoidal quadrature rule to approximate the integral in (96), we obtain the numerical scheme:

$${\mathcal{X}}_j={\sigma }^{1-\alpha }\left({\mathcal{X}}_0+{\displaystyle \frac{h}{\mathsf{\Gamma }\left(\alpha \right)}}\sum _{k=0}^j{({\tau }_j-{\tau }_k)}^{\alpha -1}{\omega }_{jk}\mathcal{F}\left({\mathcal{X}}_k\right)+{\mathcal{E}}_j\right),$$

(97)

where ${\omega }_{jk}={\displaystyle \frac{1}{2}}$ if $k=0$ or $k=j$, and ${\omega }_{jk}=1$ otherwise, and ${\mathcal{E}}_j$ is the local truncation error satisfying:

$$|{\mathcal{E}}_j|\le C{h}^2,\mathrm{with}C={\displaystyle \frac{T}{12}}\underset{\tau }{max}\left|{\mathcal{X}}^{\prime \prime }\left(\tau \right)\right|.$$

(98)

The system (97) is solved recursively using Newton’s method due to its nonlinear dependence on ${\mathcal{X}}_k$. To solve the optimal control problem, we apply the forward-backward sweep algorithm. First, the state equations in (80) are integrated forward in time using the above numerical method and initial conditions. Then, the adjoint equations in (88)—derived from the Hamiltonian—are solved backward in time using the transversality conditions. At each iteration, the control variables $u_1\left(t\right)$ to $u_4\left(t\right)$ are updated by minimizing the Hamiltonian function, according to the characterizations provided by the necessary optimality conditions in (91). The overall procedure is summarized below:

Algorithm 1 (Forward–backward sweep method for solving the fractional optimal control problem)

1: Input: Initial conditions ${\mathcal{X}}_0$, terminal time T, time step h, fractional order $\alpha$, initial guesses for controls $u_1^{\left(0\right)},u_2^{\left(0\right)},u_3^{\left(0\right)},u_4^{\left(0\right)}$. 2: Initialize: Set iteration index $m\leftarrow 0$. 3: repeat 4:     Step 1 (Forward Sweep): 5:     Solve the state system (80) over $[0,T]$ using the fractional trapezoidal rule and current controls $u_i^{\left(m\right)}\left(t\right)$ to obtain ${\mathcal{X}}^{\left(m\right)}\left(t\right)$. 6:     Step 2 (Backward Sweep): 7:     Solve the adjoint system (88) backward in time using the same numerical method, with transversality conditions at $t=T$, to obtain ${\lambda }^{\left(m\right)}\left(t\right)$. 8:     Step 3 (Update Controls): 9:     for all time nodes $t_j$ do 10:         Compute updated control $u_i^{(m+1)}\left(t_j\right)$ using the characterization from Pontryagin’s maximum principle (Equation (91)). 11:         Project updated $u_i^{(m+1)}\left(t_j\right)$ onto admissible control set $[0,1]$ and apply budget constraints if needed. 12:     end for 13:     Step 4 (Check Convergence): 14:     Compute relative errors between $({\mathcal{X}}^{(m+1)},{\lambda }^{(m+1)},u_i^{(m+1)})$ and their previous iterates. 15: until convergence criterion is met (e.g., $max\parallel u_i^{(m+1)}-u_i^{\left(m\right)}\parallel <ϵ$). 16: Output: Optimal state trajectory $\mathcal{X}\left(t\right)$, adjoint variables $\lambda \left(t\right)$, and control functions $u_i\left(t\right)$ for $i=1,\dots ,4$.

Figure 2. Flowchart of the forward–backward sweep algorithm.

Figure 2. Flowchart of the forward–backward sweep algorithm.

5.1. Convergence Analysis

The primary goal of the convergence analysis in our study is to demonstrate that the numerical method used to approximate the solution of the fractional-order optimal control system converges to the true solution as the step size $h\to 0$. Recall the integral approximation error in the fractional trapezoidal method, given by:

$$\underset{h\to 0}{lim}\parallel {\mathcal{E}}_j\parallel =\underset{h\to 0}{lim}\parallel {\mathcal{X}}_j-\mathcal{X}\left({\tau }_j\right)\parallel =0,$$

(99)

where ${\mathcal{X}}_j$ is the numerical approximation of the state vector $\mathcal{X}\left({\tau }_j\right)$ at node ${\tau }_j$. Subtracting the exact integral representation from its numerical approximation yields the error:

$$\parallel {\mathcal{E}}_j\parallel \le ∥h\sum _{k=0}^j{\omega }_{jk}{\mathcal{K}}_{jk}\left[\mathcal{F}\left({\mathcal{X}}_k\right)-\mathcal{F}\left(\mathcal{X}\left({\tau }_k\right)\right)\right]+R_j\left(\mathcal{G}\right)∥,$$

(100)

where ${\mathcal{K}}_{jk}={({\tau }_j-{\tau }_k)}^{\alpha -1}$, ${\omega }_{jk}$ are trapezoidal weights, and $R_j\left(\mathcal{G}\right)$ is the local truncation error of the quadrature rule. Assuming that $\mathcal{F}$ satisfies a Lipschitz condition with constant L, and letting $W={max}_{j,k}\left|{\omega }_{jk}\right|$, $M={max}_{j,k}\left|{\mathcal{K}}_{jk}\right|$, we have:

$$\parallel {\mathcal{E}}_j\parallel \le LhWM\sum _{k=0}^j\parallel {\mathcal{E}}_k\parallel +\parallel R_j\left(\mathcal{G}\right)\parallel .$$

(101)

For sufficiently small h, the following recursive bound applies:

$$\parallel {\mathcal{E}}_j\parallel \le {\displaystyle \frac{\parallel R_j\left(\mathcal{G}\right)\parallel }{1-LhWM}},$$

(102)

which leads to:

$$\parallel {\mathcal{E}}_j\parallel \le {\displaystyle \frac{\parallel R_j\left(\mathcal{G}\right)\parallel }{1-LhWM}}exp\left({\displaystyle \frac{LWM}{1-LhWM}}jh\right).$$

(103)

Therefore, as $h\to 0$, we conclude that $\parallel {\mathcal{E}}_j\parallel \to 0$, proving the convergence of the numerical scheme for the fractional-order epidemic model with optimal control.

6. Numerical Results

In the numerical simulations, we applied the method described in Section 5, implemented using the Python software. Furthermore, we consider various simulation scenarios, which we will discuss hereinafter.

6.1. Simulation Scenarios

To assess the effectiveness of various intervention strategies in mitigating the spread of ACL, we define a set of control-based simulation scenarios. Each scenario corresponds to a distinct combination of the four available control functions in our model: use of bed nets and repellents ($u_1$), treatment of symptomatic individuals ($u_2$), insecticide spraying ($u_3$), and awareness campaigns ($u_4$). These controls are applied subject to resource constraints, and their configurations are chosen to reflect both practical public health policies and theoretical significance. The following scenarios are considered:

• Scenario S₁ (Vector Avoidance Only): Figure 3 illustrates the temporal evolution of all ten state variables under Scenario $S_1$, where only the personal protection control $u_1\left(t\right)$ (e.g., use of bed nets and repellents) is active, and compares it to the uncontrolled baseline. The fractional order is fixed at $\alpha =0.95$ to account for memory effects inherent in both the biological and behavioral dynamics of the system. The implementation of $u_1\left(t\right)$ significantly reduces the force of infection between humans and vectors by lowering effective contact rates. As a result, a modest yet consistent decline is observed in the number of exposed ($E_h$) and symptomatic infected humans ($I_h$) throughout the simulation period, relative to the no-control scenario. The impact on asymptomatic infections ($I_a$) is less pronounced, given their indirect dependence on exposure rather than contact intensity. For the vector population, a slight reduction in exposed ($E_v$) and infected vectors ($I_v$, $I_{vr}$) is observed due to the reduced likelihood of acquiring infection from protected human hosts. However, since vector dynamics are not directly targeted in this scenario, the changes are relatively subdued. Importantly, the awareness function $A\left(t\right)$ shows no deviation between the two cases, as no awareness-related control is employed.
  Overall, Scenario $S_1$ demonstrates that personal protection can contribute meaningfully to reducing human infections, especially during the early and middle phases of the outbreak. However, its isolated application yields limited control over the broader transmission cycle, suggesting that more comprehensive strategies are required to achieve substantial disease suppression.
• Scenario S₂ (Treatment Only): Figure 4 presents the model dynamics under Scenario $S_2$, where treatment of symptomatic individuals via control $u_2\left(t\right)$ is the sole intervention. Compared to the baseline, a marked reduction in the symptomatic infected population ($I_h$) is achieved, particularly evident during the early and mid stages of the epidemic. This decline results from a shortened infectious period, which also indirectly reduces secondary infections.
  Consequently, a delayed but observable decline in the number of exposed humans ($E_h$) emerges due to the weakened forward transmission. The asymptomatic class ($I_a$), unaffected directly by treatment, exhibits minor differences from the baseline. For the vector populations, a mild reduction in exposed and infected compartments ($E_v$, $I_v$, $I_{vr}$) occurs, driven by the lowered human-to-vector transmission intensity. No change is observed in the awareness level ($A\left(t\right)$), as no awareness-related intervention is employed. Overall, Scenario $S_2$ proves effective in rapidly reducing symptomatic infections and partially suppressing the transmission cycle, but it does not fully disrupt vector-driven dynamics without complementary controls
• Scenario S₃ (Vector Spraying Only): Figure 5 depicts the evolution of state variables under Scenario $S_3$, where insecticide spraying ($u_3\left(t\right)$) is the only active control. The most pronounced impact is observed in the vector-related compartments, with substantial reductions in susceptible ($S_v$), exposed ($E_v$), and both infected classes ($I_v$, $I_{vr}$), indicating the direct mortality effect of the intervention on the vector population. This reduction in vector density effectively lowers the force of infection toward humans, leading to a delayed but meaningful decline in exposed ($E_h$) and symptomatic infected humans ($I_h$). The asymptomatic class ($I_a$), being indirectly affected, also exhibits a slight downward shift compared to the no-control case.
  Since treatment and awareness controls are inactive, the recovery and awareness levels ($R_h$, $A\left(t\right)$) follow trajectories similar to the baseline. Overall, vector spraying in isolation demonstrates notable efficiency in disrupting vector dynamics and indirectly suppressing human transmission, though it falls short of completely halting disease propagation in the absence of host-directed interventions.
• Scenario S₄ (Awareness Campaign Only): Figure 6 shows the system trajectories under Scenario $S_4$, where only public awareness efforts ($u_4\left(t\right)$) are implemented. As expected, the awareness level $A\left(t\right)$ exhibits a sustained and significant increase over time, reflecting the direct effect of the control. This rise in awareness reduces the effective contact rate between humans and vectors, which in turn mildly suppresses transmission. The impact on epidemiological compartments is modest but consistent: small reductions are observed in the exposed ($E_h$), symptomatic ($I_h$), and asymptomatic ($I_a$) human classes, along with slight decreases in infected vectors ($I_v$, $I_{vr}$). However, due to the indirect nature of this intervention and the absence of direct treatment or vector-killing measures, the declines are not pronounced.
  Overall, Scenario $S_4$ highlights the role of behavioral modification as a supportive control. While awareness alone does not significantly alter the disease trajectory, it contributes to reducing transmission intensity, especially when integrated with other strategies.
• Scenario S₅ (Combined Treatment and Personal Protection): Figure 7 illustrates the dynamics under Scenario $S_5$, which combines personal protection ($u_1\left(t\right)$) and treatment of symptomatic individuals ($u_2\left(t\right)$). This dual intervention targets both the infection source and the transmission pathway, leading to a more pronounced suppression of the epidemic compared to single-control strategies. A significant reduction in symptomatic infections ($I_h$) is observed, resulting from both shortened infectious periods and decreased exposure risk. The exposed class ($E_h$) also declines notably, benefiting from the synergistic effect of reduced contact and diminished secondary transmission. Asymptomatic infections ($I_a$), while not directly treated, follow a downward trend as a result of reduced upstream exposure. In the vector population, moderate declines are observed in exposed and infected compartments due to the reduced infectivity of human hosts. Since no vector-targeted control or awareness campaign is applied, the susceptible vector pool ($S_v$) and awareness level ($A\left(t\right)$) remain largely unchanged.
  Overall, Scenario $S_5$ demonstrates the enhanced effectiveness of integrated host-based strategies. The simultaneous use of repellents and treatment significantly weakens transmission feedback loops and offers a balanced approach to disease mitigation
• Scenario S₆ (Combined Repellent and Spraying): Figure 8 presents the evolution of the system under Scenario $S_6$, which integrates personal protection ($u_1\left(t\right)$) with insecticide spraying ($u_3\left(t\right)$). This combined strategy simultaneously reduces human-vector contact and vector abundance, resulting in a strong and rapid suppression of disease transmission. The vector-related compartments ($S_v$, $E_v$, $I_v$, $I_{vr}$) show substantial declines due to direct vector mortality induced by spraying. Additionally, the reduction in human exposure—achieved through repellents and bed nets—further limits the replenishment of the infected vector pool. On the human side, exposed ($E_h$), symptomatic ($I_h$), and asymptomatic ($I_a$) infections are all significantly reduced relative to the baseline. These improvements stem from the combined impact of fewer infectious bites and diminished forward transmission pressure. However, since no treatment or awareness interventions are applied, recovery dynamics and the awareness level $A\left(t\right)$ remain nearly identical to the uncontrolled case.
  In summary, Scenario $S_6$ effectively disrupts both the environmental and behavioral transmission routes. It proves particularly valuable in settings where vector control and personal-level measures can be deployed in tandem, even in the absence of clinical or educational interventions.
• Scenario S₇ (Combined Treatment and Awareness): Figure 9 depicts the model behavior under Scenario $S_7$, which combines treatment of symptomatic individuals ($u_2\left(t\right)$) with public awareness campaigns ($u_4\left(t\right)$). This hybrid approach merges direct clinical intervention with indirect behavioral modification, resulting in meaningful reductions across several compartments. Symptomatic infections ($I_h$) decline sharply due to active treatment, which shortens the infectious period and interrupts human-to-vector transmission. At the same time, the awareness function $A\left(t\right)$ increases steadily, decreasing exposure risk by reducing contact rates. These complementary effects yield a visible drop in the number of exposed humans ($E_h$), as well as a moderate decline in asymptomatic infections ($I_a$), which are indirectly influenced by reduced incidence. Although the vector compartments are not targeted directly, a slight reduction in infected vectors ($I_v$, $I_{vr}$) emerges due to decreased transmission from human hosts. As expected, the susceptible vector population ($S_v$) remains stable in the absence of vector-specific controls.
  Scenario $S_7$ illustrates the added value of combining behavioral and therapeutic strategies. The synergy between heightened public awareness and timely clinical response provides a two-pronged mechanism for curbing transmission, especially in human-dominated settings.
• Scenario S₈ (Full Control Strategy): Figure 10 shows the dynamics of the full intervention scenario $S_8$, in which all four controls—personal protection ($u_1$), treatment ($u_2$), insecticide spraying ($u_3$), and public awareness ($u_4$)—are simultaneously applied. This comprehensive strategy exerts maximal pressure on the disease system by targeting both human and vector dynamics as well as behavioral responses. The result is a pronounced and rapid decline in all key epidemiological compartments. Symptomatic ($I_h$) and asymptomatic ($I_a$) infections are significantly reduced due to the combined effects of reduced exposure, enhanced recovery, and increased awareness. The exposed human class ($E_h$) also decreases sharply, reflecting strong suppression of transmission at early stages. Vector-related states—including $S_v$, $E_v$, $I_v$, and $I_{vr}$—drop substantially as a consequence of direct vector mortality (via $u_3$) and indirect transmission reductions. Moreover, the awareness level $A\left(t\right)$ rises steadily, further reinforcing the behavioral shift needed to sustain control measures.
  Overall, Scenario $S_8$ delivers the most effective outcome among all considered strategies. By synchronizing environmental, clinical, and educational interventions, it disrupts multiple feedback loops in the transmission cycle and highlights the power of integrated public health approaches in managing ACL.

6.2. Effect of Fractional Order $\alpha$

Figure 11 demonstrates the influence of the fractional derivative order $\alpha$ on the temporal evolution of symptomatic infections $I_h\left(t\right)$ and the awareness level $A\left(t\right)$ under Scenario $S_4$, where only the awareness control $u_4\left(t\right)$ is active. As the value of $\alpha$ decreases from 1.0 to 0.7, representing a stronger memory effect, a significant suppression in both the peak and the long-term prevalence of $I_h\left(t\right)$ is observed. This behavior reflects the capacity of fractional-order dynamics to capture the prolonged impact of public health campaigns, where awareness persists and continues to influence behavior over time. Furthermore, the awareness variable $A\left(t\right)$ exhibits a faster and more sustained increase at lower $\alpha$, implying that individuals respond more effectively to information and retain it longer when memory effects are stronger. Overall, the results highlight that fractional models can provide a more realistic representation of behavioral and epidemiological processes, and that awareness campaigns may become considerably more effective in the presence of population memory and behavioral inertia. Figure 12 illustrates the effect of the fractional-order parameter $\alpha$ on the dynamics of symptomatic infections $I_h\left(t\right)$ and awareness level $A\left(t\right)$ under Scenario $S_8$, where all four control strategies ($u_1$ to $u_4$) are simultaneously active. The figure shows that as $\alpha$ decreases from 1.0 to 0.7, both the peak and the duration of symptomatic infections are further reduced, with the earliest suppression occurring at $\alpha =0.7$. This enhanced mitigation highlights the memory effect embedded in the fractional derivative: a lower $\alpha$ means that the system retains more information from the past, thereby making the control interventions—particularly awareness and behavioral changes—more effective over time. Similarly, the awareness function $A\left(t\right)$ increases more rapidly and reaches higher sustained levels at lower $\alpha$, reinforcing the idea that memory-dependent behavioral responses can amplify the impact of multi-faceted control programs. Overall, these findings emphasize that the synergy between fractional memory and combined interventions can significantly accelerate disease eradication.

Figure 13 shows how the objective functional $J\left(u\right)$ varies with the fractional-order derivative $\alpha$ across different intervention scenarios. As $\alpha$ increases toward the classical case $\alpha =1.0$, the values of $J\left(u\right)$ generally increase for all scenarios. This trend reflects the reduced impact of memory effects in integer-order models, where the system responds more locally in time and lacks long-term behavioral persistence. In contrast, lower values of $\alpha$ (e.g., $\alpha =0.7$) lead to significantly lower values of $J\left(u\right)$, indicating enhanced cost-effectiveness of control strategies under fractional dynamics. The memory effect embedded in the Caputo derivative allows past control efforts—especially in awareness and behavioral response—to exert influence over extended periods, thus reducing both infection burden and intervention costs more efficiently. Notably, scenarios with multiple active controls (e.g., Scenario S8) exhibit the greatest sensitivity to $\alpha$, suggesting that fractional-order models are particularly valuable in optimizing complex, multi-layered intervention programs.

6.3. Optimal Scenario Identification

To determine the most effective intervention strategy, we compare the total cost functional J and accumulated infection burden E across all control scenarios for the fractional-order $\alpha =0.95$, as shown in Table 5. The scenario involving no control results in the highest infection burden $E=1995.7209$ with a baseline cost of $J=490.9330$. The lowest value of J is achieved by the scenario $\{u_2,u_3\}$ with $J=481.1943$, which also substantially reduces the infection burden to $E=1919.9739$. Interestingly, adding awareness control $u_4$ to this combination ($\{u_1,u_2,u_4\}$ or “All Controls”) does not significantly improve either cost or effectiveness, suggesting a redundancy or saturation effect when all interventions are simultaneously applied. Scenarios that employ only awareness ($u_4$), or only bed-net usage ($u_1$), result in nearly identical infection levels to the no-control case, but with a higher total cost. This indicates that awareness campaigns alone—although valuable from a public health perspective—are insufficient to control the disease without complementary biomedical interventions such as treatment or vector control. These findings are quantitatively confirmed by the Incremental Cost-Effectiveness Ratio (ICER) values in Table 6. The scenario $\{u_2,u_3\}$ exhibits the lowest ICER of $-0.1286$, indicating a cost-saving effect relative to no control. In contrast, scenarios with only $u_1$ or $u_4$ yield highly inefficient outcomes, with ICER values exceeding 1900 and 36000 respectively. Even the “All Controls” scenario, despite achieving a slightly reduced infection level, shows a less favorable ICER compared to $\{u_2,u_3\}$ due to higher costs. In summary, the optimal control configuration in terms of cost-effectiveness and infection mitigation is the combined treatment and spraying scenario $\{u_2,u_3\}$. This highlights the biological importance of directly targeting disease progression and vector transmission rather than relying solely on awareness or preventive measures.

6.4. Scientific Interpretation, Policy and Public Health Implications, and Comparison with Classical Control Models

The present study proposes a novel fractional optimal control model for ACL, integrating key epidemiological features such as asymptomatic transmission, vector resistance, and behaviorally-driven awareness dynamics. By employing Caputo fractional derivatives, the model captures non-local temporal effects inherent to biological and behavioral processes, such as memory of past exposures or sustained public responses to interventions [10,23].

Scientifically, one of the major innovations is the explicit modeling of the asymptomatic human class, which reflects sub-clinical carriers who sustain transmission without detection or treatment. Previous studies have shown that ignoring asymptomatic individuals can result in underestimation of disease burden and intervention failure [19,20]. Our simulations confirm that targeting only symptomatic individuals leads to limited impact unless asymptomatic transmission is also suppressed. Another key contribution is the inclusion of a resistant vector compartment, modeling the gradual emergence of insecticide resistance—a growing concern in leishmaniasis-endemic regions [4]. The model demonstrates that even a small fraction of resistant vectors can significantly erode the effectiveness of vector control measures over time, supporting the call for resistance monitoring and rotation of control tools as recommended by WHO [32].

The introduction of a dynamic awareness function allows for the simulation of behaviorally adaptive responses in the human population. This component, which evolves based on campaign effort and decays in the absence of reinforcement, aligns with current findings in behavioral epidemiology that emphasize the time-dependent and nonlinear nature of public health compliance [14]. Under fractional-order dynamics, awareness campaigns exhibit greater long-term persistence and efficacy due to memory effects, highlighting the relevance of fractional modeling for socio-epidemiological systems [11]. From a policy and public health perspective, the model provides an actionable framework for resource-aware control planning. Through optimal control theory, it identifies dynamic strategies for allocating limited resources across treatment, protection, spraying, and education. Numerical results show that combinations such as treatment + vector spraying outperform single strategies in terms of cost-effectiveness—quantified via the objective functional and ICER values—especially under high transmission scenarios. These results are consistent with prior work on cost-effective leishmaniasis control in resource-constrained settings [24].

When compared to classical integer-order control models, the proposed fractional model offers notable advantages. Classical models assume memory-less dynamics and typically fail to capture delayed or accumulated effects of interventions. As demonstrated in previous research [21,22], fractional-order systems can reflect biologically realistic features such as incubation delays, treatment response variability, and behavioral inertia. Moreover, the smoother control profiles generated by the fractional approach (as opposed to abrupt switching seen in classical settings) are more compatible with real-world health program logistics and behavioral adaptation lags. In summary, this study extends the theoretical framework of disease modeling by incorporating biologically meaningful memory mechanisms, resistance evolution, and awareness-driven behavior into a single, resource-constrained fractional optimal control model. It bridges the gap between rigorous mathematical modeling and applied public health policy, and offers valuable insight into designing robust and sustainable intervention strategies for ACL and similar vector-borne diseases.

7. Conclusion and Future Work

This study presents a novel fractional-order optimal control model for ACL, capturing the complexities of disease dynamics through the integration of asymptomatic carriers, insecticide-resistant vectors, and behaviorally adaptive awareness. The use of Caputo fractional derivatives enables the model to reflect memory-dependent processes in both biological and behavioral domains, such as incubation delays, persistent public health responses, and the long-term evolution of resistance.

Key findings from the analysis highlight the importance of incorporating memory effects into epidemiological models. Our results show that intervention strategies such as treatment and vector spraying are significantly more effective under fractional-order dynamics ($\alpha <1$), as they benefit from the sustained effects of past actions. Additionally, modeling public awareness as a dynamic state variable — rather than a static parameter — allows the model to reflect intervention fatigue and the temporal evolution of behavior change. This is in line with recent behavioral modeling frameworks which emphasize feedback mechanisms in public response [15,16].

The model also illustrates that the presence of asymptomatic individuals can perpetuate transmission even in the presence of symptomatic treatment programs, underscoring the need for broader surveillance approaches [1]. Moreover, simulations involving insecticide-resistant vectors suggest that traditional spraying interventions may lose efficacy over time unless resistance management strategies are implemented. The biological relevance of these extensions demonstrates the model’s ability to bridge mathematical rigor with epidemiological realism.

From a public health perspective, the optimal control framework provides practical insights into how limited resources can be allocated among multiple interventions. In particular, combinations of interventions — especially treatment and vector control — are shown to yield high cost-effectiveness, while awareness-only strategies, although slower, produce more sustainable effects under memory-rich dynamics. These findings complement recent theoretical and empirical research emphasizing the role of long-term planning in epidemic control [12].Despite its strengths, the model has limitations. It assumes homogeneous population mixing and does not incorporate spatial heterogeneity or network-based interactions, which have been shown to significantly affect disease spread and control outcomes [33]. Additionally, some parameters — particularly those related to behavioral response and resistance evolution — were estimated or assumed in the absence of high-quality empirical data.

Future research may address these limitations through several directions. First, integrating stochastic effects and uncertainty analysis would enhance model robustness, especially under real-world unpredictability [34]. Second, applying Bayesian inference or machine learning methods to estimate behavioral and biological parameters from field data would improve the model’s predictive accuracy [13]. Third, incorporating spatial structure and meta-population dynamics — including human mobility and network connectivity — could reveal more nuanced transmission patterns and inform targeted intervention strategies [35]. Finally, expanding the model to account for co-infections (e.g., HIV-Leishmania) or immune memory effects may increase its relevance for broader public health contexts [36]. In summary, this work introduces a biologically meaningful and mathematically robust framework for modeling ACL, offering both theoretical innovation and practical policy implications. The use of fractional calculus enriches our understanding of memory-driven dynamics in disease transmission and control, making this approach valuable not only for ACL but also for a wide range of persistent vector-borne diseases.