An ecoepidemic model with healthy prey herding and infected prey drifting away

. We introduce here a predator–prey model where the prey are affected by a disease. The prey are assumed to gather in herds, while the predators are loose and act on an individualistic basis. Therefore their hunting affects mainly the prey individuals occupying the outermost positions in the herd, which is modeled via a square root functional response. The conditions of boundedness and uniform persistence are established. Stability and bifurcation analysis of all feasible equilibrium are carried out. Conditions on the model parameters for the possible existence of limit cycles are derived, global stability analysis is also shown in proper choice of suitable Lyapunov function. Numerical simulation of the various bifurcations validate the theoretical results. It is found that the system ultimate behavior depends mainly on two crucial parameters, the force of infection and predator average handling time. A discussion of the biological signiﬁcance of the investigation concludes the paper.


Introduction
In recent times, herd behavior was introduced in interacting population systems of various nature, of symbiotic, competing and predator-prey type, by [2], via a square root functional response, which models the predator-prey interaction occurring mainly on the perimeter of the prey herd.The main result appears to be the possible onset of persistent oscillations, which are not originated by the Holling type II functional response.Thus a different mechanism, other than feeding saturation, entails as consequence the insurgence of limit cycles.It should be noted that in fact, the idea of a nonlinear response involving a power law for the prey dates back a few decades, appearing indeed already in the original work of Gause [8].Anyway, after the appearance of [2], other researchers undertook similar investigations based on the same idea.In particular, [4] considers a simple predator-prey model with interactions on the boundary, but in which also feeding saturation through a Holling type II functional response appears.Note also that the shape of the herd does not remain circular at all times so that, in general, the square root term becomes inadequate.But [5] shows that the modification to account for this change via a generic exponent in place of the square root do not lead to substantial modifications in the results.Further, a minimal model combining ecoepidemics and group behavior is introduced in [22], while an ecoepidemic model with prey herd behavior and an infection in the predators appears in [3].Ecoepidemic models with Holling type I and II responses but without the group effect are, for instance, found in [12,19].In [24], they introduce a pathogenic agent, which influences the force of infection.Also, a predator-only equilibrium exists, which not appear in our system.Biologically, this occurs due to the absence of other food sources except for prey.
In [23] the alternative cases of infection remaining harmless for predators, infected individuals being not predated and the infected prey being toxic for predators were considered.It is assumed that susceptible prey behave individually, while infected gather together.Occasional contacts among the solitary susceptibles S and the herd of infected lead to new contagions via interactions with the infected individuals I occupying the outermost positions in the group, this being modeled by the term S √ I. Our basic assumption here is to revisit the assumptions of [23], exchanging the roles of the two prey subpopulations.Thus susceptible prey group together, while infected abandon the herd and behave individually.Still, the susceptibles can be infected via interactions that occur on the herd boundary, modeled in this case via the term I √ S. The paper is organized as follows.In the next section, we formulate the model and then discuss the boundedness and the persistence of the species.Next, we analyze the different kinds of feasible equilibria and carry out their stability analysis as well as investigating the Hopf and other bifurcations.Using Bendixson's criterion, we show the conditions for nonexistence of limit cycles at E * .Applying instead the center manifold theorem, we also investigate the onset and nature of the bifurcating limit cycles.By using the results on uniform persistence, the conditions for global stability are obtained, constructing the suitable Lyapunov function.The fourth section contains the numerical validation of the results as well as the bifurcation diagrams.A final comparison with the earlier model of [23] concludes the paper.

Mathematical model formulation
We consider the situation in which prey are affected by an unrecoverable disease so that their population X is partioned among susceptibles s and infected i, X = s + i. Infected Nonlinear Anal.Model.Control, 28(2):326-364, 2023 are assumed to be too weak to reproduce and compete with susceptibles for resources.They also drift away from the herd of susceptible prey and therefore are hunted on a oneto-one basis by predators.Also, the disease is assumed not to spread to the latter.The predators furthermore hunt the susceptible prey herd by mainly capturing the individuals occupying the perimeter of the herd and are subject to feeding satiation.Thus hunting on susceptible prey obeys a law that is the combination of Holling type II and the square root functions: where T h represents the predator's prey handling time, and α is the predator's search efficiency.Finally, infected prey with occasional contacts with the susceptible's herd could transmit the disease.Again, we assume that this process occurs mainly on the boundary of the herd.With these assumptions, the model reads: The first equation models the susceptible prey dynamics: they reproduce logistically with reproduction rate r and carrying capacity K, are hunted on the perimeter of the herd by predators with the functional response (1) and become infected with transmission rate λ, again by contacts on the herd boundary.
The second equation describes the infected prey recruited by "successful" contacts with susceptibles, hunted at rate m by predators and subject to natural plus disease-related mortality µ.
Finally, the third equation contains the predators evolution.They are specialist on the modeled prey, hunt the susceptible ones by predating on the herd perimeter with conversion rate θ 1 , while capture the infected ones on an individual basis and with possibly a different conversion rate θ 2 due to the fact that the latter might be less palatable or contain different, most likely less, nutrient than the susceptible ones.Their mortality rate is δ.
Note that this model differs from other similar ones already published.In particular, with respect to [16], here infected prey are also abandoned, but we assume that they occasionally can still interact with the susceptible ones on the boundary of the herd, while in [16], such contagion process is assumed to occur inside the herd, before the infected abandon the it.A different mechanism is instead modeled in [6], where infected stay in the herd and are subject to hunting as the susceptibles are, this predation occurring once more on the perimeter of the herd.https://www.journals.vu.lt/nonlinear-analysisLetting √ s = S, then the above system reduces to which is now the subject of our investigation.
3 Mathematical analysis

Persistence
A system is persistent if a compact set Ψ ⊂ G = {(S, I, P ): S > 0, I > 0, P > 0} exists so that all the solution of system (2) eventually enter and remain in it.
Theorem 1. System (2) is persistent if the following conditions hold:

Boundedness
Proposition 1.The healthy prey population is uniformly bounded.Proof.From system (2) we have Now there exist a A 1 > 0 such that for all t > A 1 , we have S(t) < K + = D. Proposition 2. The second equation of (2) can be bounded from above as follows: Therefore when R 0 > 1, then epidemic will spread, and if R 0 1, then epidemics is eradicated.
Theorem 2. All the solutions of system (2), which initiate in R + 3 , are uniformly bounded.

Stability analysis
The system Jacobian matrix J B ≡ DF of system (2) is given by The behavior near E 0 (0, 0, 0) The eigenvalues of the Jacobian matrix J E0 are which shows that E 0 is an unstable hyperbolic critical point, namely, a saddle with the instability in the orthogonal direction to the IP -plane.

Stability
At the equilibrium point E 1 , the eigenvalues of the Jacobian J 1 are Thus E 1 will be asymptotically stable if
and g 2 represent any two nonzero real numbers.Since 2) will experience neither a saddle node nor a transcritical bifurcation at , let γ 3 and µ 3 be the eigenvectors respectively corresponding to the 0 eigenvalues of the matrix J 1 and its transpose.We obtain µ T 3 = (m 3 , 0) so that m 3 , where g 3 , m 3 are two nonzero real numbers.Then 3 , 0 rλ 2µ , 0, 0 3 m so that the system will experience neither a saddle node nor a transcritical bifurcation at , let µ 4 , γ 4 respectively denote the eigenvectors corresponding to the eigenvalue 0 of the matrix J 1 and its transpose.Then we get µ T 4 = (m 4 , 0, m γ T 5 = (0, 0, g 4 ) so that m where g 4 , m [3] 4 are two arbitrary nonzero real numbers.Then 4 , 0 for which the system have neither a saddle node nor a transcritical bifurcation at
We defer the proof to Appendix G.

3.4.3
The behavior near E 2

Stability
One eigenvalue of the Jacobian J 2 is and the other two are the roots of the equation x 2 + a 1 x + a 2 = 0, where E 2 will be locally asymptotically stable if the Routh-Hurwitz conditions are met, a 1 > 0 and a 2 > 0, which explicitly become

Bifurcations
The matrix and the values of ] and µ [te 2I ] are reported in Appendix E.
https://www.journals.vu.lt/nonlinear-analysis (i) For λ = λ [te2] , let µ 5 , γ 5 respectively be the eigenvectors corresponding to the zero eigenvalue of the matrix J 2 and its transpose.We obtain γ T 5 = (0, g and m 5 are any two positive real numbers.Since 5 = 0, the system experiences a transcritical bifurcation at , there will be a transcritical bifurcation around E 2 .
For λ = λ [te 2I ] , let µ 6 , γ 6 respectively be the eigenvectors corresponding to the zero eigenvalue of the matrix J 2 and its transpose.We obtain γ T 6 = (0, 0, g and m 6 are any two positive real numbers.Since the system will experiences neither saddle-node nor a transcritical bifurcation at = 0 and the other eigenvalues of J 2 are nonzero, we will get K = K [te 2I ] so that let µ 7 , γ 7 respectively be the eigenvectors corresponding to the zero eigenvalue of J 2 and J T 2 .Then µ T 7 = (m where m 7 are two nonzero real numbers.Since 7 m [3] 7 7 , 0 The proof is given in Appendix G.

Bifurcations
The matrix When λ , and h https://www.journals.vu.lt/nonlinear-analysis and h [te 3I ] is the positive value of The values of , there exist eigenvectors µ 8 , γ 8 corresponding to the vanishing eigenvalue of J 3 and its transpose, respectively.Then µ T 8 = (m where m 8 are nonzero real numbers.Since 8 , 0 the system will experiences neither saddle-node nor transcritical bifurcations at , there exist eigenvectors µ 9 , γ 9 corresponding to the vanishing eigenvalue of J 3 and its transpose, respectively.Then µ T 9 = (m 9 ), γ T 9 = (0, 0, g where m 9 are nonzero real numbers.Since 9 , 0 the system experiences a transcritical bifurcation at , there exist eigenvectors µ 8 , γ 8 corresponding to the vanishing eigenvalue of J 3 and its transpose, respectively.We find µ T 10 = (m 10 , 0) so that m where m 10 are any nonzero real numbers.Since 8 , 0 (µ 10 , µ 10 ) = 0, the system will experiences neither saddle-node nor transcritical bifurcations at , there exist eigenvectors µ 11 , γ 11 corresponding to the vanishing eigenvalue of J 3 and its transpose, respectively.Then µ T 11 = (m 11 , 0, m 11 , where m 11 are nonzero real numbers.The system will experiences transcritical bifurcation at E 3 for h = h [te 3I ] in view of the following results: 11 m [1] 11 m [3] 11 δ 2 θ 1 αK = 0.

Global stability analysis at E 3
Theorem 5.The equilibrium E 3 (S 3 , 0, P 3 ) is globally asymptotically stable if the following conditions hold: The proof is contained in Appendix G.
The proof is once more written in Appendix D.

Other bifurcations
We rewrite system (2) as where Let U = (γ 1 , γ 2 , γ 2 ) T , V = (µ 1 , µ 2 , µ 3 ) T respectively denote the two eigenvectors corresponding to the zero eigenvalue of J and its transpose J T , then we have with γ 1 , γ 2 , µ 1 , µ 2 being arbitrary nonvanishing real numbers.Since V T F m (E * , m 4 ) = 0 and V T DF m (E * , m 4 )(U ) = 0, there are three cases: Case 1.There exist a Turing-saddle node bifurcation node if the following conditions hold: Case 2. There is a Turing-transcritical bifurcation if conditions (12) hold along with Case 3.There is a Turing-pitchfork bifurcation if conditions ( 12) and ( 13) hold together with the critical value of m and the variational matrices being reported in Appendix E.
3.4.5.4 Global stability at E * Theorem 6.The coexistence equilibrium E * (S * , I * , P * ) is globally asymptotically stable if it satisfies conditions similar to the ones of [12,20] and [19], i.e., The proof is given in Appendix G. To show the nonexistence of periodic orbits of the system near E * (S * , I * , P * ), we apply the technique of Li and Muldowney mentioned in [19].Recalling the notation in (11), we denote by J [2] the N 2 × N 2 second additive compound matrix associated with the Jacobian J [19].Explicitly, We now use the results of [17] on Bendixson's criterion in R n to analyse closed orbits.
Theorem 7. A simple closed rectifiable curve that is invariant with respect to system (2) cannot exist if the following condition holds [12]: Proof.Let us reduce the system dimension by setting https://www.journals.vu.lt/nonlinear-analysisNow (15) becomes Sufficient conditions to satisfy (16) are In these conditions, system (2) admits no periodic solutions.

Stability of bifurcating limit cycle
We now establish the stability of the limit cycle arising from the Hopf bifurcation.To this end, we apply the center manifold theorem [15,20].Since the Jacobian J * (S * , I * , P * ) has purely imaginary eigenvalues leading to the Hopf bifurcation, we can analyze the present system just on a two-dimensional manifold, where the flow is exponentially contracting.
We translate E * to the origin by S = S − S * , Ī = I − I * and P = P − P * .The original system becomes where the Jacobian J * is and Nonlinear Anal.Model.Control, 28(2):326-364, 2023 We neglect the higher-order term such as S3 and the ones containing At the Hopf bifurcation, the characteristic equation has the eigenvalues λ 1 = −a 1 = η for some η, and λ 2,3 = ±ı √ a 2 = ±ıω, where Let the eigenvector of J * associated with λ 1 be η 1 , and the ones corresponding to λ 2,3 be η 2 ± ıη 3 , where η 1 , η 2 , η 3 denote real vectors.Then it can be shown that the matrix is nonsingular and and Next, letting Ȳ = ( S, Ī, P ) T , we use the following linear transformation: Substituting ( 19) into (17), we get Now ( 20) can be rewritten as where l = (I 1 , I 2 ) T , z = (I 3 ).H and N are the constant matrices and F 1 and G 1 are C 2 functions.System (20) can now be written as Now system (22) has a local center manifold z = f (l), l < , where f is in C 2 .The function f (l) can be approximated arbitrarily closely by a Taylor series as shown in the next theorem.
https://www.journals.vu.lt/nonlinear-analysisTheorem 8. Let φ : R n → R m be C 1 in a neighborhood of the origin, φ(0) = 0, Hence in the present case the center manifold up to a quadratic approximation is described by , which leads to Calculating φ 1 , φ 2 , φ 3 of ( 19), we have Now from ( 18) and ( 24) we get which equals the right-hand side of (23).Comparing both sides, the coefficients of I 2 1 , I 1 I 2 , I 2 2 , we get Here Then the flow on the center manifold is governed by the two-dimensional system The central manifold theorem tell us that (25) contains all the information needed to determine the asymptotic behavior of the solution of (21).
In detailed form, (25) can be written as d dt where, letting h.o.t.standing for higher-order terms, https://www.journals.vu.lt/nonlinear-analysis The stability of the limit cycle arising from a Hopf bifurcation is determined by the sign of Ω: If Ω < 0, the Hopf bifurcating limit cycle is stable, and the Hopf bifurcation is supercritical; if Ω > 0, the bifurcating limit cycle is unstable, and the Hopf bifurcation is subcritical.

Numerical settings
For the convenience of the reader, all the conditions found and the numerical data used in the following are listed in Tables 1-6.
Table 1.Conditions of equilibrium global stability.

Equilibrium Global stability condition E 1
Unconditionally globally stable --  The disease-free equilibrium E 1 = (7.1,0, 0) can be achieved with the parameters of Table 2, which satisfy the feasibility condition K < min{100, 59.2}.
The equilibrium E 2 = (8.33,3.66, 0) is obtained with the values of Table 2 satisfying the global stability condition See Appendix G with the eigenvalues in Table 6.E 3 = (7.1029,0, 21.3839) arises by the parameter values from Table 2 satisfying the global stability condition See Appendix G with the eigenvalues provided in Table 6.
The Hopf bifurcation discussed in Section 3.4.5.1 is numerically validated by the parameters h = 0.513, m = 1.2, µ = 0.01 and the other ones from Table 5.By varying the values of λ we obtain Figure 2. Increasing the value of λ to 1.8 while keeping the same values for the remaining parameters, all the feasibility conditions of the coexistence equilibrium are satisfied, giving E * (8.989, 1.151, 27.708).Here λ min = 1.8, λ cr = 3 and λ max = 4.4.E * is asymptotically stable when λ min = 1.8 < λ cr .When λ lies between λ cr and the maximum value of λ cr , a bifurcating limit cycle occurs from a Hopf bifurcation; Fig. 2.Here Ω = 10.49;see Section 3.4.5.6.Therefore the Hopf bifurcation is subcritical.An extensive numerical simulation shows that when the value of λ is very close to λ cr , the three populations S, I and P take long time to stabilize whereas when λ crosses the value λ cr coming very close to λ max , the three populations become unstable.Thus in summary we have the following proposition.
Proposition 5.The interval [λ min , λ max ] contains a critical value λ cr , where a subcritical Hopf bifurcation occurs.For the Hopf bifurcation in the interval [λ min , λ max ], the interior equilibrium point is asymptotically stable between λ min and λ cr , and for λ lying between λ cr and λ max , a limit cycle occurs.On the other hand, for λ > λ max , the equilibrium does not exist.
The basic reproduction number depends on the ratio between infection and prey mortality multiplied by the square root of the system carrying capacity K.The disease will not die out if the prey mortality increases due just to infection rather than the natural predation rate.

Interpretation
Note that the trivial equilibrium is always unstable.This indicates that each subpopulation cannot disappear, thereby preserving the ecosystem.
The equilibrium E 1 is locally asymptotically stable if (8) holds, which biologically can be expressed by saying that the carrying capacity is bounded by a combination of the prey mortality rate and disease transmission.The bound is directly proportional to the former and inversely to the latter.In addition, it contains the equilibrium level of the predator-free point, which indicates the possibility of a transcritical bifurcation among the two equilibria.
Asymptotic stability of E 1 occurs for R 0 < 1; see ( 8).E 2 is feasible for R 0 > 1; see (4).These are opposite conditions so that a transcritical bifurcation occurs for which E 2 emanates from E 1 when the latter becomes unstable.Also, E 2 is stable for λ 1 , which gives a condition for disease control.
The disease-free equilibrium E 3 is feasible when the predator search efficiency h = (α √ K) −1 lies in the interval (h 1 ); see (5).We find numerically the disease-free system experiences a Hopf bifurcation with bifurcation parameter h.E 3 is stable for h < h 1 (see (10)), and the coexistence equilibrium is feasible for h > h [3] 1 (see (7)) so that there exists a transcritical bifurcation for which the coexistence equilibrium E * emanates from the disease-free equilibrium E 3 when the latter becomes unstable; recall Section 3.4.5.4.Also, for feasibility condition of the interior equilibrium, the force of infection must lie in the interval (λ 1 ); see (7).This gives the condition for disease control.Based on these results, for all the equilibrium points, we underline the relevant role played by the parameter λ for the stability of points E 1 and E 2 , while the parameter h is critical for achieving the disease-free equilibrium.The coexistence equilibrium evolution is regulated by the force of infection.Below a certain value of λ, the system may become extinct; Section 4.3.
In addition, the equilibria global stability conditions reported in Table 3 are obtained using Lyapunov and Lasalle theorem of Section 3.4.5.4.https://www.journals.vu.lt/nonlinear-analysis The crucial role of the parameters λ, µ, θ 2 , δ, h and r in the control of the dynamical behavior of the system is therefore apparent.A future direction of the present work therefore can be well extended by introducing recovery of disease incorporating the prey refuge.

Bifurcations
We now investigate the bifurcations and corresponding limit point diagrams with respect to the disease transmission λ.
In each bifurcation diagram in the λS, λI, λP -planes, there exists a complete loop on the right half part.It has two branches, one of which is stable (blue), and the other is unstable (red) branch.The loop joins two different equilibria, thereby, it is a heteroclinic loop.This can be interpreted by saying that the heteroclinic point has different past and future.Also, the existence of a heteroclinic orbit for the critical parameter value related to the force of infection implies biological overexploitation by the disease [21].This shows that the force of infection affects the ultimate behavior of both system populations, determining their survival or extinction.

A comparison
The model presented here is closely related to the one of [23], where instead of the infected, it is the susceptible prey that behave individually.In both cases, the ecosystem cannot disappear, a result that from the biodiversity point of view is good.It is mainly due to the fact that environment has always means to support the prey, in particular healthy prey, by providing them enough feeding resources.
In [23] the predator-free point harbors always the disease endemically.In addition, we find the coexistence equilibrium in which also predators thrive, but where again the disease is not eradicated.
In this system, however, the disease can be eradicated, while the healthy prey is preserved at equilibrium E 1 .Its feasibility and stability conditions provide the theoretical tools to achieve such goal if needed.
On the other hand, if the predators constitute a nuisance and should be eliminated, both here and [23] contain the predator-free point.The role of the relevant parameters in controlling the possible system outcomes have been elucidated and can be obtained from the tables provided in the previous sections.
The predator-free equilibrium becomes unstable for a very low predator mortality rate δ < δ [1] , where δ [1] = θ 1 α √ K(1 + hα √ K) −1 ; see (8).A similar result is found for the disease in predators [4,9], where environmental carrying capacity K and the predation rate play an essential role.We find a supercritical pitchfork bifurcation around E 1 at K. In addition, the predator's average handling time h influences the threshold level.In [3], it is shown that the stability of the prey-only equilibrium in a predator-prey model with disease in predator changes when the predator mortality rate exceeds a threshold value.We have also similar results but for a large enough prey mortality rate, µ > µ [1] , where λ influences the threshold value, which represents a different result from [3].
The predator-free point becomes unstable in the presence of the disease in [23].But in this system the predator-free point is stable for 1 < R 0 < R [3] 0 , where R (9).In [2,4,9] the disease-free equilibrium is stable if the deaths of infected predator lie above a threshold value, where only susceptible predators thrive.A similar situation arises in [23], where the stability of the disease-free system required the predator mortality rate must fall below a threshold so that the predators invade the environment permanently.In this system, we find both upper and lower threshold value for predator mortality rate δ.For the feasibility of the equilibrium point, it must lie in the interval (δ 1 ); see (5).The stability of E 3 holds for small enough predator's mortality rate, δ < δ [3] 1 ; see (10).The predator average handling time h plays a role in both δ [3] 0 and δ [3] 1 .A healthy predatorprey system becomes stable for predation rate above a threshold value in [3], but we get an opposite condition: the disease-free system is feasible for low prey search efficiency α < α [3] , where α [3] = δ[ √ K(θ 1 − δh)] −1 .Both prey-only and predator-free equilibria are unconditionally unstable in [24].Instead, we get suitable stability conditions.
The different roles of the various parameters are shown in different equilibrium points.However, λ has the most crucial role for which a heteroclinic orbit in the coexistence equilibrium arises.It leads to extinction when the force of infection falls below the threshold λ [ * ] 0 ; see (7).From the ecological point of view such a situation can arise due to overexploitation of force of infection discussed in Section 4.3.Both the theoretical and numerical analysis are found important to draw conclusions on a general level.

Appendix B: Proof of Proposition 3
The Jacobian eigenvalues at E * are the roots of the cubic equation .
By the Routh-Hurwitz criteria asymptotic stability is achieved whenever b if we choose j 0 = j 1 = j 2 for all cases when Since the value of Θ(S, I, P ) is positive at all boundary points, system (2) is persistent.

Table 2 .
Changes of various parameters values in the different equilibrium points.

Table 3 .
Analytic conditions for feasibility.

Table 5 .
Variables and parameters used during simulations.