Prey-predator interactions play a central role in ecological systems and form an important area of research in mathematical ecology. Various mathematical models have been developed to understand the dynamics of interacting species and the factors afecting their coexistence and stability [1, 2].

Models involving multiple predator species have received considerable attention because they help explain competition and resource sharing mechanisms among predators. Tansky [3] investigated the efects of predator switching, while Hsu and Hubbell [4] analyzed competition between predators sharing common prey resources. Kirlinger [?] established conditions for the permanence of species in a two-predator two-prey system. Subsequently, Bhattacharyya and Mukhopadhyay [5] incorporated spatial migration and predator switching into prey-predator dynamics and Hadziabdic et al. [6] further examined the behavior of Lotka-Volterra systems with two predators and one prey.

Prey refuge is another important ecological mechanism that allows prey populations to avoid predation by utilizing protected habitats. Ma et al. [7] demonstrated that prey refuges significantly influence predator-prey interactions and functional responses. Ko and Ryu [8] showed that the incorporation of a constant prey refuge can enhance the stability of predator-prey systems. Li et al. [9] further analyzed a predator-prey model with prey refuge and discussed its efects on the long-term dynamics and stability of populations. Recent studies have also explored the combined efects of refuge and other ecological factors on predator-prey interactions [10, 11].

Harvesting is another factor that strongly influences population dynamics. In natural and managed ecosystems, harvesting rates often vary with environmental conditions, seasonal changes and human activities. Dieci and Rebaza [12] investigated dynamical transitions in nonlinear systems, while Leard et al. [13] studied ratio-dependent predator-prey models with non-constant harvesting and demonstrated its impact on population stability and persistence.

Motivated by these studies, the present work proposes a prey-predator model consisting of two predator species, a prey refuge for the first prey population, and non-constant harvesting of the specialist predator. The objective is to investigate the dynamical behavior of the system and analyze the stability characteristics of the interacting populations.

The prey populations follow logistic growth and grow rapidly when small but slow down near their carrying capacities. They interact with predators based on prey availability. A prey refuge protects a fixed proportion of the population, leaving only exposed to predation. The system includes a specialist predator , which feeds on a limited prey type, and a generalist predator , which consumes a wider range of prey. It also incorporates non-constant harvesting, where harvesting rates vary with time, population size, or environmental conditions. The corresponding mathematical formulation and its interaction represention are given below:

The positivity and boundedness properties ensure that all species populations remain non-negative and finite over time. Ecologically, this reflects that populations cannot be negative and are limited by resources and environmental constraints, keeping the model biologically realistic.

Theorem 3.1. All solution of the system (1) initiating from is positive for .

Proof. The first equation of (1) gives

0. \end{document} ]]>

The second equation of the system (1) is

0. \end{document} ]]>

From the third equation of the system (1), we get

0. \end{document} ]]>

The last equation of system (1) gives

0. \end{document} ]]>

This proves the theorem.

Theorem 3.2. All solutions of system (1) in are uniformly bounded.

Proof. Let us define a function as follows:

Differentiate equation (1) with respect to , we get

,

.

For any 0 \end{document} ]]>.

.

we define a constant 0 \end{document} ]]> such that . It

shows that . Applying theory of differential inequality, we obtain .

As , this become . Hence, that concludes all solution of (1) with initial condition that initiate in are bounded to the region .

Various equilibrium points of the system have been examined. The system is observed to have six possible equilibrium points. and , all the equilib-rium points satisfy the condition. These equilibrium points are now presented as follows:

(i)    Trivial equilibrium represents the complete extinction of all prey and predator populations under extreme ecological stress.

(ii)  Coexistence of prey without predator equilibrium represents the persistence of both prey species in the absence of predators.

(iii)  Single prey with dual predator coexistence equilibrium

occurs when

c_2 \end{document} ]]>, \frac{l_3w_1n_2\mathcal{J}}{l_2c_2} \end{document} ]]>and \frac{c_2j}{v_1(1-\hat{m})(n_2l_1-c_2)} \end{document} ]]>ensuring that one prey species persists and sustains both predator populations while the other prey is absent.

Where

(iv)  Generalized predator absent equilibrium

exists when c_1+\hat{h}_1 \end{document} ]]>andT \end{document} ]]>under which both prey coexist with the specialist predator while the generalist predator is absent.

(v)  Coexistence equilibrium point exists when all population components are positive, representing the simultaneous persistence of both prey and both predator species.

The equilibrium points are presented below:

The dynamics behavior of the equilibrium points can be analyzed by calcu-lating the eigenvalue of the Jacobian matrix of system (1) given by:

where

,

Theorem 5.1. The trivial equilibrium point of model (1) is a saddle point.

Proof. The jacobian matrix is,

The eigenvalues are λ1 = w1, λ2 = w2, λ3 = −c1 − hˆ1 and λ3 = −c2.

Thus, being a saddle point, it is unstable.

Theorem 5.2. Equilibrium point of model (1) is stable if the condition

\dfrac{n_1 v_1 l_2}{v_2 + j} \quad \end{document} ]]>and \dfrac{n_2 l_1 (1-\hat{m}) v_1}{v_1 (1-\hat{m}) + j} + \dfrac{n_3 l_3 v_2}{v_2 + j} \quad \end{document} ]]>hold.

Proof. The jacobian matrix is,

The eigenvalue of

In the case where \dfrac{n_1 v_2 l_2}{v_2+j} \quad \text{and} \quad c_2 > \dfrac{n_2 l_1(1-\hat{m})v_1}{v_1(1-\hat{m})+j} + \dfrac{n_3 l_3 v_2}{v_2+j} \end{document} ]]>the Equilibrium point is stable.

Prey populations persist while predators cannot invade. This reflects a sit-uation where prey refuges and predator mortality prevent predator establishment, allowing prey to maintain a stable population.

Theorem 5.3. Single prey with dual predator coexistence equilibrium point  is stable if the condition and hold.

Proof. The jacobian matrix is,

Where

The characteristic equation is

Where,

Here,

In the case where and is stable.

Under these conditions, the single prey species and both predators coexist stably. The inequalities ensure that predator growth is balanced by prey availability and mortality, preventing any species from overexploiting the prey. Ecologically, this reflects a stable prey predator system where interactions are regulated and no population collapses occur.

Theorem 5.4. Equilibrium point  is stable if the condition hold.

Proof. The jacobian matrix is,

Where

The characteristic equation is

Where

Here

In the case where is stable. This point represents sustainable coexistence where all species persist and interactions among them prevent extreme population fluctuations.

Theorem 5.5. Coexistence Equilibrium point of model (1) stable if thecondition 0 \end{document} ]]> and 0. \end{document} ]]>

Proof.

Where

The characteristic equation is

Where,

In the case where 0 \end{document} ]]>

and 0 \end{document} ]]> is stable.

This equilibrium represents a fully stable coexistence of all species. Predator and prey populations are regulated by their interactions, growth and mortality (or harvesting), so no species ove exploits another. Small disturbances in populations decay over time, maintaining the balance of the ecosystem.

The proposed model is not based on real ecological data or a specific case study therefore all parameter values used in the numerical simulations are purely hypothetical and chosen for illustrative purposes. The values are selected to be positive and within reasonable ranges to ensure biologically feasible solutions and to satisfy the analytical conditions derived in this section. Different sets of parameters are employed in Figures 1, 2, and 3 to illustrate distinct dynamical behaviors of the system and to demonstrate the qualitative validity of the theoretical results. All simulations are performed using MATLAB R2024b. The values are,

Set 1 represents the case where the harvesting rate is less than the specialist predator population. Set 2 denotes values where the harvesting rate is higher than the population. Set 3 indicates no harvesting occurs in the specialist predator popuatin.

This study analyzes a mathematical model with two prey species and two distinct predator types, in which the first prey species k1 is protected by a refuge. The dynamics of the specialist predator are influenced by a non constant harvest-ing strategy. Figure 1 shows the interaction structure, while Figure 2 presents the equilibrium points of the system. When the harvesting rate is lower than the specialist predator population, the predator population increases even as the prey k2 declines (Figure 3). Conversely, a higher harvesting rate reduces predator den-sity (Figure 4). In the absence of harvesting, the specialist predator population rises while the prey population decreases (Figure 5). Non constant harvesting allows the rate to vary with population density. Simulation results align with stability theory, as population trajectories approach the predicted equilibrium points and confirming their local stability. This adaptive harvesting reduces pressure at low population levels and facilitates recovery. However, the model assumes fixed parameter values and does not consider spatial heterogeneity, which may limit its applicability. Fu-ture research could extend the model by incorporating empirical data and spatial dynamics with fear effects and anti predator behavior to enhance ecological realism and predictive power.