Dynamics of a Predator-Prey Population in the Presence of Resource Subsidy under the Influence of Nonlinear Prey Refuge and Fear Effect

In this work, our aim is to investigate the impact of a non-Kolmogorov predator-prey-subsidymodel incorporating nonlinear prey refuge and the effect of fear with Holling type II functional response. )e model arises from the study of a biological system involving arctic foxes (predator), lemmings (prey), and seal carcasses (subsidy). )e positivity and asymptotically uniform boundedness of the solutions of the system have been derived. Analytically, we have studied the criteria for the feasibility and stability of different equilibrium points. In addition, we have derived sufficient conditions for the existence of local bifurcations of codimension 1 (transcritical and Hopf bifurcation). It is also observed that there is some time lag between the time of perceiving predator signals through vocal cues and the reduction of prey’s birth rate. So, we have analyzed the dynamical behaviour of the delayed predator-prey-subsidy model. Numerical computations have been performed usingMATLAB to validate all the analytical findings. Numerically, it has been observed that the predator, prey, and subsidy can always exist at a nonzero subsidy input rate. But, at a high subsidy input rate, the prey population cannot persist and the predator population has a huge growth due to the availability of food sources.


Introduction
In the ecological system, the predator-prey interaction is one of the most significant tools which is comparatively easy to observe in the field. But fear of the predator felt by the prey (indirect effect) also plays a vital role since its effect is stronger than direct predation [1,2]. e cost of fear can reduce the reproduction rate of prey because it affects the physiological condition of prey population. As a result, the prey species may get a long-term loss. In support of this, it is mentioned that, in the Greater Yellowstone Ecosystem, wolves (Canis lupus) affect the reproductive physiology of elk (Cervus elaphus) [3]. When the prey species recognize the predator signal (chemical/vocal), they spend more time to become keenly watchful to detect danger rather than in foraging. So, the birth rate of the scared prey reduces and adopts some survival mechanisms like starvation [1,2]. For examples, some birds react to the sound of predator with antipredator defenses [1,2] and they flee from their nests at the first sign of danger [2]. is antipredator behaviour may affect survival and reproduction of the birds [2]. It has been experimentally investigated that, in the absence of direct killing, the reproduction of the offspring of song sparrows (Melospiza melodia) could be reduced by 40% as a result of impact of feeling fear created by the predator [4]. So, this reduction caused by the antipredator behaviour affects the birth rate and survival of offspring. us, the cost of fear (apart from direct predation) should be introduced in a predator-prey interaction. Mathematical formulation of the impact of fear on the two species prey-predator system has been initiated by Wang et al. [5] in 2016 introducing fear factor: f(k, y) � (1/(1 + ky)). It involves a parameter k denoting the level of fear to represent the antipredator behaviour of the prey. Some research works have already been done on the ecological system under the influence of predation fear felt by prey species [6][7][8][9][10][11][12][13][14][15][16][17]. Moreover, the impact of fear in a two-species predator-prey model with prey refuge was analyzed by many researchers [11,18,19].
In evolutionary biology, prey refuge is a concept which helps an organism to protect themselves from predation by hiding in an area inaccessible to the predator, for example, in a wolf-ungulate system, ungulates may seek refuge by migrating to areas outside the core territories of wolves. Also, it has many significant roles on the dynamics of predator-prey interactions: prey refuge may decrease the chance of extinction of prey. Researchers have mainly used the dynamic nature of predator-prey model with linear prey refuge (that is, mx amounts of prey are unavailable to the predator, where m ∈ (0, 1) is the coefficient of refuge and x is the biomass of prey species) with Volterra response [20][21][22]. Recently, Mondal and Samanta have studied the dynamics of the predator-prey system with prey refuge dependent on both species (that is, mxy amounts of prey are free from the predator risk, where 0 < 1 − my < 1, m is prey refuge coefficient, x is the biomass of prey population, and y is the biomass of predator population) in the presence of additional food (for details, see [23]). In 2020, Mondal and Samanta [11] have also analyzed the dynamics of predatorprey interaction having nonlinear prey refuge function Φ(x, y) � (mxy/(a + y)) which is the amount of prey that are free from predation, where a is half saturation constant and y is the biomass of predator population.
Many experimental studies suspect that the introduction of resource subsidies may disrupt otherwise stable food web linkages [24][25][26]. Such concept is significant for resource management purposes. It is learnt that reintroduced wolves in Yellowstone Park switch to bison when their preferred ungulate prey, namely, elk, are rare in the concerned ecosystem [27]. Mathematically, the influence of resource subsidy on the predator-prey model has been initiated in the work of Nevai and Van Gorder [28]. ey have discussed how different subsidy input rate may affect the prey and predator population to persist in the ecosystem.
Generalist predator can consume more than one food source: either multiple prey population or a combination of prey population and resource subsidy. ere are many rich theoretical research on ecological systems involving generalist predator [29][30][31]. Also, there are a variety of real-life applications for such systems [32][33][34]. From literature surveys, it has been shown that generalist predator can persist in an ecosystem even if one particular prey species is going towards extinction [29][30][31].
In 2012, Nevai and Van Gorder [28] first extended the Kolmogorov model to a non-Kolmogorov predator-preysubsidy model. It has been observed that the predator-preysubsidy model occurs in the arctic foxes (predator), lemmings (prey), and seal carcasses (subsidy). Motivated by the works of Das and Samanta [35], Nevai and Van Gorder [28], and Xu et al. [36], we have analyzed the dynamical behaviour of a mathematical model of non-Kolmogorov form that includes the three components (predator, prey, and subsidy) with the impacts of nonlinear prey refuge function and the fear effect felt by the prey in the presence of the predator. To the best of our knowledge, there does not exist any mathematical model to explore the impact of fear effect incorporating nonlinear prey refuge function in predator-prey-subsidy interaction. e organization of this work is structured as follows: in Section 2, a mathematical model has been formulated with the influences of nonlinear prey refuge and fear effect. Section 3 shows that the proposed model is well-behaved. In Section 4, feasibility criteria and stability of all the equilibria of the proposed system (in absence of delay) have been studied.
e equilibria can change their stability nature through transcritical and Hopf bifurcation which are also analyzed in this section. Generally, the reduction of prey's birth rate due to the effect of fear will not be an instantaneous biological process but deviated through some time lag, so the study of time-delay τ is very meaningful to obtain the more realistic dynamics. So, Section 5 deals with the dynamic behaviour of the delayed system for two equilibrium points E 3 (subsidy free) and E * (interior), respectively. Section 6 provides the numerical computations which support the analytical calculations. Section 7 provides a brief conclusion about the system dynamics.

Model Formulation
In 2020, Mondal and Samanta [11] analyzed the dynamics of a delayed predator-prey interaction incorporating nonlinear prey refuge function under the influence of fear effect and additional food. Motivated by the work of Mondal and Samanta [11], we have first considered a predator-preysubsidy model with nonlinear prey refuge function where the prey and subsidy occur in the same habitat and they are both consumed by a single generalist predator according to the following differential equations: where x is prey population, w denotes the population of the subsidy, and y is the generalist predator which exploits both the prey and subsidy. For example, wolves (predator) consume both deer (prey) and salmon carcasses (subsidy) [37]. e term (1 − (my/(a + y)))x represents the quantity of prey available to the predator, i.e., (mxy/(a + y)) amounts of prey are free from predation risk where (mxy/(a + y)) is designated as nonlinear prey refuge function. Also, we have modeled the dynamics of a generalist predator with Holling type II [38][39][40][41] response function in the presence of nonlinear prey refuge function.
All parameters are positive (except A ≥ 0) and biologically meaningful. Parameters are described in Table 1.
Apart from direct consumption, feeling of fear among the individuals of the prey species in presence of predator is very common in predator-prey interaction which changes life-history, behavioural responses, and reproduction capability of prey species. In ecology, effect of fear is a common factor, but there does not exist any considerable attention to introduce the impact of fear in the mathematical modeling.
Experimental studies indicate that the feeling of fear among the individuals of the prey species in presence of predator reduces the prey's birth rate. So, birth rate of prey species r is multiplied by a monotone decreasing function f(k, y) � (1/(1 + ky)), where k( ≥ 0) is a level of fear [5].
e fear function f(k, y) satisfies the following conditions: (1) f(0, y) � 1: when there is no fear effect on the prey species, the birth rate of the prey is not reduced (2) f(k, 0) � 1: when there is no predator, the birth rate of the prey species is not reduced in the presence of fear effect (3) (zf(k, y)/zk) < 0: when fear effect increases, the birth rate of the prey reduces (4) (zf(k, y)/zy) < 0: when predator species increases, prey population reduces Our main focus is to analyze the dynamic nature of the predator-prey-subsidy model with the influence of nonlinear prey refuge and fear effect. So, system (1) can be modified in the following aspects: with initial conditions: roughout the analysis of this work, we have taken c 1 > c 2 which is biologically meaningful.

Positivity and Uniform Boundedness
Theorem 1. Every solution of system (2) with (3) uniquely exists and is positive for all t ≥ 0.
Proof. Case 1: if r > d 1 , from the first equation of (2), Let us take P � x + w + (y/c 1 ). Differentiating both sides with respect to t, we obtain Birth rate of the prey d 1 Natural death rate of the prey d 2 e subsidy decay rate a 1 Mortality rate due to intraspecific competition among the individuals of the prey population a 2 Consumption rate of the predator a 3 Maximum rate at which the predator consumes the subsidy b 1 Handling time assumed to be uniform over all food sources c 1 Conversion rate of the energy that the predator obtains from the target prey c 2 Conversion rate of the energy that the predator obtains from the subsidy A Subsidy input rate m Coefficient of prey refuge (m ∈ (0, 1)) d 3 Mortality rate of the predator a Half saturation constant for refuge function 4 Complexity Let en, Using the Gronwall inequality, we obtain us, all solutions of system (2) enter into the region: Case 2: if r < d 1 , from the first equation of (2) we obtain lim t⟶∞ x(t) � 0. Now, from the second and third equations of (2), we have For large t, en, Using Gronwall inequality, we obtain Hence, the theorem.

Equilibrium Points and Stability Analysis
and we get where x is a positive root of the equation: Here, and ) exists if c 2 a 3 > d 3 and y > 0.

Interior (Coexistence) Equilibrium
Point. Solving the following system of equations, 6 Complexity we can obtain E * (x * , w * , y * ) using the software MATH-EMATICA with the following existence conditions: ) (otherwise, predator population goes into extinction)

Theorem 5. Axial equilibrium point (prey only)
e Jacobian matrix J 3 (x, 0, y) at E 3 (x, 0, y) is as follows: where e characteristic equation corresponding to J 3 (x, 0, y) is expressed as where Theorem 6. Subsidy-free equilibrium point E 3 (x, 0, y) is locally asymptotically stable if b 11 < 0, b 33 < 0, and b 13 < 0. e Jacobian matrix J 4 (0, w, y) corresponding to E 4 (0, w, y) is given by where 8 Complexity e characteristic equation corresponding to J 4 (0, w, y) is expressed as where where e characteristic equation corresponding to where Complexity e Jacobian matrix J * (x * , w * , y * ) corresponding to E * (x * , w * , y * ) is as follows: where e characteristic equation corresponding to J * (x * , w * , y * ) is expressed as where

Local Bifurcations of Codimension 1 4.3.1. Transcritical Bifurcation
Proof. We apply Sotomayor's theorem [43] to prove the occurrence of a transcritical bifurcation around E 1 with d 1 as bifurcation parameter. For applicability of Sotomayor's theorem, exactly one of the eigenvalues of the Jacobian matrix at E 1 must be zero and other eigenvalues must have negative real parts. So, we need to fulfill the condition Compute Δ 1 , Δ 2 , and Δ 3 as follows: where F � (F 1 , F 2 , F 3 ) T and F 1 , F 2 , and F 3 are given by erefore, by Sotomayor's theorem [43], system (2) undergoes a transcritical bifurcation at Proof. Let us apply Sotomayor's theorem [43] to prove the occurrence of a transcritical bifurcation around E 2 with d 3 as bifurcation parameter. For applicability of Sotomayor's theorem, exactly one of the eigenvalues of the Jacobian matrix at E 2 must be zero and other eigenvalues must have negative real parts. e eigenvectors of J( Compute Δ 1 , Δ 2 , and Δ 3 as follows: where F � (F 1 , F 2 , F 3 ) T and F 1 , F 2 , and F 3 are given by 12 Complexity erefore, by Sotomayor's theorem [43], system (2) exhibits a transcritical bifurcation at d 3 � d [

Complexity 13
Proof. Proof is the same as in eorem 10.

Theorem 13. System (2) undergoes a transcritical bifurcation around
Proof. Proof is the same as in eorem 11.  (2) where the characteristic equation at E * is en, Hopf bifurcation theorem is stated as follows.
Theorem 14 (Hopf bifurcation theorem [44]). If A 1 (k), ∈ R for which the characteristic equation (50) has the following: Theorem 15. If E * (x * , w * , y * ) is locally asymptotically stable, then a Hopf bifurcation is exhibited around E * (x * , w * , y * ) when k passes through its critical value k [H] provided , we can write equation (50) as e roots of equation (51) are where p i (k) are real functions of k in an open neighborhood of k [H] for i � 1, 2, 3. Next, we verify the transversality condition: Putting λ(k) � p 1 (k) + ip 2 (k) in (59), we get Differentiating both sides with respect to k, we have Comparing real and imaginary parts from both sides, we obtain where and X 4 � 2 _ A 1 p 1 p 2 + _ A 2 p 2 . Multiplying (55) by X 1 and (56) by X 2 and then adding, we get At k � k [H] , Hence, this theorem is proved by virtue of eorem 14.  (2) where the characteristic equation of E 3 is and then Hopf bifurcation theorem is stated as follows.

Delayed Dynamical System
In biological point of view, many processes, both natural and man-made, include time-delay. e study of delay factor makes our system much more realistic than non-delayed system. Also, a delay differential equation reveals much more complicated dynamics than an ordinary differential equation (for details, see [10-13, 23, 45-49]).
In reality, after sensing the vocal cue, individuals of prey species take some time for assessing the predation risk. So, the effect of fear (felt by prey) of predator does not respond spontaneously on the birth rate of prey population; some time lag must be needed. In view of this fact, the predatorprey-subsidy interactions (2) can be modified as follows: e initial conditions are assumed as (i � 1, 2, 3) For biological feasibility: Let us linearize (60) using the following transformations: It leads to (63) e characteristic equation corresponding to (62) is where If τ ≠ 0, E 3 of system (60) is LAS provided equation (64) has no purely imaginary roots and it is also LAS for τ � 0. Further, it has been shown that stability nature of E 3 switches at τ � τ ′ * . Already, it has been derived that E 3 is LAS provided B 1 > 0, B 3 > 0, and B 1 B 2 > B 3 for τ � 0 (nondelayed system). Let us discuss if the real part of the roots of equation (64) gradually increases to reach zero and eventually turns to a positive value when τ increases.

Complexity
Equating respective real and complex parts from both sides, we obtain Now, let us examine whether equation (64) has purely imaginary roots or not. For this purpose, let us take q 1 ′ � 0. en, equations (67) and (68) become Eliminating τ from (69) and (70) (squaring and adding), we get Putting q ′2 2 � β, we have is is a cubic equation of β. It is noticed that L(∞) � ∞. So, equation (72) has exactly one positive real root if L(0) < 0, i.e., if L 2 3 < M 2 2 .
□ Now, linearize system (60) using the transformations Complexity e characteristic equation corresponding to (81) is where If τ ≠ 0, E * of system (60) is LAS provided equation (83) has no purely imaginary roots and it is LAS for τ � 0. Furthermore, it has to be noted that changes of stability occur at τ � τ * . Already, it has been discussed that E * is LAS when τ � 0 provided A 1 > 0, A 3 > 0, and A 1 A 2 > A 3 . Here, equation (83) is a transcendental equation, so it contains infinitely many eigenvalues. In this situation, we cannot apply the Routh-Hurwitz criteria to determine the stability of system (60). To understand the stability behaviour, our necessity is to check the sign of the real parts of the eigenvalues of equation (83). Now, putting λ � q 1 + iq 2 in equation (83), we have Equating respective real and complex parts from both sides, we get To check whether (83) has purely imaginary roots or not, set q 1 � 0; then, (86) and (87) become Eliminating τ from (88) and (89) (squaring and adding), we get Putting q 2 2 � σ, we have is is a cubic equation of σ. It is noted that R ′ (∞) � ∞. So, equation (91) has exactly one positive root if R ′ (0) < 0, i.e., if R 2 3 < S 2 2 . Let σ � σ + be a positive root of (91); then, q 2 � �� σ + √ .
Let us study the existence of Hopf bifurcation around E * with τ as bifurcation parameter.

Theorem 19.
Suppose E * exists and is locally asymptotically stable for system (2) when τ � 0. If R 2 3 < S 2 2 , then there exists a critical value τ * such that E * of system (60) is LAS when τ ∈ [0, τ * ) and unstable when τ > τ * , where and τ * � τ (0) + (minimum value). Also, a supercritical Hopf bifurcation is exhibited around E * at τ � τ * provided Proof. Proof is similar to that in eorem 18.  Figure 10: Nature of steady state E * when subsidy input rate A varies from 1 to 10 and other parameters are fixed as in Figure 9.
then there is a positive integer k such that and there are k switches from stability to instability to stability; that is, when then E * is locally asymptotically stable and when

Numerical Computations
Here, we have illustrated numerical simulations to verify the analytical findings of the proposed system (2). We select a parameter set: Under this set of parametric values, the stable nature of E 0 (0, 0, 0) is shown in Figure 1. If we take subsidy input rate A � 2 and other parametric values are chosen from the data set of Figure 1, then the subsidy only equilibrium E 1 (0, (A/d 2 ), 0) ≡ E 1 (0, 0.5, 0) exists and stable nature of E 1 (0, 0.5, 0) with time t is depicted in Figure 2. Now, we choose another parameter set: {r � 5.5, Under this set of parametric values, the prey only equilibrium E 2 ((r − d 1 )/a 1 , 0, 0) ≡ E 2 (8, 0, 0) exists and stable behaviour of E 2 (8, 0, 0) is presented in Figure 3. Let us choose the parameters as follows: If we take subsidy input rate A � 0 and other parameters are taken from set (100), then subsidy-free equilibrium point E 3 (x, 0, y) ≡ E 3 (1.20053, 0, 6.49267) exists and is locally asymptotically stable. Stable time series and stable phase diagram are represented in Figure 4. In the same manner, if we change the value of the parameter k(� 0.01) and others are the same as in the data set of Figure 4, then it is observed that E 3 (1.20053, 0, 6.49267) is unstable accompanied with a limit cycle (see Figure 5). From Figures 4 and 5, it can be easily noted that there must exist a threshold value of k, say k * � 0.019 for which unstable behaviour of E 3 changes to stable spiral. Since the vector fields for k < k * and k > k * are qualitatively different, a Hopf bifurcation is created around E 3 taking k as bifurcation parameter (see Figure 6). For the set of parameter values {r � 5.5, d 1 � 6.5, d 2 � 0.3, k � 0.2, d 3 � 0.2, a 1 � 0.6, a 2 � 0.98, a 3 � 0.8, b 1 � 2.5, c 1 � 0.85, c 2 � 0.7, A � 2, a � 1.1, m � 0.01}, prey free equilibrium point E 4 (0, w, y) ≡ E 4 (0, 1.3889, 5.5417) exists and is stable (see Figure 7). Next, let us take a different set of parameters of system (2): {r � 5.5, , 0) ≡ E 5 (8.3333, 6.6667, 0) is locally asymptotically stable. e stable behaviour with time t is shown in Figure 8.
If we take subsidy input rate A � 2 and others are fixed as in the data set of Figure  4, then E * (x * , w * , y * ) ≡ E * (0.385717, 0.950363, 8.34509) exists and is locally asymptotically stable. Figure 9 depicts the stable behaviour of E * . Comparing Figures 4 and 9, it is observed that subsidy input rate A enhances the value of y component of E * and decreases the value of x component of E * . Also, from Figure 10, it is noticed that the prey population is leading towards extinction and the predator population has enormous growth (due to huge supply of food source) at high subsidy input rate (when A ∈ (2, 10]) in the presence of fear felt by prey population. So, it can be concluded that it is not possible to control prey population from extinction in presence of nonlinear prey refuge because they cannot get enough time to protect themselves from  predation risk. After extinction of prey species, predator can easily survive with the help of resource subsidy. us, the parameter A has great importance in the proposed population dynamics.
Moreover, Figure 11 represents the unstable nature of E * when k � 0.02 and other parametric values are the same as in Figure 9. So, the parameter k has an interesting nature because there exists a threshold value k [H] � 0.025 of k for which unstable nature (limit cycle) of E * switches to stable behaviour (stable spiral) when k passes through its critical value k [H] ; i.e., the vector fields for k > k [H] and k < k [H] are topologically different. Hence, a Hopf bifurcation occurs around E * and Figure 12 depicts the corresponding bifurcation diagram taking k as bifurcation parameter. Also, it has to be noted from Figure 13 that, in the absence of fear effect, the oscillatory behaviour of E * changes to stable nature when subsidy input rate A crosses its critical value A * � 7.9 (approximately) and since predator population has huge growth rate at very large value of subsidy input rate A, the prey population cannot persist in ecosystem in presence nonlinear prey refuge. is phenomenon is very interesting because the prey refuge cannot control the prey population from extinction due to enormous growth of predator when subsidy input rate is very high. In this manner, system (2) is not persistent.
Further, Figure 14 depicts that the nature of steady states E 3 and E * when m ∈ (0, 1). Here, the predator population cannot go extinct for large value of coefficient of prey refuge parameter. Also, Figure 15 shows the changes of predator's growth at the steady states E 3 and E * for three different fear levels k when m varies from 0 to 1. Here, also the predator can persist for large m. From here, it may be concluded that system (2) is always persistent for small subsidy input rate in the presence of nonlinear prey refuge function. Again, Figures 16 and 17, respectively, show that, in the absence of fear effect (k � 0), the equilibria E 3 and E * are approaching towards stable state by excluding the existence of oscillatory behaviour taking m as the bifurcating parameter. In this manner, predator population also survives in ecosystem for large coefficient of prey refuge parameter m.
For the parameter set {r � 5.

Conclusion
We have analyzed a system for generalist predator which utilizes more than one food source: predator-prey-subsidy model of non-Kolmogorov form introducing nonlinear prey refuge function and the effect of fear felt by prey population. Our main interest is to find the situations such that dynamical stability and instability appear so as to make out more fully how subsidy may influence the predator and their prey. It has been shown that the solutions of system (2) remain positively invariant always and they are asymptotically uniformly bounded. ese, in turn, imply that system (2) is biologically well-behaved. Existence Complexity criteria and stable behaviour of all the biologically meaningful equilibria have been discussed. It has to be noted that Hopf bifurcations are exhibited around E 3 (subsidy free) and E * (interior) of system (2) considering k as a bifurcating parameter (see 9,11,and 12). Also, observing Figures 6 and  12, it can be concluded that high levels of fear can stabilize system (2) by excluding the existence of periodic solutions. ese phenomena are biologically significant because prey species are aware after a certain level of fear; i.e., after a certain level of fear, they are not affected as they are aware and show signs of habituation.
Moreover, this work derives transcritical bifurcations (local bifurcation of codimension 1) at the various equilibrium points E 1 , E 2 , E 4 , and E 5 , respectively (see . Also, we have discussed numerically the influences of coefficient of prey refuge parameter m on the nature of the equilibrium points E 3 (zero subsidy input rate) and E * (fixed small subsidy input rate) irrespective of fear level k. Noting Figures 14-17, it is observed that both the prey and predator species always persist in ecosystem due to continuous increment of coefficient of prey refuge. But Figures 10 and 13 depict that, irrespective of fear level, a highly subsidized predator should indeed drive the prey population towards extinction regardless of whether the prey and subsidy arise in the same habitat. is phenomenon is ecologically meaningful because the prey population cannot get enough time to protect themselves from predation risk for enormous growth of predator at high subsidy input rate. So, the prey population is leading towards extinction, but the predator species can easily survive in ecosystem with the help of resource subsidy. us, the study of system (2) is ecologically very significant.
In reality, fear effect does not instantaneously reduce the birth of a prey population, but some time lag should be needed to create an impact on the birth rate of the prey population. We have considered that there is a time-delay on the impact of fear to the birth rate of prey, from the instance it perceives the fear of predator through any means. So, the incorporation of time-delay makes system (60) more realistic. It is noted that delay parameter τ has a significant role because there exists a threshold value τ′ * such that stable behaviour of planer equilibrium point E 3 (in the absence of subsidy input rate) switches to oscillatory nature when τ passes through its threshold value τ ′ * ; i.e., the vector fields for τ < τ ′ * and τ > τ ′ * are qualitatively different. So, system (60) exhibits a supercritical Hopf bifurcation around E 3 considering τ as bifurcation parameter (see Figures 21-23). Also, a rigorous study of the stability and bifurcation of interior equilibrium point E * has been performed. Our analysis describes that the delay within a certain specified range could maintain the stable behaviour of E * . On the other hand, the delay could drive the system into an unstable state. Hence, the study of the timedelay parameter has a regulatory impact on the whole system.

Data Availability
e data used to support the findings of the study are available within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.