Treatment of infected predators under the influence of fear-induced refuge

In this research, we delve into the dynamics of an infected predator–prey system in the presence of fear and refuge, presenting a novel inclusion of treatment for infected individuals in this type of model. Through our analytical efforts, we establish a significant reproduction number that holds a pivotal role in determining disease extinction or persistence within the system. A noteworthy threshold value for this reproduction number delineates a boundary below which the infected population cannot endure in the system. It’s important to note that a range of reproduction numbers leads to both disease-free and endemic scenarios, yet the stability of these situations is contingent upon the initial population sizes. Furthermore, our investigation extends to the exploration of various types of bifurcation-namely, Backward, Saddle-node, and Hopf bifurcations. These findings unravel the intricate and diverse dynamics of the system. Of particular significance is the derivation of an optimal control policy for treatment, augmenting the practical utility of our work. The robustness of our analytical findings is fortified through meticulous verification via numerical simulations. These simulations not only bolster the credibility of our analytical results but also enhance their accessibility. Our study unveils that fear, refuge, and treatment possess individual capabilities to eradicate the disease from the system. Notably, increasing levels of fear and refuge exert a passive influence on the elimination of the infected population, whereas treatment wields an active influence-a crucial insight that bolsters the foundation of our model. Furthermore, our investigation uncovers a spectrum of system dynamics including bistability, one-period, two-period, and multi-period/chaotic behavior. These discoveries contribute to a profound enrichment of the system’s dynamic landscape.


Treatment of infected predators under the influence of fear-induced refuge
Bapin Mondal 1 , Abhijit Sarkar 2 & Nazmul Sk 3* In this research, we delve into the dynamics of an infected predator-prey system in the presence of fear and refuge, presenting a novel inclusion of treatment for infected individuals in this type of model.Through our analytical efforts, we establish a significant reproduction number that holds a pivotal role in determining disease extinction or persistence within the system.A noteworthy threshold value for this reproduction number delineates a boundary below which the infected population cannot endure in the system.It's important to note that a range of reproduction numbers leads to both disease-free and endemic scenarios, yet the stability of these situations is contingent upon the initial population sizes.Furthermore, our investigation extends to the exploration of various types of bifurcation-namely, Backward, Saddle-node, and Hopf bifurcations.These findings unravel the intricate and diverse dynamics of the system.Of particular significance is the derivation of an optimal control policy for treatment, augmenting the practical utility of our work.The robustness of our analytical findings is fortified through meticulous verification via numerical simulations.These simulations not only bolster the credibility of our analytical results but also enhance their accessibility.Our study unveils that fear, refuge, and treatment possess individual capabilities to eradicate the disease from the system.Notably, increasing levels of fear and refuge exert a passive influence on the elimination of the infected population, whereas treatment wields an active influence-a crucial insight that bolsters the foundation of our model.Furthermore, our investigation uncovers a spectrum of system dynamics including bistability, one-period, two-period, and multi-period/chaotic behavior.These discoveries contribute to a profound enrichment of the system's dynamic landscape.
The study of population interactions through mathematical models has become an integral part of ecology.Over the last few decades, mathematical models have been extensively studied to explore the rich dynamics they offer.Among the most important mathematical population models, predator-prey models play a crucial role in understanding the dynamics of different species.These models are developed to capture volatile natural phenomena and have been extensively studied and widely investigated [1][2][3][4][5] .In nature, it's common for predators to hunt prey.However, predators cannot capture all prey because the latter can exhibit anti-predator behaviors like seeking refuge, allowing them to evade predation [6][7][8] .This indicates that refuge behavior protects and prevents the prey from extinction, thus playing a pivotal role in the interaction between prey and predator populations.Recent research on prey refuges in prey-predator models has captured the attention of mathematicians, leading to the discovery of interesting dynamics [9][10][11][12] .Moreover, most predator-prey models assume direct predation as the only way prey are killed.However, this does not align with the real scenario, where prey alter their behavior in the presence of predators [13][14][15] .Recent studies have demonstrated that prey exhibit anti-predator behavior due to predation risks, including habitat changes, vigilance, foraging adjustments, and internal changes.While these anti-predator behaviors enhance short-term survival, they significantly impact reproduction costs in the long term 2 .Furthermore, frightened prey tend to forage less, leading to survival mechanisms like starvation, subsequently affecting birth and death rates [16][17][18] .Whenever it is not always possible to experiment, recently, mathematical modeling has gained popularity as an effective tool for interpreting and studying infectious disease transmission.Mathematical modeling can offer insights relevant to such circumstances.This area which deals with ecological situations with an epidemic, is known as eco-epidemiology in mathematical biology.The effect of disease on the ecological system is a very important issue.Recently, mathematicians have been exploring ecological dynamics from the perspective of epidemiology.A SIRS system was first described by Kermack and Mckendrick 19 , in which the disease becomes transmissible through contact.First and foremost, Anderson and May 20 looked at disease factors in prey-predator model.Most of the previous reviews discussed the dynamics of the predator-prey system characterized by disease in the prey population.These are cited, such as [20][21][22][23][24][25][26][27][28] .These studies were conducted with the aim of exterminating the infection in prey through adequate predation methods or other ecological effects.Thus, ecological models with infected predators are of the utmost importance when addressing the issue of predator conservation.Few cases have been reported where disease spread among predators, rather than among prey.
In 2002, Venturino 29 investigated a predator-prey model involving disease within the predator population.Subsequently, in 2009, Haque and Pal 30,31 examined a prey-predator model where an SIS parasitic infection exclusively affected the predator species.Their findings revealed a unique insight: despite the basic reproduction number for the prey species allowing invasion of the predator-only equilibrium, the presence of infection in the predator population could prevent its extinction.In 2014, Pal and Chakrabarti 31 proposed a predator-prey model with diseases exclusively in predators.Their study demonstrated that, under specific predation levels, both prey and predator species could be rescued from extinction, as long as the disease did not spread within the predator population.Building upon this line of inquiry, Juneja et al. 32 introduced an eco-epidemiological model with delays, assuming infection's impact solely on predators.Their work factored in differential predation rates for infected and susceptible predators, along with distinct harvesting rates.Further developments emerged in 2019, when Huang et al. 33 formulated a delayed SIS model, considering disease transmission only among predator species, while accounting for varied predator responses.Kumar et al. 34 responded to a robust Allee effect in prey populations in 2020 by presenting an eco-epidemiological model.This model assumed disease propagation within predators exclusively and incorporated a ratio-dependent response function.Recently, in 2022, Dutta and Paul 35 crafted a prey-predator model involving disease, fear factors, and additional food for efficient predators.Their model adopted a disease incidence rate inspired by the Beddington-DeAngelis type.Furthermore, Zhang et al. 36 explored an eco-epidemic model with predator-borne diseases influenced by environmental noise.Despite these advancements, the field of mathematical eco-epidemiology has not extensively explored this topic.Consequently, mathematical epidemiology should dedicate more attention to the intricate dynamics of predator-prey interactions, especially when infections affect the predator population.
Defining the disease transmission term is a crucial aspect of epidemic modeling.In many epidemic models, disease spread is assumed to follow the law of mass action.If we denote the densities of susceptible and infectious populations at time t as S(t) and I(t) respectively, then, adhering to the mass action law (also known as bilinear law), the rate of new infections (or incidence rate) at any time t is characterized by β(t) = g(I)S , where g(I) = I .The coefficient is termed the disease transmission coefficient.However, this mass action law presents certain unrealistic characteristics.Notably, the function g(I) becomes unbounded as I grows larger 37 .To address this limitation, Liu et al. 38 proposed a nonlinear incidence rate and introduced a saturated nonlinear function for g(I), specifically g(I) = I p 1+bI q , where p, q are positive constants, and b is a nonnegative constant.In this formulation, I p represents the force of infection, while 1 1+bI q accounts for inhibitory effects stemming from behavioral changes in infectious individuals as their numbers increase or due to the crowding impact caused by the presence of more infected individuals.Consequently, the incidence rate can be expressed as β(t) = I p S 1+bI q .In their study of the cholera epidemic, Capasso and Serio 39 considered an incidence rate of β(t) = SI 1+bI assuming p = q = 1 .It's important to note that the crowding effect is negligible when I is small, and in such cases, the two infection rates become equivalent.Subsequently, various other researchers 38,40,41 have adopted this incidence rate to analyze the dynamics of diverse epidemic models.
According to our knowledge, there are very few published studies on the effects of the treatment function on an eco-epidemiological model.However, none yet consider treatment for an eco-epidemic model where predators get infected.So, this gap needs to explore.Also, there are few eco-epidemic models where the effects of fear and refuge were considered but in those cases, infection was considered among prey species only 28,42 .But, none yet studied these effects by considering disease among predators.Here, is also a gap which need to explore.To our knowledge however, effects of fear and refuge in an infected predator-prey model with treatment has not yet been considered.Motivated by this discussion, in the present paper, we build a predator-prey model considering disease among the predator population only and due to fear of predator prey can take refuge.Additionally, we provide medical resources for treatment of infected individuals.The novelty of our work consists in the study of combined effects of fear and refuge in a prey-predator model where disease can spread among predators along with treatment.This is new model and none yet studied such type of model.This is our main contribution.The following ecological issues, we would like to explore in detail to describe the impact of our work: 1. What are the impacts of fear and refuge on system's dynamics, specifically on infected species? 2. Whether treatment can control the disease or not? 3.Under which conditions disease can be eradicated from the system? 4. Is there any benefit of optimal control policy? 5. Are simultaneously disease-free and endemic steady states possible?If so, under which conditions do they exchange their stability?6.When the ecosystem is in a coexistence state, what types of dynamics occur?
Though it is enough difficult to analyze or discuss all these questions because of it's complexity, but, this investigation can result in a huge eco-epidemiological gain, and we will try to answer these questions by 1. obtaining system's equilibria and complete stability analysis, Our work is organized as follows: we developed an eco-epidemic model in "Model formulation" section by considering the fear effect and prey refuge.The treatment of infected populations was also considered.In "Preliminaries results" section, we have presented some preliminary results related to the mathematical analysis, including the positivity and boundedness of the solutions, the local stability of the system's equilibria, and find the basic reproduction number.In "Extinction and persistence of the disease" section, we discussed the extinction and persistence of the disease such as Backward bifurcation and global stability of DFE (disease-free equilibrium), and the global stability of EE (Endemic equilibrium) points are analyzed.In "Bifurcation results, saddle-node and Hopf bifurcations of the system were discussed.The optimal control policy of the system is discussed in "Optimal control" section.In "Numerical analysis" section, all theoretical results are validated using numerical simulations.Finally, "Discussion and Conclusion" section concludes the paper with some discussion and conclusion.

Model formulation
In order to formulate the mathematical model, assumptions are outlined below.A1: In the considered predator-prey ecosystem, X and Y denote the total prey and predator populations.A2: It has been assumed that predators are susceptible to some transmissible diseases.Disease in the predator population divides the whole predator population into two categories: (i) susceptible predators ( Y s ) and (ii) infected predators ( Y i ).A3: We have assumed that when predators are absent, prey population grows logistically.A4: Prey feels predation fear in the presence of predators.As a result, a linear portion of prey population takes refuge.A5: We have considered that disease spreads among predator species only and the disease is not genetically inherited.The susceptible predator becomes infected under the attack of many viruses or parasites, respectively.Assume that the contact process follows a nonlinear interaction given by Y s Y i 1+bY i .There are compelling justifications for favoring a nonlinear incidence rate over the bilinear formulation 43 .The mass action law lacks saturation and exhibits unbounded behavior as the infectious population, denoted as Y i , increases significantly.Moreover, it neglects the inhibitory impact stemming from the behavioral changes in susceptible individuals and the crowding effect caused by the presence of infectious individuals 39 .A6: Further, we assume that the infected predator population does not have the ability to reproduce as they are vulnerable to the disease.Also, we assume that infected predators can recover from the disease by their self-immunity.A7: It is assumed that both predator populations predate the prey population following Holling type II functional form.But, the impact of disease ( e : 0 < e < 1 ) acts on the predation of infected predator.A8: Furthermore, we provide medical resources for treatment to the infected population.Taking into account the limitations of treatment capacity in different regions or communities, Wang 44 introduced a staged treatment function that better represents real-world scenarios.This function states that the rate of treatment is proportional to the number of infectious individuals when the treatment capacity is not yet reached, and once the capacity is reached, the treatment rate saturates at its maximum level.This approach is more realistic compared to the conventional linear function.Additionally, it's important to consider that treatment effectiveness can be compromised if infected individuals experience delays in receiving treatment.
To capture the saturation effect in treatment, we adopt the treatment function Keeping the above assumption in mind, we draw a schematic diagram, Fig. 1, which clearly shows all the interaction terms among the species.Based on this diagram, we formulate the following prey-predator model where predators get infection: Here, obviously d 4 ≥ d 3 .Moreover, we assume > d 4 , a biologically consistent phenomenon.The ecologi- cal significance of the model parameters is provided in Table 1.Moreover, to analyze system (1), the following initial conditions must be met: (1) (2) X(0) ≥ 0, Y s (0) ≥ 0, Y i (0) ≥ 0.

Preliminaries results
Theorem 1 Solutions of the system (1) with initial condition (2) are non-negative for t ≥ 0.
Proof Let (X(t), Y s (t), Y i (t)) be a solution of the system (1) with initial condition (2).We can derived from the first and third equations of system (1) that In order to demonstrate the positivity of Y s (t) on the interval [0, ∞) , we consider the case as follows.
We first assume that there exists a t 1 > 0 such that Y s (t 1 ) = 0, Ẏs (t 1 ) < 0 and Y s (t) > 0 for t ∈ [0, t 1 ) .There- fore, from the second equation of system (1), we get We get contradiction to the fact that Ẏs (t 1 ) < 0 .Thus, our assumption is not true.Hence, Y s (t)≥ 0 ∀ t ≥ 0 .Therefore, solutions of the system (1) with non-negative initial state are always non-negative.
Theorem 2 All solutions of system (1) initiating in R 3 + are uniformly bounded.
, differentiating Z with respect to t along the solutions of system (1), we have Therefore, . Thus, every solution of system (1) initiating in R 3 + are contained in the region

The predator-free equilibrium
> 0 and Y s is (are) the positive solution(s) of the following quadratic equation where For the parameters value mentioned in Table 1, we can easily obtain the value of X = 0.1868 and for Y s , we plot the polynomial in Fig. 2a.From the figure, we find the value of Y s = 2.18 .Thus, for Table 1, we get a unique disease-free equilibrium E = (0.1868, 2.18, 0).

The endemic equilibrium
X and Y i are the positive solution(s) of the following isoclines Mathematically, it is difficult to analyze the behavior of isoclines ( 4) and ( 5).Hence, we draw these isoclines in X − Y i plane to ensure the existence of endemic equilibrium (see Fig. 2b).It is depicted in the figure that there exist two intersecting points of the isoclines ( 4) and ( 5) for the parameters value mentioned in Table 1.These intersecting points give the endemic equilibrium points E = (0.277, 1.785, 0.601) and E = (0.577, 1.259, 1.452).

Reproduction number
The reproduction number, often denoted as R 0 , is a fundamental epidemiological concept employed to quantify the potential transmission of an infectious disease.It signifies the average count of new infections initiated by a single infected individual within a wholly susceptible population.To illustrate this concept with numerical values, following cases arise: Case-I: R 0 = 1 : This signifies that, on average, every infected individual will pass the virus to exactly one other person.In this scenario, the number of infections remains constant, indicating that the disease is endemic but not spreading rapidly.
Case-II: R 0 < 1 : Each infected person, on average, will transmit the virus to less than one other person.Consequently, the number of infected individuals diminishes over time, leading to the eventual decline and extinction of the disease.
Case-III: R 0 > 1 : In this situation, each infected individual, on average, will transmit the virus to more than one other person.As a result, the number of infected individuals increases, resulting in the rapid spread of the disease within the population.A higher R 0 value indicates a faster rate of disease propagation.
It is crucial to regulate the reproduction number to prevent exponential growth in infections during an outbreak (i.e., maintain it below 1 to bring the outbreak under control) 45,46 .The basic reproduction number ( R 0 ) can be computed using the next-generation matrix formula.In this context, the following theorem is pertinent.

Theorem 3 System (1) has the following basic reproduction number
Proof There is only one infected population ( Y i ) in system (1).In the system, third equation represents the dynamics of infected population.Infection term and remaining terms of that equation, respectively given by Now, at the disease-free equilibrium E , the next generation matrix calculated as where

Local stability analysis
The Jacobian matrix of system (1) is given by J = [a ij ] 3×3 , where 2. The eigenvalues of the Jacobian matrix J( E) are obtained as − d 3 and ) , otherwise it is unstable.3. Calculating at the Jacobian matrix at disease-free equilibrium E , we obtain one eigenvalue as ) , whereas other two are the roots of the following quadratic equation: with Therefore, E will be stable if Otherwise it is unstable.4. The Jacobian matrix at the equilibrium point E is given by J E = [b ij ] 3×3 , where b ij = a ij at E .The characteristic equation of J E is given by where Using Routh-Hurwitz criteria, roots of ( 6) are negative or negative real part if B 2 > 0 , B 0 > 0 and B 1 B 2 − B 0 > 0 .Thus, for these conditions, E will be locally asymptotically stable, otherwise unstable.

Theorem 4
The system (1) undergoes backward bifurcation at R 0 = 1 whenever the sign of the coefficient a ′ is positive where a ′ is defined in (8).
Proof We redefine the system (1) by changing variables as X = x 1 , Y s = x 2 and Y i = x 3 .Applying the vector notation x = (x 1 , x 2 , x 3 ) T , the system (1) can be written as dx dt = f (x) where f = (f 1 , f 2 , f 3 ) T as follows: (6) (7) Vol:.( 1234567890 The disease free-equilibrium of the above system is given by E ′ (x 1 , x 2 , 0) , where x 1 = (1−m) > 0 and x 2 is (are) the positive solution(s) of the following equation Taking γ as a bifurcation parameter we found R 0 (γ = γ * ) = 1 which gives that the Jacobian of the transformed system (7) at the disease free equilibrium at γ = γ * has a simple zero eigenvalue, while all other eigenvalues have nega- tive real part.The Jacobian matrix of the above system at disease free equilibrium is given by where Now, we use the centre manifold theorem discussed by Castillo-Chavez and Song 47 to study the dynamics of the system (7) near γ = γ * .The necessary computation for the theorem given below: The right eigenvector of the Jacobian matrix J(E ′ ) associated with zero eigenvalue at γ = γ * is given by and w 3 = 1 .It is also easy to see that the left eigenvector of the Jacobian matrix J(E ′ ) associated with zero eigenvalue at γ = γ * is given by V where, v 1 = 0, v 2 = 0 and v 3 = 1 satisfying the condition V .W = 1.
The bifurcation coefficient a ′ and b ′ at the disease free equilibrium and at γ = γ * are given by Hence, following Theorem 4 in 48 , the system (1) undergoes backward bifurcation at R 0 = 1 whenever a ′ > 0 .
Theorem 5 For R 0 < 1 , the disease-free equilibrium E is globally asymptotically stable in the interior of the positive Proof The global asymptotically stable of the boundary equilibrium point E(X, Y s , 0) is shown with the help of Lyapunov function.Let us consider a Lyapunov function V as, where σ is a positive constant which will be chosen later.Taking trajectory derivative of this scalar valued func- tion along the solutions of the disease-free subsystem of the system (1), we get , which yields M = {E} .Note that M is invariant.By applying LaSalle's invariance principle for delay differential equations 49 , we establish the global asymptotic stability of the disease-free equilibrium E under the specific conditions delineated in the theorem.
Theorem 6 For R 0 > 1 , the endemic equilibrium E is globally asymptotically stable if following conditions are satisfied Proof In order to prove this theorem, we consider a Lyapunov function V that is positive definite Now calculating the time derivative of V along the solution of the system (1), we get where T and V * is a symmetric matrix given by with Note that under given conditions specified in the theorem, V * is positive definite matrix, then dV dt will be negative definite.Hence, for R 0 > 1 , E is globally asymptotically stable according to the Lyapunov stability theo- rem.
Remark From a biological perspective, the global stability of disease-free equilibrium points and endemic equilibrium points are of the utmost importance.As a result of the analysis, we obtain conditions under which diseases Y s ](X − X) 2 < 0 (by the condition given in the theorem).
, www.nature.com/scientificreports/may or may not survive in the system.Moreover, in these situations, regardless of population size, either disease will persist or not.

Saddle-node bifurcation analysis
System (1) have at most two endemic equilibrium points, depending on m.Assume that the two endemic equilibrium points coincide at m = m [SN] .On one side of m = m [SN] , there are two endemic equilibrium points, while on the other there are no endemic equilibrium points.Consequently, the saddle-node bifurcation may occur around the endemic equilibrium point when the system parameter m crosses its bifurcation threshold.In the following theorem, we shall establish the conditions for the occurrence of saddle-node bifurcation concerning m as the bifurcation parameter.
Theorem 7 System (1) undergoes saddle-node bifurcation at the coincident endemic equilibrium E when it goes through the critical value m = m [SN] and satisfy following conditions Proof We must use Sotomayor's theorem in order to establish dynamics of the system in the neighbourhood of the equilibrium point E since the eigenvalue of J( E) is zero at m = m [SN] .The following function is considered for this purpose: and F 1 , F 2 and F 3 represent right hand sides of first, second and third equations of system (1), respectively.Then, we have and T are the eigenvectors corresponding to the eigenvalue zero of the matrix DF(Z = E, m = m [SN] ) and its transpose, then we have Now, Thus,

Now,
In other words, if the conditions stated in the theorem are met, system (1) experiences saddle-node bifurcation around the coincident endemic equilibrium point E .

Theorem 8 System (1) possesses Hopf bifurcation if
when m crosses the critical value m = m [HB] .
Proof Hopf bifurcation in our model arises when one of the eigenvalues of the variation matrix J( E) has nega- tive real part with Re( dξ dm ) m=m [HB] � = 0 and other two eigenvalues are purely imaginary.The characteristic Eq. ( 6) must satisfy the conditions B 2 B 1 − B 0 = 0 for purely imaginary eigenvalues.A purely imaginary eigenvalue is assumed to be iω for m [HB] .
Differentiating the characteristic Eq. ( 6) with respect to m, we obtain Now, putting ξ = iω and comparing with real part, we have Thus, system (1) experiences a Hopf bifurcation with respect to the parameter m if condition satisfies at m = m [HB] .
As it is very difficult to follow the Hopf bifurcation result analytically, so, we verify it numerically.We have taken the values of all the parameters from Table 1 except m.Calculating B 1 B 2 − B 0 = 0 for m, we find the critical value of m as m [HB] = 0.23058 .Then at this critical value Eq. ( 6) becomes as follows: Solving this equation, we get one of the eigenvalues as iω (where ω = 0.0673 ), which is purely imaginary.Now, we calculate the expression Re dξ dm at ξ = iω and get 7.36 ( = 0) .Thus, the transversality condition verified at m = m [HB] .Therefore, this numerical example confirms the existence of Hopf bifurcation of system (1) at m = m [HB] .

Optimal control
Controlling the predator-prey system is possible if humans are able to reach a certain limit.To prevent infection on predator, control is applied in this model.The purpose of the optimal control is to optimize the size of susceptible predator and minimize the cost of treatment of infected individuals.A problem of optimal control is considered in which the objective function should be minimized as follows: where A and B are respectively weighted constants per capita less due to presence of infected population and treatment of infected individuals.Our problem is to find optimal control function (γ * (t)) such that The Lagrangian of the given problem (1) is to defined as L = AY i + Bγ 2 .Now, we form the Hamiltonian H for the problem is given by Re dξ dm Vol:.( 1234567890 www.nature.com/scientificreports/Using Pontryagin's Maximum Principle, the adjoint equations are given by with the transversality conditions i (T) = 0, i = 1, 2, 3. Thus, we have with the transversality conditions Now, using the optimality condition The value of γ is in the interval [0, 1] so that some possibilities are obtained below Thus, the following optimal control variable is obtained ) is solution of the system (12) with the condition (13).In the following theorem, we summarize the details.

Theorem 9
The optimal control γ * which minimize the objective function J over the region U shown in (11) are given by γ * = min max 0,

Numerical results
To explore the insight dynamics of the system (1) and verify analytical results, we perform extensive numerical simulations in this section.Simulations are performed in MATLAB and numerical continuation software MAPLE.In our numerical simulations, we have selected a collection of biologically plausible hypothetical parameter values, which are outlined in Table 1.Although these values are hypothetical in nature, they are drawn from existing literature sources 27,28,41,42,51,52 .Notably, these values closely resemble those found in these existing works.
In order to determine the influence of parameters on R 0 , we compute normalized forward sensitivity indices for R 0 53 .An index of normalized forward sensitivity measures how much the variable has changed based on changes in the parameter.In Fig. 3, a plot is displayed showing the normalized forward sensitivity indices for each parameter in the expression of R 0 .This figure illustrates that the value of R 0 increases as the parameters , 1 .
Vol.:(0123456789) www.nature.com/scientificreports/and β have positive indices.On the other hand, parameters µ , γ and d 4 possess negative indices with R 0 .Figure 3, reveals that the parameters , β and γ are most sensitive to R 0 .It is important to note that we prefer low values of R 0 as they enhance the possibility of eradicating the disease.The increment of parameters and β should therefore be prevented at all costs.On the other hand, the increment of µ , γ , and d 4 should be encouraged.Therefore, in order to eradicate disease, policymakers should focus on preventing an increase in the parameters with positive indices, while stimulating the parameters with negative indices.The justification for the Backward bifurcation can be visually grasped through a one-parameter bifurcation diagram utilizing R 0 as a basis.As evident from Theorem 4, the exploration of backward bifurcations is contingent upon the parameter γ .This relationship is represented graphically in terms of R 0 in Fig. 4. The diagram clearly indicates that the disease-free equilibrium point is stable if R 0 < 1 and unstable if R 0 > 1 .Conversely, for R 0 > 1 , only one endemic equilibrium ( E ) exists and is stable.Furthermore, there exists a range of R * 0 < R 0 < 1 (with  www.nature.com/scientificreports/R * 0 < 1 ), where a significant phenomenon arises.In this interval, the system exhibits two interior (endemic) equilibrium points.Among these points, the one characterized by a higher infected density is stable, while the other is unstable.This scenario leads to a state of bistability, where both the DFE and the EE stable simultaneously.The second row of Fig. 4 visually illustrates these scenarios through phase portraits.In the first figure (where R 0 < R * 0 ), all trajectories converge to the disease-free state; in the second figure (where R 0 ∈ [R * 0 , 1] ), the trajectories' convergence depends on initial population size, with some reaching the disease-free state and others the endemic state; in the last figure (where R 0 > 1 ), all trajectories converge to the endemic state.Further, insights can be gleaned by considering the specific ranges of R 0 .For R 0 < R * 0 , only the DFE exists and is globally asymptotically stable, as depicted in Fig. 5a.Conversely, for R 0 > 1 , one endemic equilibrium ( E ) exists and is globally asymptotically stable, as evident from Fig. 5b.The Backward bifurcation thus plays a pivotal role in comprehending the dynamics of disease within an ecosystem.Now, we see the effect of important model parameters on the system's dynamics by varying those parameters in an interval of values.It is clear from Fig. 6a that for lower values of r, the system exhibits a disease-free stable state up to a certain value after that system oscillates around the disease-free state.Note that in this case fluctuation level is very high.After that disease affected the system and the system shows oscillatory behavior around the endemic state.Here the level of fluctuation is very low.Further higher values of r give a stable endemic state.It is apparent from Fig. 6b, that opposite dynamics could be observed for the parameter K with respect to the parameter r.That is as the value of K changing system exhibits a stable endemic state to an oscillating disease-free state to a stable disease-free state.The refuge parameter stabilizes the endemic state from the oscillatory state, Fig. 6c.An oscillatory disease-free state with higher amplitude is observed for lower values of and moderate values of , system exhibits an endemic state with oscillating behavior of lower amplitude.Finally, the system falls into the stable endemic state for higher values of , Fig. 6d.Clearly, a totally opposite dynamics can be seen for the parameters µ and γ to the parameter (see Fig. 6e,f).Therefore, one can easily conclude from these bifurcation results that by manipulating these parameters, we can easily remove the disease from the system.
Next, to see the effect of fear and refuge on the system's dynamics, we draw two parametric bifurcation diagram by varying these parameters (m and K) at once.The whole region is divided into 6 subregions ( R 1 − R 6 ) with different dynamics, Fig. 7.It is clear from the figure that there is a large region R 1 in which a disease-free state is stable (see Fig. 8R 1 ), which is biologically most preferable region.In region R 2 , disease-free equilibrium oscillates, Fig. 8R 2 .Note that in these two regions diseases can not persist in the system.It indicates that higher values of both fear and refuge give a disease-free state.Next, as the system enters into region R 3 , it shows bista- bility, i.e., here our system shows oscillatory behavior around both endemic and disease-free equilibrium (see Fig. 8R 3 ).That is in R 3 whether disease will persist or not persist, depends on the initial population size.Further, for lower values of K, the system enters into the region R 4 .In this case, also bistability occurs.But, here endemic equilibrium is stable and disease-free equilibrium behaves periodically, Fig. 9R 4 .It is apparent from Fig. 9R 5 that in region R 5 both disease-free and endemic states are stable.Finally, in region R 6 , interesting dynamics can be observed.Here, up to a certain time system behaves in an oscillatory manner around the disease-free state and ultimately goes to the endemic stable state (see Fig. 9R 6 ).Furthermore, all the characteristics of each subregions ((R 1 − R 6 )) are explicitly given in Table 2.One can easily follow this table by corresponding region's figures.
Moreover, for a different value of µ other than the value mentioned in Table 1, we plot the bifurcation diagram with respect to the treatment parameter γ , Fig. 10.All the bifurcating results due to the variance of γ are depicted in phase portraits plotted in Fig. 11.It is transparent from these phase portraits that for lower values of γ , the system shows a stable endemic state.As the value of γ increases, the system exhibits one period to two periods to chaos to one period to a stable state of endemic equilibrium.Finally, for higher values of γ , the system settles to a periodic disease-free state.
Furthermore, to be more transparent about the change of density of the infected population by varying model parameters , K, m and γ see Fig. 12.It is apparent from the figures that is directly proportional to the density of the infected population whereas other parameters (K, m and γ ) are inversely proportional to infected popu- lation density.It means if the value of increases (decreases), the density of the infected population increases (decreases) accordingly.On the other hand, if the value of each of these parameters K, m, and γ increases (decreases), the density of the infected population decreases (increases) accordingly.Thus, from these results also, one can conclude that the infection can be suppressed by fear, refuge, and treatment.From the biological point of view, these results are very important to conserving any ecosystem.
In order to solve the optimality system numerically, the forward-backward sweep method 54 has been used.Some initial guesses about the control variable are made before the process begins.Runge-Kutta fourth order iterative scheme is used first to solve the system of state variables forward in time, then to solve the system of    1. adjoint variables backward in time.In each iteration, the control value is updated.Until the desired convergence is achieved, the process is repeated.Weight factors associated with infected predator and treatment are chosen as A = 0.01 and B = 0.05 , whereas the other parameter values are same as in Table 1 with the initial population sizes X(0) = 4 , Y s = 1 and Y i = 2 .Moreover, we consider the time interval as [0, 35].The simulations for the density of prey, susceptible predator, and most important infected predator populations under the optimal control parameter γ (t) and without control strategies are shown in Fig. 13.The drop in the infected predator density and rise in the susceptible predator density in presence of time-dependent optimal control is apparent in the figure .A plot of γ (t) 's profile is shown in Fig. 13b,c along with the transversality condition i (t) = 0 .In this figure, it can be seen that it is optimal to exert efforts until 16.7 units of time, then lower them down.As time increases, i converges to zero as we observe in the adjoint variables.From these figures, one can say that optimal control in treatment is very effective to reduce the number of infected predators and eradicate the infection from the system.

Discussion and conclusion
This research proposes a novel predator-prey model incorporating the effects of fear and fear-induced refuge, where the predator population can become infected by pathogens or viruses.The study introduces a treatment strategy for the infected predator individuals.To the best of our knowledge, this type of predator-prey model involving treatment is innovative and previously unexplored.The primary objective of this study is to investigate the impacts of fear, fear-induced refuge, and treatment on an infected predator-prey system.To validate the model, we first analyze the positivity and boundedness of the system's solutions.Subsequently, we examine the biologically feasible equilibria of the system.Interestingly, we derive a reproduction number that plays a crucial role in discussing the possibilities of disease elimination and persistence within the system.This analysis involves rigorous exploration of Backward bifurcation, as well as local and global stability.Additionally, we delve into the dynamics of the system by studying Saddle-node and Hopf bifurcations.In addition to exploring the system's dynamic behavior, we also delve into the dynamics of the infected population.Here, our goal is to mitigate the density of the infected population by exerting control over infected individuals and the medical  1.
Table 2. Equilibrium points and their characteristics in each region as illustrated in Fig. 7.

Region
Equilibrium point Characteristic of equilibrium point Saddle, saddle, unstable spiral (stable limit cycle) Saddle, saddle, unstable spiral (stable limit cycle), unstable, unstable spiral (stable limit cycle) Saddle, saddle, unstable spiral (stable limit cycle), unstable, stable resources allocated for treatment.To this end, we establish the framework of an optimal control problem and leverage Pontryagin's maximum principle to arrive at a solution.Optimal control emerges as significantly more potent in diminishing the density of the infected population while concurrently augmenting the density of susceptible predators.
Our findings suggest that when the basic reproduction number, denoted as R 0 , is less than a critical threshold R * 0 , the disease is unable to sustain itself within the system.However, within the range R 0 ∈ [R * 0 , 1] , the persistence of the disease becomes contingent upon the initial population sizes.Conversely, once R 0 surpasses unity, the  7 are presented for the system (1).Remaining parameters adopted from Table 1.The complete phase portrait along with the corresponding time series for the various regions outlined in Fig. 7 are presented for the system (1).Remaining parameters adopted from Table 1.
disease becomes a permanent fixture within the system.Furthermore, our bifurcation analyses indicate that the individual growth rate of prey and the disease prevalence rate independently contribute to the escalation of the infected population's density.In contrast, factors such as fear, fear-induced refuge, recovery rate, and treatment each play a distinct role in diminishing the density of the infected population.Remarkably, higher values of these influencing factors ultimately lead to the eradication of infected individuals from the system.Moreover, our findings indicate that when fear is at a very low level, the infected population persists within the system for all future time intervals.Increasing the level of fear introduces a nuanced dynamic: the persistence or extinction of the disease depends on the initial population size.Nevertheless, it's important to note that the infected population cannot sustain itself within the system when fear and refuge are elevated to higher levels.Furthermore, through manipulation of the recovery rate, intriguing dynamics emerge.As the treatment intensity increases, the system transitions from a single periodic behavior to two periodic cycles, further evolving into multiple periodic or even chaotic states, and then reverts back to a single periodic pattern.Ultimately, the system stabilizes into an endemic state before eventually converging to a disease-free periodic state.This exploration unveils an interesting phenomenon: when medical resources are limited, the disease persists, accompanied by chaotic behavior.Conversely, an ample supply of medical resources can effectively manage the chaos and lead to disease eradication.The overarching conclusion drawn from this investigation is that the study's findings hold significant real-world implications.Moreover, it sheds fresh light on the intricate dynamics inherent in a predator-prey system afflicted by infections within the predator population.
The proposed model finds relevance in diverse biological scenarios where predator-prey dynamics and disease interactions are influenced by ecological factors such as prey refuge and fear-induced behaviors.The model can be applied in fields like wildlife conservation, agricultural systems, epidemiology, disease ecology, ecological restoration, and behavioral ecology.The adaptability of the model to simulate different scenarios underscores its capacity to address a wide range of ecological contexts where predator-prey dynamics and disease interactions play a pivotal role.
In wildlife conservation setting, where the predators could represent natural predators like wolves or big cats, and the prey could be animals like ungulates or small mammals.The inclusion of disease dynamics and ecological effects like refuge and fear could provide insights into disease spread within wildlife populations and how these interactions impact conservation efforts.Agricultural Systems: In agricultural settings, the predator-prey dynamics could be relevant to pests and their natural predators.Understanding how disease transmission and fear-induced behaviors impact the dynamics of pest control can have implications for sustainable agriculture practices.1.Moreover, we take A = 0.01 and B = 0.05. https://doi.org/10.1038/s41598-023-43021-0www.nature.com/scientificreports/

Figure 1 .e 5 bβ. 2 the delay in treatment for the infected individuals d 4
Figure 1.A complete schematic diagram describing growth, death and the dynamical interactions among prey (X), susceptible predator ( Y s ) and infected predator ( Y i ) in the context of disease, fear, refuge and treatment.

Figure 2 .
Figure 2. In figures, (a) red curve represents polynomial (3) and blue dot is the root of the polynomial; (b) magenta curve represents isocline (4) and blue curve represents isocline(5).From Table1, the parameters' values are selected.

Figure 4 .
Figure 4. First row: Backward bifurcation with respect to R 0 .Second row: Phase portraits for different values of R 0 .Parameters have the same values as in Table 1 except m = 0.7.

Figure 5 .
Figure 5. Solution trajectories of system (1) in X − Y s − Y i space for different initial values.Parameters are at the same values as in Table1except m = 0.7 .In the figures, green dots represent the equilibrium points.Circles in each trajectory describe the initial start positions of the trajectories.

Figure 6 .
Figure 6.Bifurcation diagrams of system (1) with respect to the parameters (a) r, (b) K, (c) m, (d) , (e) µ and (f) γ .Rest of the parameters are at the same values as in Table1.

Figure 7 .
Figure 7. Two parametric bifurcation diagram with respect to m and K. Remaining parameters have the same values as in Table1.

3 Figure 8 .
Figure 8.The complete phase portrait along with the corresponding time series for the various regions outlined in Fig.7are presented for the system (1).Remaining parameters adopted from Table1.

6 Figure 9 .
Figure 9.The complete phase portrait along with the corresponding time series for the various regions outlined in Fig.7are presented for the system (1).Remaining parameters adopted from Table1.

Figure 10 .Figure 11 .
Figure10.Bifurcation diagrams of system (1) with respect to the treatment parameter γ .Rest of the parameters are at the same values as in Table1except µ = 0.075.

Figure 12 .
Figure 12.Variations in the infected predator populations for different values of (a) , (b) K, (c) m and (d) γ .Parameters are at the same values as in Table 1 except in (b) γ = 0.51 (c) γ = 0.31.

Figure 13 .
Figure 13.Time series solutions of (a) each population with and without control; (b) treatment control ( γ ); (c) adjoint variables ( i , i = 1, 2, 3 ).Parameters are at the same values as in Table1.Moreover, we take A = 0.01 and B = 0.05.