Fractional-order analysis of a fear-induced ecoepidemiological predator–prey model with optimal control and bifurcation dynamics

Ecoepidemiological models provide a fundamental mathematical framework for studying the interplay between ecological interactions and infectious disease transmission. Since the seminal works of Anderson and May^1,2, the incorporation of disease dynamics into predator–prey systems has revealed a wide range of complex behaviors, including threshold phenomena, population persistence, oscillatory outbreaks, and extinction scenarios. Such models are particularly relevant in natural ecosystems where disease alters prey susceptibility to predation or where predators preferentially consume infected individuals^3,4,5.

Beyond direct predation, empirical evidence has demonstrated that predators exert significant nonconsumptive effects on prey populations. The mere presence of predators can induce fear, stress, and behavioral changes in prey, leading to reduced reproduction and altered foraging strategies^6,7,8. These fear-induced effects have been shown to strongly regulate population dynamics, even when actual predation rates are low. Motivated by these observations, several mathematical models have incorporated fear effects into predator–prey interactions, revealing that fear can stabilize equilibria, shift bifurcation thresholds, and fundamentally change system dynamics^{9,10,11,12,13}.

In ecoepidemiological contexts, fear effects play an even more subtle role by indirectly modifying disease transmission. Fear-induced suppression of prey reproduction reduces population density and contact rates, thereby influencing infection persistence and outbreak severity. Recent studies have highlighted that the interaction between fear, disease, and predation can generate rich dynamical behaviors, including Hopf bifurcations and complex oscillations^14,15,16,17. Nevertheless, most existing studies are formulated using classical integer-order differential equations, which may not adequately capture the memory-dependent nature of ecological and epidemiological processes.

Fractional calculus has emerged as a powerful modeling tool for incorporating memory and hereditary effects into dynamical systems. Fractional-order models, particularly those based on Caputo derivatives, have been successfully applied in epidemiology, predator–prey dynamics, and ecoepidemiological systems^18,19,20. One of the key advantages of fractional models is their ability to suppress oscillations, enlarge stability regions, and modify bifurcation structures. The Matignon stability criterion provides a natural extension of classical linear stability analysis to fractional-order systems, revealing that decreasing the fractional-order has a stabilizing effect on equilibria¹⁸. Despite these advantages, fractional-order formulations of fear-induced ecoepidemiological models remain relatively unexplored.

In parallel, optimal control theory has become an essential tool for designing cost-effective intervention strategies in epidemiological and ecological models. By balancing disease reduction against implementation costs, optimal control approaches provide valuable insights into prevention and treatment policies^21,22. Fractional optimal control problems have attracted increasing attention due to their ability to incorporate memory effects in both state dynamics and control actions^{23,24,25,26,27,28}. Recent developments have established rigorous theoretical foundations and efficient numerical methods for fractional optimal control, demonstrating significant reductions in total control costs compared to classical integer-order approaches^29,30,31,32. Recent studies have demonstrated the effectiveness of fractional optimal control in modeling complex biological and physical systems. In particular, fractional-order frameworks have been successfully applied to cancer dynamics, viral transmission, and tumor–immune interactions, providing improved accuracy and flexibility in capturing memory effects and nonlinear behaviors^33,34,35. In addition, fractional optimal control has been extended to quantum systems governed by time-fractional Schrödinger equations, highlighting its broad applicability beyond biological models³⁶.

Motivated by the above considerations, this paper develops and analyzes a fractional-order fear-induced ecoepidemiological predator–prey model with disease transmission in the prey population.

The model is formulated using Caputo fractional derivatives of order (0<nu le 1). This allows the inclusion of memory effects in the population dynamics. Fundamental analytical properties of the model, including positivity, boundedness, and existence and uniqueness of solutions, are rigorously established. Equilibrium points are derived, and their local stability is investigated using the Matignon stability criterion.

Furthermore, the influence of the fear factor is examined in detail through numerical simulations and bifurcation analysis with respect to key biological parameters, namely the conversion efficiency of predators, the hunting efficiency toward infected prey, and the infection transmission rate. An optimal control problem is then formulated to minimize the infected prey population and the cost of intervention. Necessary optimality conditions are derived using the fractional Pontryagin Maximum Principle, and numerical simulations demonstrate the effectiveness and cost-efficiency of the proposed control strategies. The results highlight the combined regulatory roles of fear effects and fractional memory in stabilizing ecoepidemiological systems and improving disease management outcomes.

Aim of the study.

The aim of this study is to develop and analyze a comprehensive fractional-order ecoepidemiological predator–prey model that simultaneously incorporates predator-induced fear, disease transmission in the prey population, and memory effects. The study seeks to quantify how fear and fractional dynamics interact to shape stability, bifurcation structure, and long-term population behavior, and to design optimal intervention strategies that efficiently suppress infection while minimizing implementation costs. A schematic illustration of the proposed model and its main mechanisms is presented in Fig. 1.

Innovation and applications

Unlike existing ecoepidemiological models that are predominantly formulated in an integer-order framework, this work introduces a fractional-order formulation in which predator-induced fear and memory effects coexist within a single unified model. The inclusion of Caputo fractional derivatives provides a realistic representation of hereditary and long-term biological processes, leading to qualitatively different dynamics that cannot be captured by classical models. In particular, the present study demonstrates that fractional memory can fundamentally alter bifurcation thresholds, suppress complex oscillations induced by fear and predation, and significantly reduce the cost of disease control. To the best of our knowledge, this is one of the first studies to integrate fear effects, fractional dynamics, bifurcation analysis, and optimal control within an ecoepidemiological predator–prey framework. Moreover, the proposed framework provides a bridge between theoretical fractional modeling and practical applications in ecological management and disease control.

Main contributions

The principal contributions of this paper are summarized as follows:

C1: A novel fractional-order ecoepidemiological predator–prey model with infected prey and predator-induced fear is formulated using Caputo fractional derivatives.

C2: Rigorous analytical results are established for the fractional system, including positivity, boundedness, and global existence and uniqueness of solutions.

C3: Equilibrium points are characterized and their local stability is analyzed using the Matignon stability criterion, highlighting the stabilizing role of fractional memory.

C4: A comprehensive numerical bifurcation analysis is performed with respect to the fear factor, conversion efficiency, hunting efficiency, and infection transmission rate, revealing transitions between stable equilibria, periodic oscillations, and complex dynamics.

C5: A fractional optimal control problem is formulated to minimize infection prevalence and intervention costs, and necessary optimality conditions are derived via the fractional Pontryagin Maximum Principle.

C6: Numerical simulations demonstrate that decreasing the fractional-order significantly reduces oscillatory behavior, infection burden, and total control costs, providing strong evidence of the practical benefits of fractional modeling in ecoepidemiological systems.

Related applications

The proposed model has several important applications in ecological and epidemiological contexts:

• Wildlife disease management: The model can be used to study disease transmission in prey populations and to design effective strategies for controlling outbreaks in natural ecosystems.

• Predator–prey interaction analysis: The inclusion of fear effects allows the model to capture non-consumptive predator impacts, which are increasingly recognized as important in ecological systems.

• Epidemiological control strategies: The optimal control framework provides insights into cost-effective intervention policies, such as prevention and treatment measures.

• Complex dynamical systems analysis: The model can be applied to investigate bifurcation phenomena and transitions between stable and oscillatory regimes in biological systems.

• Fractional modeling of biological processes: The incorporation of memory effects makes the model suitable for systems where past states influence current dynamics, such as disease progression and behavioral responses.

Organization of the paper

The remainder of this paper is organized as follows. Section “Introduction” introduces the ecological and epidemiological background of the study and outlines the motivation, aims, and main contributions. Section “Preliminaries” presents the necessary preliminaries from fractional calculus, including Caputo derivatives, Mittag–Leffler functions, and stability tools required throughout the analysis. In Section “Model formulation (baseline system)”, the classical integer-order fear-induced ecoepidemiological predator–prey model is formulated, and the biological interpretation of variables and parameters is discussed. Section “Fractional-order model formulation (Caputo derivative)” extends the baseline model to a fractional-order framework using Caputo derivatives, incorporating memory effects into the population dynamics. Section “Well-posedness of the fractional-order model” is devoted to the well-posedness analysis of the fractional-order system, where positivity, uniform boundedness, and existence and uniqueness of solutions are rigorously established. In Section “Fractional optimal control problem and Pontryagin Maximum Principle”, a fractional optimal control problem is formulated, and the existence of an optimal control pair is proved. Section “Stability analysis and fractional Lyapunov approach” derives the necessary optimality conditions using the fractional Pontryagin Maximum Principle, including the adjoint system, transversality conditions, and explicit characterizations of the optimal controls. Section “Stability analysis of the equilibrium points when susceptible prey experiences no predation” investigates the local stability of equilibrium points using the Matignon stability criterion and fractional Lyapunov methods, with special attention to scenarios where susceptible prey experiences no predation. Section “Fractional forward–backward sweep method (numerical scheme)” presents the numerical scheme based on a fractional forward–backward sweep algorithm used to solve the optimality system. Section “Influence of fear factor k₁ and bifurcation analysis with respect to θ, β, and λ” provides how the fear factor (k_1) affects the qualitative dynamics of the fractional ecoepidemiological system (4.1), and we perform bifurcation analysis with respect to the conversion efficiency (theta), the predator hunting efficiency toward infected prey (beta), and the infection rate (lambda). Section “Numerical simulations, bifurcation, and fractional-order effects”, in this section, numerical simulations are carried out to illustrate the analytical results and to investigate the complex dynamics of system (4.1). Finally, Section “Conclusion and future work” concludes the paper with a summary of the main findings and outlines possible directions for future research.

In this section, we briefly recall essential concepts and results from fractional calculus that will be used throughout the paper. For more details, the reader is referred to^18,19.

Fractional integrals and derivatives

Definition 2.1

(Riemann–Liouville fractional integral) Let (fin L^1([0,T])) and (nu>0). The Riemann–Liouville fractional integral of order (nu) is defined by

where (Gamma (cdot )) denotes the Gamma function.

Definition 2.2

(Caputo fractional derivative) Let (n-1<nu <n) with (nin mathbb {N}) and (fin C^n([0,T])). The Caputo fractional derivative of order (nu) is defined as

In particular, for (0<nu <1),

The Caputo derivative is especially suitable for modeling real-world processes, since it allows the use of classical initial conditions expressed in terms of integer-order derivatives.

Equivalent integral formulation

Fractional differential equations with Caputo derivatives are equivalent to Volterra-type integral equations.

Lemma 2.3

(Integral representation) Let (0<nu <1). The fractional differential equation

is equivalent to the integral equation

This formulation plays a crucial role in establishing existence, uniqueness, and positivity of solutions.

Mittag–Leffler functions

Definition 2.4

(Mittag–Leffler function) The one-parameter Mittag–Leffler function is defined by

The Mittag–Leffler function generalizes the exponential function and naturally appears in solutions of linear fractional differential equations.

Fractional Grönwall inequality

Lemma 2.5

(Fractional Grönwall inequality) Let u(t) be a nonnegative locally integrable function satisfying

for constants (c,dge 0) and (0<nu <1). Then

This inequality is widely used to prove boundedness and continuous dependence of solutions to fractional systems.

Stability of fractional-order systems

Consider the linear fractional system

where (Ain mathbb {R}^{ntimes n}).

Theorem 2.6

(Matignon stability criterion¹⁸) The equilibrium (x=0) of system (7.1) is locally asymptotically stable if and only if

where (lambda _j) are the eigenvalues of the matrix A.

This result shows that decreasing the fractional-order (nu) enlarges the stability region, explaining the stabilizing effect of fractional dynamics observed in numerical simulations.

Admissible controls

Let (T>0) be fixed. The set of admissible controls is defined as

This set is nonempty, convex, and closed, which is essential for the existence of an optimal control.

State variables and parameters

Let

denote the densities of susceptible prey, infected prey, and predator at time t, respectively.

Baseline fear-induced ecoepidemiological model

Under the assumptions in³⁷ (linear incidence; infected prey do not recover; Holling type II predation; predator switching; only susceptible prey experiences fear), the model is (see³⁷)

with (s(0)=s_0>0), (i(0)=i_0>0), (p(0)=p_0>0) and (a s_0+b i_0>0).

Basic reproduction number

Following³⁷, we define

where s denotes the (relevant) susceptible prey level used for threshold analysis.

State variables

Let

denote the densities of susceptible prey, infected prey, and predator at time t, respectively.

Fractional fear-induced ecoepidemiological model

Replacing the classical time derivatives in (3.1) by Caputo derivatives of order (nu in (0,1)), the fractional-order model is given by

subject to the fractional initial conditions

with (a s_0+b i_0>0), where the biological descriptions of the variables and the parameters will be described in Table 1.

A schematic representation of the main biological interactions and the fractional-memory structure of the model is shown in Fig. 2.

Remark 4.1

If different memory effects are desired for each population, one may consider distinct fractional-orders (0<nu _1,nu _2,nu _3<1):

For simplicity, we use a common order (nu) in (4.1).

Basic reproduction number

Following the reference, we define the basic reproduction number (threshold quantity) as

where s denotes the relevant susceptible prey level used for the threshold analysis.

Throughout this section, we consider the fractional-order system (4.1) with Caputo derivatives of order (0<nu <1).

Positivity of solutions

Theorem 5.1

(Positivity) If (s_0>0), (i_0>0), and (p_0>0), then the solution of the fractional-order system (4.1) satisfies

Proof

By the Caputo–Volterra equivalence^18,19, the solution satisfies

where (x(t)=(s(t),i(t),p(t))^top) and (x_0in mathbb {R}_+^3).

Step 1: Nonnegativity is invariant. We check that the vector field does not point outside (mathbb {R}_+^3) on the boundary:

because each (f_j) contains the factor of its state variable (s in (f_1), i in (f_2), and p in (f_3)) in every term. Therefore, once any component reaches 0, its Caputo derivative becomes 0 at that boundary state, so the component cannot cross into negative values. This establishes positive invariance of (mathbb {R}_+^3) for the fractional system (a standard argument in fractional predator–prey eco-models; see²⁰).

Step 2: Strict positivity for strictly positive initial data. Assume (s_0,i_0,p_0>0). If, for instance, s(t) attained 0 at some first time (t_1>0), then by continuity (s(tau )ge 0) for (tau in [0,t_1]) and (5.1) yields

The integral term is finite (since (f_1) is continuous on the invariant nonnegative region and ((t_1-tau )^{nu -1}) is integrable for (0<nu <1)), hence the right-hand side cannot equal 0 because (s_0>0). This contradiction shows (s(t)>0) for all (tin [0,T]). The same argument applies to i(t) and p(t).

Hence s(t), i(t), p(t) remain strictly positive on [0, T]. (square)

Uniform boundedness

Theorem 5.2

(Uniform boundedness) All solutions of the fractional-order system (4.1) with nonnegative initial data are uniformly bounded on [0, T] for any finite (T>0).

Proof

Let (s(t), i(t), p(t)) be a nonnegative solution on [0, T].

Step 1: Bounding s (t). From the first equation of (4.1) and using (p(t)ge 0), (i(t)ge 0),

Let y(t) solve the fractional logistic equation

By the fractional comparison principle^18,19 (applied in fractional predator–prey contexts in²⁰),

Moreover, (y(t)le max {s_0,k}) for all (tge 0), hence

Step 2: Aggregate function W (t). Define

By linearity of the Caputo derivative,

Substituting (4.1), the infection terms (pm lambda s i) cancel, and the predation terms cancel as well, giving

Using (frac{1}{1+k_1p}le 1) and (1-frac{s+i}{k}le 1), we obtain

By (5.2), (r s(t)le M:=rS_{max }). Also, since (i(t)le W(t)) and (p(t)le theta W(t)), we have

Therefore,

Step 3: Mittag–Leffler bound via fractional Grönwall. Consider the linear comparison equation

Its solution is

where (E_nu) is the Mittag–Leffler function^18,19. By the fractional comparison (or Grönwall-type) inequality^18,19, (5.3) yields (0le W(t)le Z(t)). Since (0<E_{nu }(-delta t^nu )le 1) for (delta>0), we conclude

Step 4: Componentwise boundedness. Because (0le s(t)le W(t)), (0le i(t)le W(t)), and (0le p(t)le theta W(t)), boundedness of W implies boundedness of s, i, p on [0, T].

Hence all solutions are uniformly bounded on [0, T]. (square)

Existence and uniqueness

Theorem 5.3

(Existence and uniqueness) For any initial condition ((s_0,i_0,p_0)in mathbb {R}_+^3) such that (a s_0+b i_0>0), the fractional-order system (4.1) admits a unique solution (s(t), i(t), p(t)) on [0, T] for any finite (T>0).

Proof

Let (x(t)=(s(t),i(t),p(t))^top) and write (4.1) as

where (F=(f_1,f_2,f_3)^top) is the nonlinear vector field defined by the right-hand side of (4.1). Define the domain

Since (a,b>0) and (a s_0+b i_0>0), we have (x_0in Omega).

Step 1: Volterra integral formulation. For Caputo derivatives, (5.4) is equivalent to the Volterra integral equation

see^18,19. (This reduction is also standard in fractional predator–prey modeling; cf²⁰..)

Step 2: Local Lipschitz property of F on bounded subsets of (Omega). Fix (R>0) and let

On the compact set (Omega _R), the quantity (D=a s+b i) satisfies (Dge delta _R>0) for some (delta _R). Hence all rational terms of the form (frac{s p}{a s+b i}) are (C^1) on (Omega _R), and so F is continuously differentiable on (Omega _R). Therefore, F is Lipschitz on (Omega _R): there exists (L_R>0) such that

Step 3: Contraction mapping. Let (X=C([0,T_0],mathbb {R}^3)) with (Vert xVert _infty =sup _{tin [0,T_0]}Vert x(t)Vert). Define the operator (mathcal {T}) on X by

Using (5.6), for x, y with values in (Omega _R) we obtain

Choose (T_0>0) such that (frac{L_R T_0^nu }{Gamma (nu +1)}<1); then (mathcal {T}) is a contraction, hence has a unique fixed point on ([0,T_0]). This fixed point solves (5.5), thus (4.1) has a unique local solution^18,19.

Step 4: Extension to [0, T].

By Theorems 5.1 and 5.2, solutions remain in a bounded subset of (Omega) for all (tin [0,T]). Thus the local solution can be continued stepwise uniquely up to any finite (T>0) (standard continuation argument for Caputo IVPs; see^18,19).

Hence, (4.1) admits a unique solution on [0, T]. (square)

Remark 5.4

The proof above follows the standard Caputo (leftrightarrow) Volterra reduction and contraction-mapping strategy widely used in fractional eco-dynamical models; see, for instance, the fractional predator–prey formulation and analysis framework in²⁰. If your thesis requires an explicit citation to a dedicated fractional-calculus monograph for the Banach fixed point framework, you may add one classical reference (e.g., Diethelm/Podlubny) to the bibliography; however, the present proof uses only items already included in your reference list.

In this section we extend the fractional-order model (4.1) to an optimal control setting. Let (0<nu <1) and denote the Caputo derivative by ({}^{C}D_t^nu).

Controlled fractional-order system

We introduce two bounded controls:

• (u_1(t)): prevention effort reducing effective transmission, (lambda mapsto lambda (1-u_1(t))),

• (u_2(t)): treatment/removal effort increasing removal of infected prey, (-u_2(t)i).

The admissible control set is

The controlled fractional system is

with initial conditions (s(0)=s_0>0), (i(0)=i_0>0), (p(0)=p_0>0) and (a s_0+b i_0>0).

Objective functional

We minimize infection burden and intervention costs:

where (A,B_1,B_2>0).

Problem 6.1

(Fractional optimal control problem) Find ((u_1^*,u_2^*)in mathcal {U}) such that

subject to (6.2).

Existence of an optimal control

Theorem 6.1

(Existence of an optimal pair) Problem 6.1 admits at least one optimal control pair ((u_1^*,u_2^*)in mathcal {U}).

Proof

We prove existence by the direct method of the calculus of variations adapted to fractional dynamical constraints.

Step 1: Nonemptiness, convexity, and compactness of the admissible set. The admissible control set

is nonempty (e.g., (u_1equiv 0), (u_2equiv 0)), convex, and bounded in (L^infty ([0,T])^2). Hence, by the Banach–Alaoglu theorem, (mathcal {U}) is weak-(*) compact in (L^infty ([0,T])^2).

Step 2: Existence of state trajectories for each admissible control. Fix any ((u_1,u_2)in mathcal {U}). The controlled fractional system (6.2) can be written as

with (x(t)=(s(t),i(t),p(t))^top). Since ((u_1,u_2)) are essentially bounded and F is continuous in t and locally Lipschitz in x on the biologically feasible region (because (a s+b i>0) is preserved), it follows from standard existence and uniqueness theory for Caputo fractional differential equations that there exists a unique state solution (x(cdot )in C([0,T],mathbb {R}_+^3)); see^18,19. Moreover, by Theorems 5.1 and 5.2, the solution remains nonnegative and bounded on [0, T].

Step 3: Existence of a minimizing sequence. Let

Since the integrand is nonnegative, (J(u_1,u_2)ge 0) for all admissible controls, hence (J^*) is finite. Therefore, there exists a minimizing sequence ({(u_1^n,u_2^n)}_{nge 1}subset mathcal {U}) such that

Step 4: Weak-(*) convergence of controls.

Because (mathcal {U}) is weak-(*) compact in (L^infty ([0,T])^2), there exists a subsequence (still denoted the same) and a limit control ((u_1^*,u_2^*)in mathcal {U}) such that

Step 5: Convergence of associated state trajectories. Let (x^n(t)=(s^n(t),i^n(t),p^n(t))^top) be the unique solution of (6.2) corresponding to ((u_1^n,u_2^n)). By uniform boundedness (Theorem 5.2) the family ({x^n}) is uniformly bounded in (C([0,T],mathbb {R}^3)). Using the Volterra integral formulation of the Caputo system (see^18,19),

and the boundedness of F on bounded sets, one obtains equicontinuity of ({x^n}) on [0, T]. Thus, by the Arzelà–Ascoli theorem, there exists a further subsequence (not relabeled) and (x^*in C([0,T],mathbb {R}^3)) such that

Passing to the limit in the integral equation (using dominated convergence and the weak-(*) convergence of (u_j^n) together with the uniform convergence of (x^n)) shows that (x^*) is the state trajectory corresponding to the control ((u_1^*,u_2^*)).

Step 6: Lower semicontinuity of the cost functional. The running cost

is convex in ((u_1,u_2)) and continuous in all variables. By weak lower semicontinuity of convex integral functionals and the convergences above,

Hence (J(u_1^*,u_2^*)=J^*), so ((u_1^*,u_2^*)) is an optimal control pair.

Step 7: Fractional optimal control existence references. The above argument is consistent with the existence framework widely used for fractional optimal control problems with Caputo dynamics; see, for example, Agrawal’s foundational formulations and existence arguments^23,24,25 and related Caputo fractional control studies^26,27,28,29. This completes the proof. (square)

Fractional Pontryagin Maximum Principle (PMP)

For fractional problems of Caputo type, the adjoint system typically involves the right-sided Riemann–Liouville (or Caputo) fractional derivative. We denote the right-sided Riemann–Liouville derivative of order (nu) by ({}^{RL}D_{t,T}^{nu }):

Hamiltonian

Let (psi _1(t),psi _2(t),psi _3(t)) be adjoint variables. Define the Hamiltonian

where ((f_1,f_2,f_3)) are the right-hand sides of (6.2).

Theorem 6.2

(Fractional PMP: necessary optimality conditions) Let ((u_1^*,u_2^*)in mathcal {U}) be an optimal control with corresponding state ((s^*,i^*,p^*)) solving (6.2). Then there exist adjoint functions (psi _1,psi _2,psi _3) such that:

1. (i)

   State system: ((s^*,i^*,p^*)) satisfies (6.2).

2. (ii)

   Adjoint system:

   $$begin{aligned} {}^{RL}D_{t,T}^{nu }psi _j(t)=frac{partial mathcal {H}}{partial x_j}big (t,s^*,i^*,p^*,u_1^*,u_2^*,psi (t)big ), quad (x_1,x_2,x_3)=(s,i,p), end{aligned}$$

   (6.6)

   where ({}^{RL}D_{t,T}^{nu }) is the right-sided Riemann–Liouville fractional derivative of order (nu).

3. (iii)

   Transversality:

   $$begin{aligned} psi _1(T)=psi _2(T)=psi _3(T)=0. end{aligned}$$

   (6.7)
4. (iv)

   Optimality (pointwise minimization): for a.e. (tin [0,T]),

   $$begin{aligned} mathcal {H}big (t,s^*,i^*,p^*,u_1^*,u_2^*,psi big ) =min _{(u_1,u_2)in [0,1]^2}mathcal {H}big (t,s^*,i^*,p^*,u_1,u_2,psi big ). end{aligned}$$

   (6.8)

Proof

The proof follows the classical variational approach but uses the fractional integration-by-parts formula that links the left-sided Caputo derivative in the state equations to a right-sided Riemann–Liouville derivative in the adjoint equations. This is standard in fractional optimal control formulations; see Agrawal^23,24,25 and related developments^{26,27,28,29,38}.

Step 1: Perturbation of controls and the induced state variation. Let ((u_1^*,u_2^*)in mathcal {U}) be optimal and let ((v_1,v_2)in L^infty ([0,T])^2) be an admissible variation direction. For sufficiently small (varepsilon), define the perturbed controls

and (if necessary) project them back into [0, 1] so that ((u_1^varepsilon ,u_2^varepsilon )in mathcal {U}). Let (x^varepsilon (t)=(s^varepsilon ,i^varepsilon ,p^varepsilon )^top) be the corresponding state solution of (6.2).

Define the state variation

which exists in (C([0,T])^3) under the differentiability assumptions on the dynamics (standard for Caputo systems with locally Lipschitz right-hand side; see^18,19).

Step 2: Linearized fractional variational system. Let f(t, x, u) denote the right-hand side of (6.2) (vector form). Linearizing around ((x^*,u^*)) yields the fractional variational equation

where (delta u(t)=(v_1(t),v_2(t))^top).

Step 3: First variation of the cost. The objective is

Hence the first variation is

where (L_x) and (L_u) are evaluated along ((x^*(t),u^*(t))). Because ((u_1^*,u_2^*)) is optimal, we have (delta Jge 0) for all admissible variations (delta u) consistent with the control constraints.

Step 4: Introduction of adjoints and fractional integration by parts. Introduce adjoint functions (psi (t)=(psi _1,psi _2,psi _3)^top) and consider the identity

which holds because (6.9) is satisfied.

Using the fractional integration-by-parts formula that couples the left-sided Caputo derivative with the right-sided Riemann–Liouville derivative (see^18,19 and its use in fractional optimal control^23,24), we obtain

where (mathcal {B}(t)) is the boundary term. Because (delta x(0)=0) and the terminal state is free, we enforce the transversality condition

so that the boundary contribution vanishes, yielding

Step 5: Derivation of the adjoint system. Add (6.11) to (6.10) and use (6.13). We obtain

Choose (psi) so that the coefficient of (delta x) vanishes for all (delta x), i.e.,

Recalling the Hamiltonian

we have (frac{partial mathcal {H}}{partial x}=L_x+f_x^top psi), hence the above equation is equivalent to

which is exactly (6.6) componentwise. The terminal condition (psi (T)=0) is the transversality condition (6.7).

Step 6: Optimality condition (Hamiltonian minimization). With the adjoint choice above, the first variation reduces to

For unconstrained controls, optimality implies (frac{partial mathcal {H}}{partial u}=0) a.e. For bounded controls (0le u_jle 1), this yields the variational inequality

which is equivalent to the pointwise minimization statement (6.8). (These steps follow the standard control-constrained fractional optimality derivations in^{23,24,25,26,27,29}.)

Thus, (i)–(iv) hold and the theorem is proved. (square)

Characterization of the optimal controls

Compute the partial derivatives:

Thus the unconstrained minimizers are

Projecting onto [0, 1] gives

where (operatorname {Proj}_{[0,1]}(z)=min {1,max {0,z}}).

Linearization and local stability criterion

Let (x(t)=(s(t),i(t),p(t))^top) and consider the fractional system

Let (x^*) be an equilibrium, (F(x^*)=0), and define (y=x-x^*). The linearized system is

where (J(x^*)) is the Jacobian of F evaluated at (x^*).

Theorem 7.1

(Matignon stability criterion) Consider the linear fractional system

where (Ain mathbb {R}^{ntimes n}) is a constant matrix. The equilibrium (yequiv 0) is (locally) asymptotically stable if and only if

Consequently, for the nonlinear system ({}^{C}D_t^nu x=F(x)), an equilibrium (x^*) is locally asymptotically stable if and only if

where (J(x^*)) is the Jacobian of F at (x^*).

Proof

Step 1: Reduction to Jordan form. Let (A=PJP^{-1}) be the Jordan decomposition of A over (mathbb {C}), where J is the Jordan canonical form. Setting (z=P^{-1}y), system (7.2) becomes

Thus stability of (yequiv 0) is equivalent to stability of (zequiv 0).

Step 2: Scalar fractional equation and Mittag–Leffler solution. Consider first the scalar Caputo equation

with (lambda in mathbb {C}). The unique solution is given by the Mittag–Leffler function (see^18,19)

Step 3: Asymptotic decay sector for (E_nu (lambda t^nu )).

A fundamental asymptotic property of the Mittag–Leffler function states that, for (0<nu <1),

where (arg (lambda )in (-pi ,pi ]). This result is classical in fractional differential equations; see^18,19. Therefore, the scalar solution (7.4) converges to 0 if and only if (|arg (lambda )|>frac{nu pi }{2}).

Step 4: Jordan blocks. If J contains a Jordan block corresponding to an eigenvalue (lambda), the solution contains terms involving (t^{m}E_{nu ,m+1}(lambda t^nu )) (generalized Mittag–Leffler functions). These terms also decay to 0 as (trightarrow infty) provided (|arg (lambda )|>frac{nu pi }{2}), while they do not decay if the condition fails. Hence, the equilibrium of (7.2) is asymptotically stable if and only if every eigenvalue of A lies outside the sector

Step 5: Application to nonlinear systems.

Let (x^*) be an equilibrium of ({}^{C}D_t^nu x=F(x)) and set (y=x-x^*). Then

Thus, local asymptotic stability of (x^*) is determined by the linearized system (7.2) with (A=J(x^*)), which yields the stated criterion.

This completes the proof. (square)

Remark 7.2

For (nu =1), this reduces to the classical condition (Re (lambda _j)<0). For (0<nu <1), the stability region is a sector in the complex plane, larger than the left half-plane.

Corollary 7.3

(Practical stability checks for the (3times 3) Jacobian) Let (x^*) be an equilibrium of the nonlinear fractional system ({}^{C}D_t^nu x=F(x)), (0<nu <1), and let (J^*:=J(x^*)in mathbb {R}^{3times 3}). Denote the eigenvalues of (J^*) by (lambda _1,lambda _2,lambda _3).

1. (a)

   If all eigenvalues are real, then (x^*) is locally asymptotically stable if and only if

   $$lambda _1<0,quad lambda _2<0,quad lambda _3<0.$$
2. (b)

   If one eigenvalue is real, (lambda _1in mathbb {R}), and the other two form a complex conjugate pair (lambda _{2,3}=rho pm iomega) with (omega ne 0), then (x^*) is locally asymptotically stable if and only if

   $$lambda _1<0 quad text {and}quad |arg (rho +iomega )|>frac{nu pi }{2}.$$

   Equivalently, since (|arg (rho +iomega )|=arctan !big (frac{|omega |}{|rho |}big )) for (rho <0), the complex-pair condition can be written as

   $$begin{aligned} rho <0 quad text {and}quad frac{|omega |}{|rho |}>tan !left( frac{nu pi }{2}right) . end{aligned}$$

   (7.6)
3. (c)

   Let the characteristic polynomial of (J^*) be

   $$begin{aligned} chi (lambda )=lambda ^3+c_1lambda ^2+c_2lambda +c_3, end{aligned}$$

   (7.7)

   where

   $$c_1=-operatorname {tr}(J^*),qquad c_2=text {sum of principal }2times 2text { minors of }J^*,qquad c_3=-det (J^*).$$

   If the eigenvalues are all real, then (c_1>0), (c_2>0), (c_3>0) is necessary and (together with (Delta ge 0)) sufficient for stability. If a complex conjugate pair exists, then (c_1>0), (c_2>0), (c_3>0) remain necessary, while the additional sector condition (7.6) must be checked for the complex pair.

Proof

(a) If (lambda _jin mathbb {R}), then (arg (lambda _j)=0) for (lambda _j>0) and (arg (lambda _j)=pi) for (lambda _j<0). Since (0<nu <1), we have (frac{nu pi }{2}in (0,frac{pi }{2})), hence

Thus Theorem 7.1 reduces to the classical condition (lambda _j<0) for all real eigenvalues.

(b) For a conjugate pair (rho pm iomega) with (omega ne 0), the stability region in Theorem 7.1 is the complement of the sector (|arg (lambda )|le frac{nu pi }{2}). When (rho <0), the argument satisfies

Thus (|arg (rho +iomega )|>frac{nu pi }{2}) is equivalent to

Since (0<nu <1), (pi -frac{nu pi }{2}in (frac{pi }{2},pi )) and the inequality above reduces to (7.6). The real eigenvalue must also satisfy (lambda _1<0) by part (a).

(c) The coefficient relations in (7.7) are standard: (c_1=-sum lambda _j), (c_2=sum _{i<j}lambda _ilambda _j), (c_3=-lambda _1lambda _2lambda _3). Hence positivity of (c_1,c_2,c_3) is necessary for all eigenvalues to have negative real parts when they are real, and remains a necessary constraint in the mixed real/complex case. The decisive fractional feature is the sector condition from Theorem 7.1, which must be checked directly for the complex pair (equivalently using (7.6)). The underlying decay sector for the Mittag–Leffler solution and its link to stability follows from fractional differential equation theory^18,19. (square)

Remark 7.4

(How decreasing (nu)enlarges the stability region) The stability requirement is (|arg (lambda )|>frac{nu pi }{2}). As (nu) decreases, the threshold angle (frac{nu pi }{2}) decreases, so the forbidden sector (|arg (lambda )|le frac{nu pi }{2}) becomes narrower. Therefore, for smaller (nu), more eigenvalues satisfy the stability condition, meaning that memory effects ((nu <1)) tend to stabilize equilibria compared with the integer-order case (nu =1)^18,19.

In this section we investigate the stability of equilibrium points under the special scenario in which the susceptible prey is not subject to predation. Mathematically, this corresponds to setting the susceptible-predation term to zero, i.e.

which is achieved by taking (alpha =0) in the model. Then the fractional-order system (4.1) reduces to

with (0<nu <1) and nonnegative initial conditions.

Equilibrium points

An equilibrium (E=(s^*,i^*,p^*)) satisfies the algebraic system obtained by setting the right-hand sides of (8.1) to zero.

(i) Trivial equilibrium

The extinction equilibrium is

(ii) Disease-free prey-only equilibrium

Setting (i^*=0) and (p^*=0) in the first equation yields

Thus the disease-free prey-only equilibrium is

(iii) Predator-free endemic equilibrium

Assume (p^*=0) and (i^*>0). Then the infected equation gives

Substituting into the susceptible equation (with (p^*=0)) gives

Hence

Therefore, an endemic predator-free equilibrium exists if

which is the usual epidemiological threshold^1,2. We denote this equilibrium by

(iv) Coexistence equilibrium

Assume (p^*>0) and (i^*>0). From the predator equation of (8.1):

Using this in the infected equation:

Substituting (i^*=frac{mu _2}{theta beta }) into the susceptible equation gives an implicit equation for (s^*):

Any positive solution (s^*>frac{mu _1}{lambda }) of (8.3) yields a positive coexistence equilibrium

The coexistence equilibrium therefore requires (lambda s_3^*>mu _1) (equivalently (p_3^*>0)).

Remark 8.1

Equation (8.3) captures the impact of fear (through (k_1)) on coexistence. Larger fear intensity increases the denominator (1+k_1p^*), reducing prey reproduction and potentially eliminating coexistence. This is consistent with fear-effect modeling observations in predator–prey systems^9,10.

Local stability via Matignon criterion

Let (F=(f_1,f_2,f_3)) be the right-hand side of (8.1). The Jacobian matrix J(s, i, p) is

Local asymptotic stability of an equilibrium E is determined by the eigenvalues of J(E) through the Matignon criterion (Theorem 7.1), whose proof is based on Mittag–Leffler asymptotics^18,19.

Stability of (E_0=(0,0,0))

At (E_0), the linearization is dominated by prey growth:

and other terms vanish. Thus (J(E_0)) has an eigenvalue (r>0), so (|arg (r)|=0le nu pi /2), and (E_0) is unstable for all (0<nu le 1).

Stability of (E_1=(k,0,0))

Compute the Jacobian at (E_1). Since (p^*=0) and (i^*=0),

For the infected component,

For the predator component,

Hence the eigenvalues of (J(E_1)) are

By Theorem 7.1, (E_1) is locally asymptotically stable if and only if all eigenvalues satisfy (|arg (lambda _j)|>nu pi /2). Since (lambda _1,lambda _3<0), their arguments are (pi) and automatically satisfy the criterion. Thus stability reduces to (lambda _2<0), i.e.

which matches the standard threshold interpretation^1,2.

Stability of (E_2=big (frac{mu _1}{lambda },i_2^*,0big ))

At (E_2) we have (p^*=0) and (i^*=i_2^*>0). The predator linearization is

Hence, a necessary condition for stability is

The remaining (2times 2) prey–disease block of (J(E_2)) can be written explicitly; the corresponding eigenvalues must also satisfy the Matignon sector condition. In particular, when both eigenvalues of this block are real and negative, (E_2) is locally asymptotically stable.

Stability of the coexistence equilibrium (E_3)

For (E_3=(s_3^*,i_3^*,p_3^*)) with all components positive, (J(E_3)) generally has one real eigenvalue and a complex pair. Local asymptotic stability is verified by computing the eigenvalues of (J(E_3)) and applying Theorem 7.1:

Equivalently, if a conjugate pair (rho pm iomega) occurs, one may use the practical inequality

as stated in Corollary 7.3. The stabilizing influence of smaller (nu) follows from Remark 7.4^18,19.

Fractional Lyapunov direct method

Definition 8.2

A function (V:Omega subset mathbb {R}^3rightarrow mathbb {R}) is called a Lyapunov function for the Caputo system ({}^{C}D_t^nu x=F(x)) if V is continuous, positive definite, and its Caputo derivative along solutions satisfies ({}^{C}D_t^nu V(x(t))le 0).

Lemma 8.3

(Fractional Lyapunov inequality) Let V(x) be continuously differentiable. Along solutions of ({}^{C}D_t^nu x=F(x)),

Theorem 8.4

(Fractional Lyapunov stability) Let (x^*) be an equilibrium of ({}^{C}D_t^nu x=F(x)). Assume there exists a continuously differentiable function V(x) and class-(mathcal {K}) functions (alpha _1,alpha _2,alpha _3) such that

and

in a neighborhood of (x^*). Then (x^*) is (locally) asymptotically stable.

A practical Lyapunov candidate (template for global analysis)

For equilibrium points with positive coordinates, a common choice is the logarithmic Lyapunov function

where (c_1,c_2,c_3>0). It satisfies (Vge 0) and (V=0) iff ((s,i,p)=(s^*,i^*,p^*)). Using Lemma 8.3,

and the right-hand side can be estimated to obtain sufficient conditions for global stability (by selecting (c_1,c_2,c_3) appropriately). The detailed algebra depends on which equilibrium ((E_1), disease-free coexistence, or endemic coexistence) is targeted.

Remark 8.5

For the disease-free equilibrium ((s^*,i^*,p^*)=(k,0,0)), a quadratic Lyapunov candidate in ((s-k,i,p)) is often more convenient because the logarithmic form requires strictly positive coordinates.

Explicit fractional adjoint system (fully expanded)

For convenience, define

The controlled fractional system (6.2) can be written as

where

The Hamiltonian is

Useful partial derivatives of G.

Derivatives of the ratio term.

Since (D=a s+b i),

Also,

Right-sided adjoint equations.

For (0<nu <1), the fractional PMP gives (right-sided RL derivative):

Now compute (partial mathcal {H}/partial s), (partial mathcal {H}/partial i), (partial mathcal {H}/partial p) explicitly.

(1) Equation for (psi _1)

(2) Equation for (psi _2)

(3) Equation for (psi _3)

Control characterization (same projection laws).

The optimal controls remain

Let (T>0) and choose a uniform grid (t_n=n h), (n=0,1,dots ,N), with (h=T/N). We present an L1-type discretization for the Caputo state equations and a compatible right-sided L1 discretization for the adjoint equations.

L1 scheme for the Caputo derivative (forward in time)

For (0<nu <1), the Caputo derivative at (t_{n+1}) is approximated by the L1 formula:

Apply (9.1) to s, i, p and solve sequentially for ((s_{n+1},i_{n+1},p_{n+1})) from the nonlinear algebraic system

etc. A practical implementation uses a one-step fixed-point iteration (or Newton) at each time step.

Right-sided L1 scheme for the adjoint (backward in time)

The right-sided RL derivative can be discretized on the same grid (backward sweep). Define (psi _{j,n}approx psi _j(t_n)). A compatible discrete right-sided operator is

with terminal condition (psi _{j,N}=0).

Then, for (n=N-1,dots ,0), compute (psi _{1,n},psi _{2,n},psi _{3,n}) from

where the right-hand sides are evaluated using (8.8)–(8.10) at ((s_n,i_n,p_n,u_{1,n},u_{2,n},psi _n)).

Algorithm (forward–backward sweep with relaxation)

1. 1.

   Initialization: Choose initial guesses (u_1^{(0)}(t_n),u_2^{(0)}(t_n)) (e.g., zeros). Set iteration counter (m=0) and tolerance (varepsilon>0).

2. 2.

   Forward step (states): Using (u^{(m)}), solve (6.2) on (t_0rightarrow t_N) with the L1 scheme (9.1) to obtain ({s_n^{(m)},i_n^{(m)},p_n^{(m)}}_{n=0}^N).

3. 3.

   Backward step (adjoints): Set (psi _{1,N}=psi _{2,N}=psi _{3,N}=0). Using the discrete right-sided operator (9.2) and the explicit formulas (8.8)–(8.10), compute ({psi _{1,n}^{(m)},psi _{2,n}^{(m)},psi _{3,n}^{(m)}}_{n=0}^{N}) backward from (n=N-1) to 0.

4. 4.

   Control update: For each grid point (t_n), update

   $$tilde{u}_{1,n}=operatorname {Proj}_{[0,1]}!left( frac{lambda s_n^{(m)} i_n^{(m)}}{B_1}big (psi _{2,n}^{(m)}-psi _{1,n}^{(m)}big )right) , qquad tilde{u}_{2,n}=operatorname {Proj}_{[0,1]}!left( frac{psi _{2,n}^{(m)} i_n^{(m)}}{B_2}right) .$$

   Use relaxation ((0<omega <1)) to stabilize iterations:

   $$u_{1,n}^{(m+1)}=omega u_{1,n}^{(m)}+(1-omega )tilde{u}_{1,n},qquad u_{2,n}^{(m+1)}=omega u_{2,n}^{(m)}+(1-omega )tilde{u}_{2,n}.$$
5. 5.

   Stopping rule: If

   $$max _nleft( |u_{1,n}^{(m+1)}-u_{1,n}^{(m)}|+|u_{2,n}^{(m+1)}-u_{2,n}^{(m)}|right) <varepsilon ,$$

   stop; otherwise set (mleftarrow m+1) and repeat.

Remark 9.1

Because the L1 scheme yields implicit nonlinear algebraic equations at each time step, a one- or two-iteration fixed-point solver per step is often sufficient in practice. If stiffness occurs, Newton’s method or smaller h is recommended.

We study how the fear factor (k_1) affects the qualitative dynamics of the fractional ecoepidemiological system (4.1), and we perform bifurcation analysis with respect to the conversion efficiency (theta), the predator hunting efficiency toward infected prey (beta), and the infection rate (lambda). The nonconsumptive fear effect acts through the prey reproduction reduction factor (frac{1}{1+k_1p}); see^6,7,9,10.

How fear enters the model and why it matters

In (4.1), the fear factor (k_1) appears only in the susceptible prey growth term:

Hence, for any fixed s, i and (p>0),

so increasing (k_1) decreases the effective prey reproduction rate. This reduces prey availability, indirectly affecting (i) infection transmission through (lambda s i) and (ii) predator recruitment through the predation gains in the p-equation.

Qualitative consequences of increasing (k_1).

• Reduced susceptible prey recruitment (Rightarrow) smaller prey biomass available to sustain predators.

• Reduced disease incidence (lambda s i) (through reduced s).

• Possible predator extinction (coexistence equilibrium may lose feasibility when prey growth becomes too weak).

• Shift of bifurcation thresholds in (lambda ,beta ,theta) because equilibria and Jacobian eigenvalues depend on (k_1).

Equilibria and (k_1)-dependent feasibility

Equilibria (E=(s^*,i^*,p^*)) satisfy (f_1=f_2=f_3=0), where ((f_1,f_2,f_3)) are the right-hand sides of (4.1).

(i) Extinction equilibrium

(ii) Disease-free prey-only equilibrium

Setting (i^*=0) and (p^*=0) gives (s^*=k), thus

(iii) Predator-free endemic equilibrium

Assume (p^*=0) and (i^*>0). From (f_2=0),

Substituting into (f_1=0) with (p^*=0) yields a positive (i^*) if and only if (lambda k>mu _1) (equivalently (mathcal {R}_0=lambda k/mu _1>1)), consistent with standard epidemic thresholds^1,2. Denote this equilibrium by

(iv) Interior coexistence equilibrium

Let (E_3=(s_3^*,i_3^*,p_3^*)) with (s_3^*,i_3^*,p_3^*>0). Using (f_2=0) and (f_3=0), we obtain the useful identities

Equation (10.2) shows explicitly how (theta) and (beta) regulate predator persistence at coexistence.

How (k_1) constrains coexistence.

At (E_3), the susceptible equilibrium equation becomes

The left-hand side is strictly decreasing in (k_1) whenever (p_3^*>0). Therefore, increasing (k_1) can destroy feasibility of (E_3) beyond a critical fear level (k_1^textrm{crit}) (where (E_3) collides with a boundary equilibrium or ceases to exist). This is a fear-induced feasibility loss.

Local stability framework for fractional systems (Matignon)

Let J(E) be the Jacobian of the vector field at equilibrium E. For the fractional-order system with (0<nu <1), local asymptotic stability is governed by Matignon’s criterion:

based on Mittag–Leffler asymptotics^18,19. Hence, bifurcations occur when eigenvalues meet the sector boundary.

Hopf-type boundary in fractional-order.

In integer order ((nu =1)), Hopf occurs at (Re (lambda )=0) for a complex pair. In fractional-order ((0<nu <1)), the corresponding oscillatory boundary is

Thus, varying (k_1), (theta), or (beta) can induce/destroy oscillations when the complex pair crosses this sector boundary.

Bifurcation in infection rate (lambda) (disease invasion/transcritical)

Consider dynamics near the disease-free manifold (iapprox 0). From the infected equation,

Hence, linearizing at a disease-free equilibrium (E=(bar{s},0,bar{p})) gives the infected-direction eigenvalue

Therefore, the invasion threshold occurs when (lambda _I=0), i.e.

This describes a transcritical-type bifurcation between a disease-free equilibrium and an endemic equilibrium branch, analogous to classical host–parasite thresholds^1,2.

Effect of fear on (lambda _c).

Because ((bar{s},bar{p})) depends on (k_1) through the term (frac{1}{1+k_1p}) in prey growth, fear shifts (lambda _c). Numerically, stronger fear typically reduces (bar{s}) and often reduces (bar{p}), which can increase (lambda _c) and make invasion harder, highlighting a stabilizing role of fear.

Bifurcation in conversion efficiency (theta) (predator persistence threshold)

At coexistence, (10.2) implies

Thus, for fixed ((s^*,i^*)) the predator can persist only if

Although ((s^*,i^*)) depends on parameters, (10.6) shows that increasing (theta) increases predator recruitment, which can create a predator persistence bifurcation (often transcritical with a predator-free equilibrium).

Interaction with fear.

Since larger (k_1) reduces prey reproduction and tends to reduce (s^*) (and may reduce (i^*)), it modifies the denominator in (10.6), thereby shifting (theta _c) upward. Hence, strong fear can require larger conversion efficiency for predator survival.

Bifurcation in hunting efficiency (beta) (predator gain vs infection suppression)

The parameter (beta) affects the model in two competing ways:

• it increases removal of infected prey via (-beta i p) (reducing infection),

• it increases predator growth via (+theta beta i p) (supporting predators).

At equilibrium, (10.1) implies

so increasing (beta) decreases the predator level required to balance infection at fixed ((lambda ,s^*)). Meanwhile (10.2) shows predator persistence gains directly through (beta i^*).

Thus, varying (beta) may produce:

• a predator persistence bifurcation (predator-free (leftrightarrow) coexistence),

• an oscillatory (Hopf-type) bifurcation through the Jacobian complex pair crossing the Matignon sector boundary.

Fear modifies these thresholds because it changes (s^*,i^*,p^*) through (10.3).

Recommended numerical bifurcation procedure (one- and two-parameter)

Because closed-form coexistence equilibria are generally intractable, we recommend numerical continuation:

Step 1 (equilibria).

For each chosen parameter value (e.g. (k_1) as primary), solve the steady-state equations (f_1=f_2=f_3=0), retain only feasible equilibria ((s^*,i^*,p^*ge 0) and (a s^*+b i^*>0)).

Step 2 (stability).

Compute the Jacobian J(E) numerically and its eigenvalues. Apply Matignon:

to classify each equilibrium^18,19.

Step 3 (one-parameter diagrams).

Construct bifurcation diagrams for:

Track equilibrium branches and detect:

• transcritical points via (lambda _I=0) in (10.4),

• Hopf-type points via (arg (lambda _{2,3})=pm nu pi /2),

• saddle-node points (turning points/feasibility loss).

Step 4 (two-parameter planes).

To highlight fear interactions, compute two-parameter bifurcation sets:

marking stability regions of disease-free, endemic, predator-free, and coexistence states. This approach aligns with modern fear/cost-of-avoidance bifurcation studies^11,12,13,39.

Remark 10.1

For (nu =1), classical continuation packages (MATCONT/XPPAUT) can be used. For (0<nu <1), equilibria and Jacobian eigenvalues are computed as usual, and stability is decided by Matignon’s criterion; time-domain verification can be carried out using fractional numerical solvers^19,20.

In this section, numerical simulations are carried out to illustrate the analytical results and to investigate the complex dynamics of system (4.1). In particular, the influence of the fear factor (k_1) and the effects of varying the conversion efficiency (theta), the hunting efficiency (beta), and the infection rate (lambda) are examined see³⁷.

All simulations are performed using the parameter values adopted in the final section of³⁷, with initial conditions

Unless otherwise stated, the parameter set is

Effect of fear factor (k_1)

Figure 3 illustrates the bifurcation structure of the infected prey population with respect to the fear factor (k_1). For small values of (k_1), the system exhibits stable periodic oscillations. As (k_1) increases, the extrema of the infected prey density split and spread, indicating the onset of more complex dynamics.

This behavior confirms that predator-induced fear plays a crucial role in regulating ecoepidemiological interactions by suppressing prey reproduction and indirectly modifying infection persistence.

Phase portraits and dynamical transitions

To further illustrate the effect of the fear factor, three-dimensional phase portraits in the (s, i, p)-space are depicted in Fig. 4. For (k_1=0.75), the system approaches a regular periodic orbit. As (k_1) increases to 0.82 and 1.2, the trajectories exhibit larger oscillations and more intricate geometric structures, indicating stronger nonlinear feedback between the prey, infected prey, and predator populations.

Bifurcation with respect to (theta), (beta), and (lambda)

Figure 3 also shows bifurcation sketches with respect to the conversion efficiency (theta), the hunting efficiency (beta), and the infection rate (lambda). Increasing (theta) enhances predator biomass conversion, which can destabilize steady states and induce oscillatory behavior. Similarly, changes in (beta) significantly alter predator–infection interactions, leading to qualitative transitions in the infected prey dynamics.

Variation in the infection rate (lambda) demonstrates a threshold-type behavior: small values of (lambda) lead to disease extinction, while larger values result in persistent infection and oscillatory outbreaks. These numerical findings are consistent with the analytical bifurcation conditions derived in the previous sections.

Remark 11.1

The observed transitions are influenced by the combined effects of fear-induced prey suppression and predator efficiency. In particular, strong fear effects can delay or suppress oscillatory instabilities that would otherwise arise from increased (theta), (beta), or (lambda).

Fractional-order dynamics ((0<nu <1))

In this subsection, we investigate how the fractional-order (nu) influences the qualitative behavior of system (4.1). All parameters and initial conditions are kept identical to those used in the integer-order simulations, except that the Caputo derivative order is taken in the interval (0<nu <1).

Fractional-order systems incorporate memory effects, which are known to suppress fast oscillations and alter stability boundaries^18,19. Numerical integration is performed using a predictor–corrector Adams–Bashforth–Moulton scheme adapted to Caputo derivatives.

Impact of decreasing (nu).

Numerical experiments indicate that decreasing (nu) has a pronounced stabilizing effect on the system:

• oscillation amplitudes of all populations decrease as (nu) decreases,

• periodic solutions observed at (nu =1) may become asymptotically stable equilibria for smaller (nu),

• complex or chaotic-like oscillations are significantly damped.

In particular, parameter sets that produce sustained oscillations or complex dynamics in the integer-order case often lead to regular periodic behavior or convergence to equilibrium when (nu <1). This behavior is consistent with the Matignon stability criterion, since decreasing (nu) enlarges the stability sector in the complex plane¹⁸.

Interaction between fear and memory.

The combined presence of fear-induced prey suppression ((k_1>0)) and memory effects ((nu <1)) further enhances stability. Strong fear reduces prey availability, while fractional memory smooths transient responses, leading to faster convergence and lower oscillatory intensity.

These observations confirm that fractional dynamics provide an additional regulatory mechanism in ecoepidemiological systems beyond classical predator–prey feedback.

Cost-effectiveness analysis

To quantify the cost-effectiveness of the fractional optimal control strategy, we evaluate the objective functional

for different values of the fractional-order (nu in (0,1]). For each (nu), the optimality system (state–adjoint–control) is solved using the fractional forward–backward sweep method, and then (J(nu )) is computed numerically using a quadrature rule (e.g., composite trapezoidal rule).

Decomposition of the objective.

It is useful to split (11.1) into interpretable parts:

The term (J_I) measures the total infection burden over [0, T], while (J_{u_1}) and (J_{u_2}) measure the energetic/economic costs of prevention and treatment.

Observed effect of the fractional-order.

Numerical results typically show that smaller values of (nu) lead to smaller total costs (J(nu )). This can be interpreted as a memory-driven stabilization: decreasing (nu) increases the influence of past states, which tends to damp oscillations and reduce peak infection levels. Consequently, the infection burden term (J_I(nu )) decreases, and the control costs may also decrease because less aggressive interventions are needed to suppress infection.

This highlights the importance of incorporating fractional dynamics in ecoepidemiological optimal control problems, since the fractional-order (nu) acts as an additional mechanism that can change both the epidemiological outcome and the required intervention effort.

Figure 5 illustrates the effect of the fractional-order on the total cost and its components.

In this work, a fear-induced ecoepidemiological predator–prey model was extended to a fractional optimal control framework by incorporating intervention strategies aimed at reducing disease transmission and enhancing the removal of infected prey. The inclusion of predator-induced fear follows recent developments showing that nonconsumptive predator effects can substantially alter population dynamics and disease persistence^6,7,9,10. By formulating the model using Caputo fractional derivatives, memory effects inherent in ecological and epidemiological processes were explicitly accounted for^18,19,20.

A fractional optimal control problem was formulated, and necessary optimality conditions were derived using the fractional Pontryagin Maximum Principle. The resulting adjoint system and explicit characterizations of the optimal controls were obtained in a projected form consistent with existing fractional optimal control theory^{23,24,25,26,27}. These analytical results provide a rigorous foundation for designing cost-effective intervention strategies in ecoepidemiological systems with memory.

Numerical implementation of the optimality system using a forward–backward sweep algorithm allows the generation of time series, phase portraits, and optimal control profiles, and enables systematic comparisons between controlled and uncontrolled dynamics. The simulations demonstrate that, under parameter regimes where the baseline system exhibits sustained oscillations or complex behavior, the proposed control strategies can effectively suppress infection and stabilize population dynamics. Moreover, decreasing the fractional-order was observed to reduce oscillation amplitudes and overall control effort, confirming the stabilizing and cost-reducing role of fractional memory effects^18,20,30,32.

This study is based on a deterministic fractional-order mathematical model with assumed parameter values. Real ecological systems are subject to environmental variability, stochastic disturbances, and measurement uncertainty, which are not included in the present formulation. Future work may incorporate stochastic effects, spatial diffusion, and empirical parameter estimation to enhance biological realism. Future work may extend the present framework in several directions. Possible generalizations include the consideration of variable-order or non-singular kernel fractional derivatives^31,40, time-delay effects in disease transmission, spatial diffusion leading to fractional reaction–diffusion systems, and the incorporation of state or mixed control constraints. Additionally, data-driven parameter estimation and sensitivity analysis could enhance the biological relevance of the model. Such extensions would further improve the applicability of fractional optimal control approaches to realistic ecoepidemiological and environmental management problems.

Take-home message.

The present study demonstrates that incorporating predator-induced fear and fractional-order memory into ecoepidemiological models fundamentally alters system dynamics and control outcomes. Fractional dynamics not only enhance stability and suppress complex oscillations, but also significantly reduce the cost and intensity of optimal intervention strategies. These findings emphasize that fractional-order modeling provides a more realistic and effective framework for understanding and managing disease dynamics in predator–prey systems.

Future research directions

An interesting extension of the present work is the study of discrete-time fractional-order ecoepidemiological predator–prey models see⁴¹. Discrete predator–prey systems are known to exhibit complex dynamical behaviors, including bifurcation structures and stability transitions. Incorporating fractional memory effects into discrete frameworks may lead to richer dynamics and new theoretical challenges. Developing such models and analyzing their stability and bifurcation properties remains an important direction for future research.