Eco-Epidemiological Model of Predator-Prey with Two-Strain Infections: The Impact of Herd Behavior

1. Introduction

Predator-prey dynamics are critical in maintaining ecological balance, influenced by various factors such as resource availability, environmental conditions, and predation efficiency [1,2,3,4]. However, the impact of diseases on these interactions [5,6,7,8,9,10,11], especially multi-strain infections [13] within prey populations, has been less explored [14,24,25,26]. When prey are infected, their ability to evade predators declines, not only due to physical weakening but also through behavioral changes such as the formation of herds. These herds serve as protective mechanisms [21,22,23,30], but in the presence of infection, they can introduce additional complexity to the predator-prey relationship [12,29,31,32,34].

In this study, we develop an eco-epidemiological model that investigates how two distinct strains of pathogens interact with herding behavior in prey populations, while predators remain unaffected by the infections themselves. This multi-strain dynamic introduces nontrivial challenges in understanding population stability and predator efficiency. Through this model, we explore how herd shape and predator mortality rates, driven by the protective behavior of herds, influence disease dynamics and population survival.

The real-world applicability of this model is significant. Wildlife populations such as deer, wild boars, and fish frequently experience multi-strain infections that affect both prey behavior and predator efficiency [35,36]. Similarly, agricultural systems [39], where livestock herds [37,38] encounter various pathogens, offer another layer of complexity. In marine ecosystems, schooling fish serve as prey for larger predators, where infection can compromise their defense mechanisms [41,42,43]. Furthermore, urban wildlife and zoonotic diseases [33,47,48,49,50] provide an additional domain where multi-strain infections and human interaction complicate traditional predator-prey dynamics [44,45,46]. The model thus offers critical insights for managing disease outbreaks, ensuring wildlife conservation, and predicting shifts in population stability under varying environmental conditions.

Some researchers have considered two strains of infection either in prey or predator [24,25,26,27,28] but could not show the coexistence of both strains of infections. Martcheva M. [51] considered two types of model with susceptible and two strain infected prey and one with generalist predator, second with specialist predator. She showed that increasing predation leads to eradication of disease without extinction of prey population and also that the presence of specialist predators may lead to the coexistence of both strains. To our knowledge, no one has studied two strains of infection in prey with prey herd and predator mortality due to prey herd. Here we have also assumed both the strains of infections don’t affect each other. That means there is neither coinfection nor superinfection,i.e. when an individual species gets infected by one strain, it cannot develop the second strain and have both strains of infection at the same time nor the second strain can replace the first strain.

We focus on analyzing the stability of populations through well-posedness and boundedness evaluations. Using bifurcation analysis, we reveal the critical thresholds at which the prey population undergoes significant dynamic shifts—highlighting the occurrence of supercritical and subcritical Hopf bifurcations, and their roles in population oscillations or collapse. This investigation of bifurcations provides valuable insights into how changes in herd formation or predator mortality due to herd protection can trigger these complex behaviors. We show the presence of the limit point bifurcation of limit cycles, emerging from the Generalized Hopf bifurcation point. This study contributes to a deeper understanding of eco-epidemiological interactions involving multiple infection strains and predator-prey dynamics, laying the groundwork for future models that can incorporate additional environmental and behavioral complexities.

2. Mathematical Model Formulation

The system model is formed based on the following assumptions:

• The susceptible prey grows logistically in the absence of infected prey and predator and prey infected with both strains neither contribute to the reproduction nor to the carrying capacity thus we have,

  $$\begin{array}{c}\hfill S^{\prime }=rS(1-{\displaystyle \frac{S}{K}})\end{array}$$

  (1)
• Infection of both strains within the prey population occurs horizontally via direct contact following the mass action incidence i.e.

  $$\begin{array}{c}\hfill H(S,I_j)={\beta }_j{I}_{j}S\mathrm{for}j=1,2\end{array}$$

  (2)

  where ${\beta }_j$ are the rate of transmission of strain j infection.
• Only susceptible prey participates in the herding group where the rate of herd shape is $0<k<1$
  • if the herd shape is circle or square, $k=0.5$
  • if the herd shape is cube or sphere, $k=0.67$.
• The infected prey dies out either due to disease by strain 1 or 2 or naturally (${\mu }_j\mathrm{for}j=1,2$ denotes the total death rate)and upon recovery(${\gamma }_j\mathrm{for}j=1,2$ denotes the recovery rates) from the infection they re-enter the susceptible population.
• Predators attack the susceptible prey according to the Holling-II interaction functional type created by S.Djilali [30] given as

  $$\begin{array}{c}\hfill H(S,P)={\displaystyle \frac{a{S}^{k}P}{1+ah{S}^k}}\end{array}$$

  (3)

  and the prey infected with strains 1 and 2 according to the Holling type-I functional response given as

  $$\begin{array}{c}\hfill H(I_j,P)={\beta }_j{I}_{j}P\mathrm{for}j=1,2.\end{array}$$

  (4)
• Predators attack the boundary of the herd and the prey on the outer herd can injure the predators leading to their death given as $\delta S^{k}P$.

Using the above postulates, we construct a four-dimensional eco-epidemiological model:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}& \hfill =rS\left(1-{\displaystyle \frac{S}{K}}\right)-[{\beta }_1{I}_1+{\beta }_2{I}_2]S+{\gamma }_1{I}_1+{\gamma }_2{I}_2-{\displaystyle \frac{a{S}^{k}P}{1+ah{S}^k}}={\varphi }_1(S,I_1,I_2,P)\hfill \\ \hfill {\displaystyle \frac{d{I}_1}{dt}}& \hfill ={\beta }_1{I}_{1}S-({\gamma }_1+{\mu }_1)I_1-{\alpha }_1{I}_{1}P={\varphi }_2(S,I_1,I_2,P)\hfill \\ \hfill {\displaystyle \frac{d{I}_2}{dt}}& \hfill ={\beta }_2{I}_{2}S-({\gamma }_2+{\mu }_2)I_2-{\alpha }_2{I}_{2}P={\varphi }_3(S,I_1,I_2,P)\hfill \\ \hfill {\displaystyle \frac{dP}{dt}}& \hfill ={\displaystyle \frac{ae{S}^{k}P}{1+ah{S}^k}}+e_1{\alpha }_1{I}_{1}P+e_2{\alpha }_2{I}_{2}P-\delta S^{k}P-\mu P={\varphi }_4(S,I_1,I_2,P)\hfill \end{array}$$

(5)

we have considered a system of equations describing a prey-predator model with two strains of infection running in prey species and prey group forming a herd against the predator.

In the first equation we assume the logistic growth rate of susceptible prey with r as the intrinsic growth rate and K as carrying capacity. The susceptible prey gets infected by strain 1 and strain 2 with rates of transmission ${\beta }_1,{\beta }_2$ respectively. The susceptible prey forms a herd and the predator attacks on the boundary of the herd with a Holling type-II interaction functional response.

In the second and third equations, we have the interaction of susceptible prey and the prey infected with strain 1 and strain 2 assuming that the infected can be recovered with rates ${\gamma }_1,{\gamma }_2$ and join the class of susceptible prey and die with total death rates ${\mu }_1,{\mu }_2$ respectively. Also the predator attacks both the prey infected with strain 1 and strain 2 with linear rates ${\alpha }_1,{\alpha }_2$ respectively.

In the fourth equation, we have the growth of predators by consuming susceptible and infected prey with conversion efficiency rates $e,e_1,e_2$. Also, we have considered the mortality of predators due to prey herd($\delta$) and their natural death($\mu )$.

Figure 1. Schematic flowchart for the model.

Figure 1. Schematic flowchart for the model.

3. Description of Parameters

Parameters	Description
S	Susceptible Prey
$I_1$	Prey Infected with strain 1
$I_2$	Prey Infected with strain 2
P	Predator
r	Intrinsic growth rate of susceptible prey
K	carrying capacity of susceptible prey
${\beta }_1$	rate of transmission of strain 1 infection
${\beta }_2$	rate of transmission of strain 2 infection
${\gamma }_1$	recovery rate from strain 1
${\gamma }_2$	recovery rate from strain 2
a	rate of predation on susceptible prey
k	rate of prey herd shape
h	handling time spent by a predator on healthy prey
${\alpha }_1$	rate of predation on prey infected with strain 1
${\alpha }_2$	rate of predation on prey infected with strain 2
e	conversion rate of susceptible prey to predator
$e_1$	conversion rate of prey infected with strain 1 to predator
$e_2$	conversion rate of prey infected with strain 2 to predator
$\delta$	mortality rate of a predator due to prey herd
$\mu$	natural death rate of predator
${\mu }_1$	total death rate of prey infected with strain 1
${\mu }_2$	total death rate of prey infected with strain 2

4. Well-Posedness of the Formulated Model

Theorem 1.

All the solutions of the system (5) under the given initial condition are unique and positive in R.

Proof of Theorem 1. 

The function ${\varphi }_1,{\varphi }_2,{\varphi }_3$ on the right hand side of system of equation(5) are locally Lipschitz Continuous in the region ${\rm Y}=\{(S,I_1,I_2,P)\in [\xi ,A]\times [0,B]\times [0,B]\times [0,C]\}$ with $\xi$ a positive real constant. This can be seen $\forall S_1,S_2\in {\rm Y}$ we have,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill |{\varphi }_1(S_1,I_1,I_2,P)-{\varphi }_1(S_2,I_1,I_2,P)|& \hfill \le r|S_1-S_2|+{\displaystyle \frac{r}{K}}|S_1-S_2\left|\right|S_1+S_2|+({\beta }_1|I_1|+{\beta }_2|I_2\left|\right)|S_1-S_2|\hfill \\ & \hfill +{\displaystyle \frac{a\left|P\right||S_1^k-S_2^k\left|\right|S_1-S_2|}{|1+ah{S}_1^k\left|\right|1+ah{S}_2^k\left|\right|S_1-S_2|}}\hfill \\ \hfill & \hfill \le |S_1-S_2\left|\right(.r+{\displaystyle \frac{r\left(\right|S_1|+|S_2\left|\right)}{K}}+{\beta }_1|I_1|+{\beta }_2\left|I_2\right|\hfill \\ & \hfill +{\displaystyle \frac{a\left|P\right|}{|1+ah{S}_1^k\left|\right|1+ah{S}_2^k|}}\left|{\displaystyle \frac{S_1^k-S_2^k}{S_1-S_2}}\right|.)\hfill \\ \hfill & \hfill \le L|S_1-S_2|\hfill \end{array}\end{array}$$

where $L=r+{\displaystyle \frac{2rA}{K}}+B({\beta }_1+{\beta }_2)+akC{q}^{k-1}>0$ is a Lipschitz constant.

Thus ${\varphi }_1(S,I_1,I_2,P)$ is locally Lipschitz-continuous with respect to S. Similarly, it can be shown for $I_1,I_2,P$. Hence the system has a unique solution in the region ${\rm Y}$.

Now, to show we have positive solutions, we will show $S\left(t\right)>0\forall t>0$ by contradiction. Let us assume there exists a $\tau >0$ with $\tau =inf\{t:S(t)=0,t>0\}$,such that ${\displaystyle \frac{dS}{dt}}{|}_{t=\tau }<0\mathrm{and}S\left(t\right)>0\forall t\in [0,\tau )$. From the first equation of (5) we have,

$$\begin{array}{cc}\hfill {\displaystyle \frac{dS}{dt}}{|}_{t=\tau }& \hfill =rS\left(\tau \right)\left(1-{\displaystyle \frac{S\left(\tau \right)}{K}}\right)-[{\beta }_1{I}_1\left(\tau \right)+{\beta }_2{I}_2\left(\tau \right)]S\left(\tau \right)+{\gamma }_1{I}_1\tau )+{\gamma }_2{I}_2\left(\tau \right)-{\displaystyle \frac{a{S}^k\left(\tau \right)P\left(\tau \right)}{1+ah{S}^k\left(\tau \right)}}\hfill \\ & \hfill ={\gamma }_1{I}_1\left(\tau \right)+{\gamma }_2{I}_2\left(\tau \right)\hfill \end{array}$$

we can write the second, third, fourth equations of the model (5) as:

${\displaystyle \frac{d{I}_1}{dt}}=I_1{f}_1(S,I_1,I_2,P),{\displaystyle \frac{d{I}_2}{dt}}=I_2{f}_2(S,I_1,I_2,P),{\displaystyle \frac{dP}{dt}}=P{f}_3(S,I_1,I_2,P)$

where

$f_1(S,I_1,I_2,P)={\beta }_{1}S-{\alpha }_{1}P-{\gamma }_1-{\mu }_1$

$f_2(S,I_1,I_2,P)={\beta }_{2}S-{\alpha }_{2}P-{\gamma }_2-{\mu }_2$

$f_3(S,I_1,I_2,P)={\displaystyle \frac{ea{S}^k}{1+ah{S}^k}}+e_1{\alpha }_1{I}_1+e_2{\alpha }_2{I}_2-\delta S^k-\mu$

thus it follows,

$I_1\left(t\right)=I_1\left(0\right)e^{{\int }_0^t{f}_1(S\left(u\right),I_1\left(u\right),I_2\left(u\right),P\left(u\right))du}\ge 0$

$I_2\left(t\right)=I_2\left(0\right)e^{{\int }_0^t{f}_2(S\left(u\right),I_1\left(u\right),I_2\left(u\right),P\left(u\right))du}\ge 0$

$P\left(t\right)=P\left(0\right)e^{{\int }_0^t{f}_3(S\left(u\right),I_1\left(u\right),I_2\left(u\right),P\left(u\right))du}\ge 0$

Now we have four cases:

$\left(1\right)I_1\left(0\right)=I_2\left(0\right)=0$ then

$I_1\left(t\right)=I_2\left(t\right)=0\forall t>0\Rightarrow I_1\left(\tau \right)=I_2\left(\tau \right)=0$

$\Rightarrow {\displaystyle \frac{dS}{dt}}{|}_{t=\tau }=0$ which contradicts the condition ${\displaystyle \frac{dS}{dt}}{|}_{t=\tau }<0$.

$\Rightarrow S\left(t\right)>0\forall t\ge 0$

$\left(2\right)I_1\left(0\right)>0and{I}_2\left(0\right)>0$ then

$I_1\left(t\right)>0and{I}_2\left(t\right)>0\forall t>0\Rightarrow I_1\left(\tau \right)>0and{I}_2\left(\tau \right)>0$

$\Rightarrow {\displaystyle \frac{dS}{dt}}{|}_{t=\tau }>0$ which contradicts the condition ${\displaystyle \frac{dS}{dt}}{|}_{t=\tau }<0$.

$\Rightarrow S\left(t\right)>0\forall t\ge 0$

Similarly, we can prove for the rest two cases i.e. $I_1\left(0\right)=0,I_2\left(0\right)>0\mathrm{and}{I}_1\left(0\right)>0,I_2\left(0\right)=0$.

Thus all solutions $(S\left(t\right),I_1\left(t\right),I_2\left(t\right),P\left(t\right))$ with the stated initial condition remain positive for $\forall t>0$. □

Theorem 2.

All solutions of the system (5) beginning in $R_{+}^4$ stay enclosed in the region ${\rm Y}=\{(S,I_1,I_2,P):0<S+I_1+I_2\le {\displaystyle \frac{{(r+{\zeta }_1)}^{2}K}{4r{\zeta }_1}},0<S+{\displaystyle \frac{1}{e}}P\le {\displaystyle \frac{{(r+{\zeta }_2)}^{2}K}{4r{\zeta }_2}}+{\gamma }_1{\displaystyle \frac{{(r+{\zeta }_1)}^{2}K}{4r{\zeta }_1{\zeta }_2}}+{\gamma }_2{\displaystyle \frac{{(r+{\zeta }_1)}^{2}K}{4r{\zeta }_1{\zeta }_2}}$, where ${\zeta }_1\mathrm{and}{\zeta }_2$ are some real numbers satisfying $0<{\zeta }_1<min\{{\mu }_1,{\mu }_2\},0<{\zeta }_2<\mu -{\displaystyle \frac{(e_1{\alpha }_1+e_2{\alpha }_2){(r+{\zeta }_1)}^{2}K}{4r{\zeta }_1}}$

Proof of Theorem 2. 

From the first equation of the system, we have $\underset{t\to \infty }{lim}\sup S\left(t\right)\le K.$

Let ${\varphi }_1=S+I_1+I_2$ and ${\zeta }_1>0$. Then

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{d{\varphi }_1}{dt}}+{\zeta }_1{\varphi }_1& \hfill =rS(1-{\displaystyle \frac{S}{K}})-{\displaystyle \frac{a{S}^{k}P}{1+ah{S}^k}}-{\alpha }_1{I}_{1}P-{\alpha }_2{I}_{2}P-{\mu }_1{I}_1-{\mu }_2{I}_2+{\zeta }_1(S+I_1+I_2)\hfill \\ \hfill & \hfill \le rS(1-{\displaystyle \frac{S}{K}})-{\mu }_1{I}_1-{\mu }_2{I}_2+{\zeta }_1(S+I_1+I_2)\hfill \\ \hfill & \hfill \le (r+{\zeta }_1)S-{\displaystyle \frac{r{S}^2}{K}}-({\mu }_1-{\zeta }_1)I_1-({\mu }_2-{\zeta }_1)I_2\hfill \\ \hfill & \hfill \le (r+{\zeta }_1)S-{\displaystyle \frac{r{S}^2}{K}},\mathrm{for}{\zeta }_1<min\{{\mu }_1,{\mu }_2\}\hfill \end{array}\end{array}$$

since $\underset{S\ge 0}{max}(r+{\zeta }_1)S-{\displaystyle \frac{r{S}^2}{K}}={\displaystyle \frac{{(r+{\zeta }_1)}^{2}K}{4r}}$.

Therefore, ${\displaystyle \frac{d{\varphi }_1}{dt}}+{\zeta }_1{\varphi }_1\le {\displaystyle \frac{{(r+{\zeta }_1)}^{2}K}{4r}}\Rightarrow \underset{t\to \infty }{lim}\sup (S\left(t\right)+I_1\left(t\right)+I_2\left(t\right))\le {\displaystyle \frac{{(r+{\zeta }_1)}^{2}K}{4r{\zeta }_1}}$.

Now, to check boundedness of $P\left(t\right)$ we assume ${\varphi }_2=S+{\displaystyle \frac{1}{e}}P$ and ${\zeta }_2>0$. Then

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle \frac{d{\varphi }_2}{dt}}+{\zeta }_2{\varphi }_2& \hfill =rS(1-{\displaystyle \frac{S}{K}})-{\beta }_1{I}_{1}S-{\beta }_2{I}_{2}S+{\gamma }_1{I}_1+{\gamma }_2{I}_2+{\displaystyle \frac{e_1}{e}}{\alpha }_1{I}_{1}P+{\displaystyle \frac{e_2}{e}}{\alpha }_2{I}_{2}P\hfill \\ & \hfill -{\displaystyle \frac{\delta }{e}}{S}^{k}P-{\displaystyle \frac{\mu }{e}}P+{\zeta }_{2}S+{\displaystyle \frac{{\zeta }_2}{e}}P\hfill \\ \hfill & \hfill \le (r+{\zeta }_2)S-{\displaystyle \frac{r{S}^2}{K}}+{\gamma }_1{I}_{1max}+{\gamma }_2{I}_{2max}-({\displaystyle \frac{\mu }{e}}-{\displaystyle \frac{{\zeta }_2}{e}}-{\displaystyle \frac{e_1}{e}}{\alpha }_1{I}_{1max}-{\displaystyle \frac{e_2}{e}}{\alpha }_2{I}_{2max})P\hfill \end{array}\end{array}$$

where $I_{1max}=I_{2max}={\displaystyle \frac{{(r+{\zeta }_1)}^{2}K}{4r{\zeta }_1}}$.

Now for $0<{\zeta }_2<\mu -e_1{\alpha }_1{I}_{1max}-e_2{\alpha }_2{I}_{2max}$, we have

${\displaystyle \frac{d{\varphi }_2}{dt}}+{\zeta }_2{\varphi }_2\le (r+{\zeta }_2)S-{\displaystyle \frac{r{S}^2}{K}}+{\gamma }_1{I}_{1max}+{\gamma }_2{I}_{2max}$

which implies $\underset{t\to \infty }{lim}\sup (S\left(t\right)+{\displaystyle \frac{1}{e}}P\left(t\right))\le {\displaystyle \frac{{(r+{\zeta }_2)}^{2}K}{4r{\zeta }_2}}+{\gamma }_1{\displaystyle \frac{{(r+{\zeta }_1)}^{2}K}{4r{\zeta }_1{\zeta }_2}}+{\gamma }_2{\displaystyle \frac{{(r+{\zeta }_1)}^{2}K}{4r{\zeta }_1{\zeta }_2}}$.

Hence the theorem follows. □

5. Equilibrium Points

The model system (5) has following feasible equilibrium points:

• $E_1=(0,0,0,0)$
• $E_2=(K,0,0,0)$
  The basic reproduction number at the disease-free equilibrium point is
  $R_0=max\{R_0^1,R_0^2\}$ where $R_0^1={\displaystyle \frac{{\beta }_{1}K}{{\gamma }_1+{\mu }_1}},R_0^2={\displaystyle \frac{{\beta }_{2}K}{{\gamma }_2+{\mu }_2}}$
• $E_3=(S_{SP},0,0,P_{SP})$
  where $S_{SP}\mathrm{and}{P}_{SP}$ are solutions of equations:
  $ah\delta S^{2k}+S^k(\delta +\mu ah-ae)+\mu =0$ and $rS(1-{\displaystyle \frac{S}{K}})-{\displaystyle \frac{a{S}^{k}P}{1+ah{S}^k}}=0$
  For positive solutions, the following conditions should hold:
  ${(\mu ah+\delta -ae)}^2\ge 4ah\delta \mu$ and $\delta +\mu ah<ae$
• $E_4=(S_1,I_1,0,0)$
  where $S={\displaystyle \frac{{\mu }_1+{\gamma }_1}{{\beta }_1}},I_1={\displaystyle \frac{rS}{{\mu }_1}}(1-{\displaystyle \frac{S}{K}})$
• $E_5=(S_2,0,I_2,0)$
  where $S={\displaystyle \frac{{\mu }_2+{\gamma }_2}{{\beta }_2}},I_2={\displaystyle \frac{rS}{{\mu }_2}}(1-{\displaystyle \frac{S}{K}})$
• $E_6=(S_{I_{1}P},I_1,0,P)$
  where $I_1={\displaystyle \frac{1}{e_1{\alpha }_1}}\left(\delta S^k+\mu -{\displaystyle \frac{ea{S}^k}{1+ah{S}^k}}\right)$, $P={\displaystyle \frac{{\beta }_{1}S-{\gamma }_1-{\mu }_1}{{\alpha }_1}}$
  $f_1\left(S\right)=A_1{S}^{2k+1}+A_2{S}^{2k}+A_3{S}^{k+2}+A_4{S}^{k+1}+A_5{S}^k+A_6{S}^2+A_{7}S+A_8$

  $$\begin{array}{c}\hfill \mathrm{where}{A}_1=\delta ahK{\beta }_1,A_2=-\delta ah{\gamma }_{1}K,A_3=ahr{e}_1{\alpha }_1,A_4={\beta }_1{e}_{1}aK-aeK{\beta }_1+K{\beta }_1\delta +\mu ahK{\beta }_1-r{e}_1{\alpha }_{1}ahK,\\ \hfill A_5=ae{\gamma }_{1}K-\delta {\gamma }_{1}K-\mu ah{\gamma }_{1}K-({\gamma }_1+{\mu }_1)e_{1}aK,A_6=r{e}_1{\alpha }_1,A_7=K(\mu {\beta }_1-r{e}_1{\alpha }_1),A_8=-\mu {\gamma }_{1}K\end{array}$$
• $E_7=(S_{I_{2}P},0,I_2,P)$
  where $I_1={\displaystyle \frac{1}{e_2{\alpha }_2}}(\delta S^k+\mu -{\displaystyle \frac{ea{S}^k}{1+ah{S}^k}})$, $P={\displaystyle \frac{{\beta }_{2}S-{\gamma }_2-{\mu }_2}{{\alpha }_2}}$
  $f_2\left(S\right)=A_1{S}^{2k+1}+A_2{S}^{2k}+A_3{S}^{k+2}+A_4{S}^{k+1}+A_5{S}^k+A_6{S}^2+A_{7}S+A_8$

  $$\begin{array}{c}\hfill \mathrm{where}{A}_1=\delta ahK{\beta }_2,A_2=-\delta ah{\gamma }_{2}K,A_3=ahr{e}_2{\alpha }_2,A_4={\beta }_2{e}_{2}aK-aeK{\beta }_2+K{\beta }_2\delta +\mu ahK{\beta }_2-r{e}_2{\alpha }_{2}ahK,\\ \hfill A_5=ae{\gamma }_{2}K-\delta {\gamma }_{2}K-\mu ah{\gamma }_{2}K-({\gamma }_2+{\mu }_2)e_{2}aK,A_6=r{e}_2{\alpha }_2,A_7=K(\mu {\beta }_2-r{e}_2{\alpha }_2),A_8=-\mu {\gamma }_{2}K.\end{array}$$
• $E^{*}=(S^{*},I_1^{*},I_2^{*},P^{*})$ where $P^{*}={\displaystyle \frac{{\beta }_{1}S-{\mu }_1-{\gamma }_1}{{\alpha }_1}}={\displaystyle \frac{{\beta }_{2}S-{\mu }_2-{\gamma }_2}{{\alpha }_2}}$
  $\Rightarrow S^{*}={\displaystyle \frac{{\alpha }_2({\mu }_1+{\gamma }_1)-{\alpha }_1({\mu }_2+{\gamma }_2)}{{\beta }_1{\alpha }_2-{\alpha }_1{\beta }_2}}$

  $$\begin{array}{cc}\hfill I_1^{*}& \hfill ={\displaystyle \frac{K{\alpha }_1({\gamma }_2-{\beta }_{2}S)(\delta S^k+\mu +ah\delta S^{2k}+\mu ah{S}^k-ae{S}^k)-e_2{\alpha }_2\left({\beta }_{1}aK{S}^{k+1}\right)+e_2{\alpha }_2({\mu }_1+{\gamma }_1)Ka{S}^k+e_2{\alpha }_2{\alpha }_{1}rS(1+ah{S}^k)(K-S)}{K{\alpha }_1(1+ah{S}^k)(e_1{\alpha }_1({\gamma }_2-{\beta }_{2}S)-e_2{\alpha }_2({\gamma }_1-{\beta }_{1}S))}}\hfill \\ \hfill I_2^{*}& \hfill ={\displaystyle \frac{K{\alpha }_1({\gamma }_1-{\beta }_{1}S)(\delta S^k+\mu +ah\delta S^{2k}+\mu ah{S}^k-ae{S}^k)-e_1{\alpha }_1\left({\beta }_{1}aK{S}^{k+1}\right)+e_1{\alpha }_1({\mu }_1+{\gamma }_1)Ka{S}^k+e_1{\alpha }_1^{2}rS(1+ah{S}^k)(K-S)}{K{\alpha }_1(1+ah{S}^k)(e_2{\alpha }_2({\gamma }_1-{\beta }_{1}S)-e_1{\alpha }_1({\gamma }_2-{\beta }_{2}S))}}\hfill \end{array}$$

6. Stability Conditions

• $E_1(0,0,0,0)$ is a saddle point because the eigenvalues of variational matrix are
  $(r,-{\gamma }_1-{\mu }_1,-{\gamma }_2-{\mu }_2,-m)$
• $E_2(K,0,0,0)$ is locally asymptotically stable if
  ${\displaystyle \frac{K^{k}ae}{K^{k}ah+1}}<K^{k}d+m$ and $K{\beta }_1<{\mu }_1+{\gamma }_1$ and $K{\beta }_2<{\mu }_2+{\gamma }_2$
• for $E^{*}$ the variational matrix is

  $$\begin{array}{c}\hfill J=\left[\begin{array}{cccc}{J}_{11}& J_{12}& J_{13}& J_{14}\\ J_{21}& J_{22}& J_{23}& J_{24}\\ J_{31}& J_{32}& J_{33}& J_{34}\\ J_{41}& J_{42}& J_{43}& J_{44}\end{array}\right]\end{array}$$

  where $J_{11}={\displaystyle \frac{a^{2}hk{S}^{*(2k-1)}{P}^{*}}{{(1+ah{S}^{*k})}^2}}-{\beta }_2{I}_2^{*}+r(1-{\displaystyle \frac{S^{*}}{K}})-{\displaystyle \frac{r{S}^{*}}{K}}-{\displaystyle \frac{ak{S}^{*(k-1)}{P}^{*}}{1+ah{S}^{*k}}}-{\beta }_1{I}_1^{*}$
  $J_{12}={\gamma }_1-{\beta }_1{S}^{*}$, $J_{13}={\gamma }_2-{\beta }_2{S}^{*}$, $J_{14}=-{\displaystyle \frac{a{S}^{*k}}{1+ah{S}^{*k}}}$
  $J_{21}={\beta }_1{I}_1^{*}$, $J_{22}={\beta }_1{S}^{*}-{\mu }_1-{\alpha }_1{P}^{*}-{\gamma }_1$, $J_{23}=0$, $J_{24}=-{\alpha }_1{I}_1^{*}$
  $J_{31}={\beta }_2{I}_2^{*}$, $J_{32}=0$, $J_{33}={\beta }_2{S}^{*}-{\mu }_2-{\alpha }_2{P}^{*}-{\gamma }_2$, $J_{34}=-{\alpha }_2{I}_2^{*}$
  $J_{41}={\displaystyle \frac{aek{S}^{*(k-1)}{P}^{*}}{1+ah{S}^{*k}}}-dk{S}^{*(k-1)}{P}^{*}-{\displaystyle \frac{a^{2}ehk{S}^{*(2k-1)}{P}^{*}}{{(1+ah{S}^{*k})}^2}}$
  $J_{42}={\alpha }_1{e}_1{P}^{*}$, $J_{43}={\alpha }_2{e}_2{P}^{*}$
  $J_{44}={\alpha }_1{e}_1{I}_1^{*}-d{S}^{*k}-m+{\alpha }_2{e}_2{I}_2^{*}+{\displaystyle \frac{ae{S}^{*k}}{1+ah{S}^{*k}}}$
  The characteristic polynomial of matrix J is

  $${\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0$$

  (6)

  where

  $$\begin{array}{cc}\hfill a_1& \hfill =-(J_{11}+J_{22}+J_{33}+J_{44})\hfill \\ \hfill a_2& \hfill =J_{11}{J}_{22}-J_{12}{J}_{21}+J_{11}{J}_{33}-J_{13}{J}_{31}+J_{11}{J}_{44}-J_{14}{J}_{41}+J_{22}{J}_{33}-J_{23}{J}_{32}+J_{22}{J}_{44}-J_{24}{J}_{42}+J_{33}{J}_{44}-J_{34}{J}_{43}\hfill \\ \hfill a_3& \hfill =J_{11}{J}_{23}{J}_{32}-J_{11}{J}_{22}{J}_{33}+J_{12}{J}_{21}{J}_{33}-J_{12}{J}_{23}{J}_{31}-J_{13}{J}_{21}{J}_{32}+J_{13}{J}_{22}{J}_{31}-J_{11}{J}_{22}{J}_{44}+J_{11}{J}_{24}{J}_{42}+J_{12}{J}_{21}{J}_{44}\hfill \\ & \hfill -J_{12}{J}_{24}{J}_{41}-J_{14}{J}_{21}{J}_{42}+J_{14}{J}_{22}{J}_{41}-J_{11}{J}_{33}{J}_{44}+J_{11}{J}_{34}{J}_{43}+J_{13}{J}_{31}{J}_{44}-J_{13}{J}_{34}{J}_{41}-J_{14}{J}_{31}{J}_{43}+J_{14}{J}_{33}{J}_{41}\hfill \\ & \hfill -J_{22}{J}_{33}{J}_{44}+J_{22}{J}_{34}{J}_{43}+J_{23}{J}_{32}{J}_{44}-J_{23}{J}_{34}{J}_{42}-J_{24}{J}_{32}{J}_{43}+J_{24}{J}_{33}{J}_{42}\hfill \\ \hfill a_4& \hfill =J_{11}{J}_{22}{J}_{33}{J}_{44}-J_{11}{J}_{22}{J}_{34}{J}_{43}-J_{11}{J}_{23}{J}_{32}{J}_{44}+J_{11}{J}_{23}{J}_{34}{J}_{42}+J_{11}{J}_{24}{J}_{32}{J}_{43}-J_{11}{J}_{24}{J}_{33}{J}_{42}-J_{12}{J}_{21}{J}_{33}{J}_{44}\hfill \\ & \hfill +J_{12}{J}_{21}{J}_{34}{J}_{43}+J_{12}{J}_{23}{J}_{31}{J}_{44}-J_{12}{J}_{23}{J}_{34}{J}_{41}-J_{12}{J}_{24}{J}_{31}{J}_{43}+J_{12}{J}_{24}{J}_{33}{J}_{41}+J_{13}{J}_{21}{J}_{32}{J}_{44}-J_{13}{J}_{21}{J}_{34}{J}_{42}\hfill \\ & \hfill -J_{13}{J}_{22}{J}_{31}{J}_{44}+J_{13}{J}_{22}{J}_{34}{J}_{41}+J_{13}{J}_{24}{J}_{31}{J}_{42}-J_{13}{J}_{24}{J}_{32}{J}_{41}-J_{14}{J}_{21}{J}_{32}{J}_{43}+J_{14}{J}_{21}{J}_{33}{J}_{42}+J_{14}{J}_{22}{J}_{31}{J}_{43}\hfill \\ & \hfill -J_{14}{J}_{22}{J}_{33}{J}_{41}-J_{14}{J}_{23}{J}_{31}{J}_{42}+J_{14}{J}_{23}{J}_{32}{J}_{41}\hfill \end{array}$$
  According to Routh-Hurwitz criteria $E^{*}$ is locally asymptotically stable if and only if
  ${\mathbf{a}}_{\mathbf{1}}>\mathbf{0},{\mathbf{a}}_{\mathbf{4}}>\mathbf{0},{\mathbf{a}}_{\mathbf{1}}{\mathbf{a}}_{\mathbf{2}}-{\mathbf{a}}_{\mathbf{3}}>\mathbf{0},{\mathbf{a}}_{\mathbf{1}}{\mathbf{a}}_{\mathbf{2}}{\mathbf{a}}_{\mathbf{3}}-{\mathbf{a}}_{\mathbf{3}}^{\mathbf{2}}-{\mathbf{a}}_{\mathbf{4}}{\mathbf{a}}_{\mathbf{1}}^{\mathbf{2}}>\mathbf{0}$

Theorem 3.

The necessary and sufficient conditions for Hopf bifurcation to occur at $E^{*}$ at $k=k^{*}$ are the following:

• $a_1>0,a_4>0,a_1{a}_2-a_3>0$
• $a_1{a}_2{a}_3-a_3^2-a_4{a}_1^2=0$
• $a_1^2\left(k^{*}\right)(a_2{\left(k^{*}\right)}^{\prime }{a}_3\left(k^{*}\right)+a_2\left(k^{*}\right)a_3^{\prime }\left(k^{*}\right)-a_4^{\prime }\left(k^{*}\right)a_1\left(k^{*}\right))-2a_3\left(k^{*}\right)a_3^{\prime }\left(k^{*}\right)a_1\left(k^{*}\right)+a_1^{\prime }\left(k^{*}\right)(a_2\left(k^{*}\right)a_3\left(k^{*}\right)a_1\left(k^{*}\right)+2a_3^2\left(k^{*}\right))\ne 0$

Proof of Theorem 3. 

At $k=k^{*}$ the characteristic equation (6)becomes

$$\left({\lambda }^2\left(k\right)+{\displaystyle \frac{a_3\left(k\right)}{a_1\left(k\right)}}\right)(\lambda \left(k\right)+x_1\left(k\right))(\lambda \left(k\right)+x_2\left(k\right))=0$$

(7)

where $x_1+x_2=a_1,x_1{x}_2=a_2-{\displaystyle \frac{a_3}{a_1}}$

Now to show the transversality condition for Hopf-bifurcation. Let us assume $\lambda \left(k\right)=u\left(k\right)+iv\left(k\right)$. Substituting this in (7) we have,

$$\begin{array}{cc}\hfill D\left(k\right)u{\left(k\right)}^{\prime }-E\left(k\right)v{\left(k\right)}^{\prime }+G\left(k\right)& \hfill =0\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill E\left(k\right)u{\left(k\right)}^{\prime }+D\left(k\right)v{\left(k\right)}^{\prime }+H\left(k\right)& \hfill =0\hfill \end{array}$$

(9)

where

$$\begin{array}{cc}\hfill D\left(k\right)& \hfill =4u^3\left(k\right)-12u\left(k\right)v^2\left(k\right)+3u^2\left(k\right)a_1\left(k\right)-3a_1\left(k\right)v^2\left(k\right)+2a_2\left(k\right)u\left(k\right)+a_3\left(k\right)\hfill \\ \hfill E\left(k\right)& \hfill =12u^2\left(k\right)v\left(k\right)-4v^3\left(k\right)+6u\left(k\right)v\left(k\right)a_1\left(k\right)+2a_2\left(k\right)v\left(k\right)\hfill \\ \hfill G\left(k\right)& \hfill =a_1^{\prime }\left(k\right)(u^3\left(k\right)-3u\left(k\right)v^2\left(k\right))+a_2^{\prime }\left(k\right)(u^2\left(k\right)-v^2\left(k\right))+a_3^{\prime }\left(k\right)u\left(k\right)+a_4^{\prime }\left(k\right)\hfill \\ \hfill H\left(k\right)& \hfill =a_1^{\prime }\left(k\right)(3u^2\left(k\right)v\left(k\right)-v^3\left(k\right))+a_3^{\prime }\left(k\right)v\left(k\right)+2a_2^{\prime }\left(k\right)u\left(k\right)v\left(k\right)\hfill \end{array}$$

Now at $k=k^{*},u\left(k^{*}\right)=0,v\left(k^{*}\right)=\sqrt{{\displaystyle \frac{a_3}{a_1}}}$

Hence we have,

$$\begin{array}{cc}\hfill D\left(k^{*}\right)& \hfill =a_3\left(k^{*}\right)-3a_1\left(k^{*}\right)v^2\left(k^{*}\right)=-2a_3\left(k^{*}\right)\hfill \\ \hfill E\left(k^{*}\right)& \hfill =2a_2\left(k^{*}\right)v\left(k^{*}\right)-4v^3\left(k^{*}\right)=2v\left(k^{*}\right)\left(a_2\left(k^{*}\right)-2{\displaystyle \frac{a_3\left(k^{*}\right)}{a_1\left(k^{*}\right)}}\right)\hfill \\ \hfill G\left(k^{*}\right)& \hfill =a_4^{\prime }\left(k^{*}\right)-a_2^{\prime }\left(k^{*}\right)v^2\left(k^{*}\right)=a_4^{\prime }\left(k^{*}\right)-a_2^{\prime }\left(k^{*}\right){\displaystyle \frac{a_3\left(k^{*}\right)}{a_1\left(k^{*}\right)}}\hfill \\ \hfill H\left(k^{*}\right)& \hfill =a_3^{\prime }\left(k^{*}\right)v\left(k^{*}\right)-v^3\left(k^{*}\right)a_1^{\prime }\left(k^{*}\right)=v\left(k^{*}\right)\left(a_3^{\prime }\left(k^{*}\right)-a_1^{\prime }\left(k^{*}\right){\displaystyle \frac{a_3\left(k^{*}\right)}{a_1\left(k^{*}\right)}}\right)\hfill \end{array}$$

On solving (8) and (9) using Cramer’s rule for $u^{\prime }\left(k^{*}\right)$ we get

$$\begin{array}{cc}\hfill {\left[{\displaystyle \frac{dRe\left(\lambda \right(k\left)\right)}{dk}}\right]}_{k=k^{*}}& \hfill =u^{\prime }\left(k^{*}\right)\hfill \\ \hfill & \hfill =-{\displaystyle \frac{G\left(k^{*}\right)D\left(k^{*}\right)+H\left(k^{*}\right)E\left(k^{*}\right)}{D^2\left(k^{*}\right)+E^2\left(k^{*}\right)}}\hfill \\ \hfill & \hfill =-{\displaystyle \frac{2a_3\left(k^{*}\right)}{a_1\left(k^{*}\right)}}{\displaystyle \frac{M}{D^2\left(k^{*}\right)+E^2\left(k^{*}\right)}}\ne 0\hfill \end{array}$$

provided $M\ne 0$ where

$M=a_1^2\left(k^{*}\right)(a_2{\left(k^{*}\right)}^{\prime }{a}_3\left(k^{*}\right)+a_2\left(k^{*}\right)a_3^{\prime }\left(k^{*}\right)-a_4^{\prime }\left(k^{*}\right)a_1\left(k^{*}\right))-2a_3\left(k^{*}\right)a_3^{\prime }\left(k^{*}\right)a_1\left(k^{*}\right)+a_1^{\prime }\left(k^{*}\right)(a_2\left(k^{*}\right)a_3\left(k^{*}\right)a_1\left(k^{*}\right)+2a_3^2\left(k^{*}\right))$

Hence the theorem follows. □

7. Numerical Simulation

We first search the equilibrium conditions at parameter values given in the Table 1. The corresponding interior equilibrium point $E^{*}=(16.1282,3.1825,0.9785,0.3282)$ with eigenvalues ($-0.0036\pm 0.2189i$, $-0.0754$, $-0.0132$). We have focused on the behavior of two types of infection within the prey population. First, we show the time series of the two infections under consideration in Figure 2 for the homogeneous distribution of the prey, which shows there is no impact of the herd here. However, we observe that due to the higher recovery and death rates of infection $I_1$ compared to $I_2$, the prevalence of $I_1$ decays sharply. Consequently, the other variables also stabilize, reaching a steady state. A similar behavior can be expected if the recovery and death rates of $I_2$ were higher than those of $I_1$.

(a) Effect of prey herd shape(k): At $k=1$, the prey distribution is homogeneous. The value of k reflects the herd shape of the prey, with higher values corresponding to more structured forms, such as cubes or spheres. For instance, $k=0.67$ represents a more defined shape like a sphere or cube, while $k=0.5$ depicts a 2D circular or square herd configuration. Figure 3 illustrates the occurrence of a Hopf bifurcation within this parameter space at $k=0.722467$, keeping all the other parameters as in Table 1. The first Lyapunov coefficient(FLC) at this point is $-1.914324\times {10}^{-4}$ which is negative thus confirming supercritical Hopf bifurcation. For the bifurcation analysis a numerical continuation bifurcation package MatCont has been used [52]. The periodic oscillations in the predator populations show that for lower values of k, the predator population decreases. As we increase the value of k, the population decay slows and eventually reaches stable limit cycles at the Hopf bifurcation point. However, we can observe a switching mechanism in the prey population with respect to strain 1. When $I_1$ is high, $I_2$ is low, and vice versa. Therefore, as $I_2$ decays during oscillations, $I_1$ increases and then decays, but eventually stabilizes at a higher value, while the $I_2$ infection dies out over time.

Table 1. Parameter Description.

Parameters	Numerical Values	Source
r	0.15	Assumed
K	25	Assumed
${\beta }_1$	0.01	[6]
${\beta }_2$	0.02	[6]
${\gamma }_1$	0.03	[6]
${\gamma }_2$	0.01	Assumed
a	0.5	[7]
k	0.55	[7]
h	2	[7]
${\alpha }_1$	0.15625	Assumed
${\alpha }_2$	0.8	Assumed
e	0.85	[7]
$e_1$	0.17	[7]
$e_2$	0.17	[7]
$\delta$	0.0145	[7]
$\mu$	0.5	[7]
${\mu }_1$	0.08	Assumed
${\mu }_2$	0.05	[7]

Table 1. Parameter Description.

Parameters	Numerical Values	Source
r	0.15	Assumed
K	25	Assumed
${\beta }_1$	0.01	[6]
${\beta }_2$	0.02	[6]
${\gamma }_1$	0.03	[6]
${\gamma }_2$	0.01	Assumed
a	0.5	[7]
k	0.55	[7]
h	2	[7]
${\alpha }_1$	0.15625	Assumed
${\alpha }_2$	0.8	Assumed
e	0.85	[7]
$e_1$	0.17	[7]
$e_2$	0.17	[7]
$\delta$	0.0145	[7]
$\mu$	0.5	[7]
${\mu }_1$	0.08	Assumed
${\mu }_2$	0.05	[7]

As we have observed the supercritical Hopf bifurcation, we attempt to identify the subcritical Hopf bifurcation in the $k-e_1$ parametric space. The occurrence of a subcritical Hopf bifurcation is particularly interesting due to the existence of an unstable limit cycle within the linearly stable region [15,16,17,18,19,20]. The demarcation point between the subcritical and supercritical Hopf bifurcations is the Generalized Hopf bifurcation point, which can be reached through the simultaneous variation of two parameters like $k-e_1$ here. In Figure 5, we have identified the Generalized Hopf bifurcation which is a codimension-2 bifurcation of fixed point, from which a limit point bifurcation of the limit cycle curve (LPC curve) emerges. In this bifurcation, the unstable limit cycle folds into a stable limit cycle of higher amplitude. This LPC curve is also referred to as the global stability boundary by Pandey et al. [15,16]

(b)Effect of death of predator due to prey herd($\delta$): In Figure 6, the parameter $\delta$ represents the rate of predator mortality due to the herd effect, and we also observe a Hopf bifurcation on this parameter at $\delta ={\delta }^{*}=0.0296$ where the First Lyapunov coefficient(FLC) is $-1.364586\times {10}^{-4}$ which is negative thus showing Supercritical Hopf bifurcation. At $\delta ={\delta }^{**}=0.0207$, we have a Hopf point with First Lyapunov Coefficient $2.085063\times {10}^{-4}$ which is positive thus confirming subcritical Hopf bifurcation. For an increased death rate of the predator, the oscillation amplitude of $I_1$ becomes very low because there are fewer predators to attack each population. As a result, due to the high recovery and death rates, $I_1$ diminishes to zero.

Figure 4. The time series plot for the Supercritical Hopf bifurcation showing the influence of the prey herd shape (k) on : (A) Susceptible prey,(B) prey infected with strain 1, (C) prey infected with strain 2, (D) Predator population.

Figure 4. The time series plot for the Supercritical Hopf bifurcation showing the influence of the prey herd shape (k) on : (A) Susceptible prey,(B) prey infected with strain 1, (C) prey infected with strain 2, (D) Predator population.

Figure 5. Hopf curve for $k-e_1$ parametric space where GH represents Generalized Hopf point.

Figure 5. Hopf curve for $k-e_1$ parametric space where GH represents Generalized Hopf point.

We demonstrate that the oscillation amplitude of $I_1$ is very low for higher values of $\delta$, as shown in Figure 6b. Furthermore, a switching mechanism is observed between $I_1$ and $I_2$. In Figure 7, for $\delta =0.018$, $I_1$ is higher while $I_2$ is dying, and for $\delta =0.028$, $I_1$ declines while $I_2$ becomes more prominent. However, oscillations are still present for higher values of $\delta$, though with a reduced amplitude due to the limited number of predators attacking each population. Figure 8 shows that only one infection survives, as $I_1$ dies, while the $I_2$ shows the growing oscillations for higher value of $\delta$.

Figure 6. The Supercritical Hopf-bifurcation diagram for $\delta$ shows the onset of periodic oscillations following the bifurcation point at $\delta ={\delta }^{*}=0.0296$. In the diagram, red and blue colors indicate the maximum and minimum values of the positive solution during the steady-state phase, respectively. For values of $\delta$ less than ${\delta }^{*}$, the convergence of the maximum and minimum values indicates the stability of $E^{*}$. Beyond this point, the solution begins to oscillate between these maximum and minimum values, indicating a loss of stability.

Figure 6. The Supercritical Hopf-bifurcation diagram for $\delta$ shows the onset of periodic oscillations following the bifurcation point at $\delta ={\delta }^{*}=0.0296$. In the diagram, red and blue colors indicate the maximum and minimum values of the positive solution during the steady-state phase, respectively. For values of $\delta$ less than ${\delta }^{*}$, the convergence of the maximum and minimum values indicates the stability of $E^{*}$. Beyond this point, the solution begins to oscillate between these maximum and minimum values, indicating a loss of stability.

Figure 7. The time series plot for the Subcritical hopf bifurcation showing the influence of the predator mortality due to prey herd ($\delta$) on: (A) Susceptible prey,(B) prey infected with strain 1, (C) prey infected with strain 2, (D) Predator population

Figure 7. The time series plot for the Subcritical hopf bifurcation showing the influence of the predator mortality due to prey herd ($\delta$) on: (A) Susceptible prey,(B) prey infected with strain 1, (C) prey infected with strain 2, (D) Predator population

Figure 8. The time series plot for the Supercritical hopf bifurcation showing the influence of the predator mortality due to prey herd ($\delta$) on: (A) Susceptible prey,(B) prey infected with strain 1, (C) prey infected with strain 2, (D) Predator population

Figure 8. The time series plot for the Supercritical hopf bifurcation showing the influence of the predator mortality due to prey herd ($\delta$) on: (A) Susceptible prey,(B) prey infected with strain 1, (C) prey infected with strain 2, (D) Predator population

8. Conclusions

In this article, we studied the dynamics of two strains of infection in the prey population. The disease transmission is not vertical,i.e. not transferable to offspring. The susceptible prey is assumed to form a herd against predators. In 2D, prey herd for a circle or square is $0.5$, and in 3D, for a sphere, it is $0.67$. We have generalized the rate of prey herd between 0 and 1. Also, we have assumed predator mortality due to prey herd. According to [21], prey on the circumference of the herd can injure predators. The presented model is well-posed and all the solutions are positive and bounded. We identified the biologically feasible equilibrium points and discussed the stability criteria for the interior equilibrium point where all populations coexist. Further, it was observed that the ecosystem cannot disappear since the trivial equilibrium point $E_1(0,0,0,0)$ is a saddle point i.e. inherently unstable which is a positive point for the given model. We also show the only healthy prey equilibrium point $E_2(K,0,0,0)$ is locally asymptotically stable under certain conditions and if the basic reproduction number $R_0<1$. If $R_0^1>1$ then prey infected with the first strain will dominate and if $R_0^2>1$ then prey infected with the second strain will dominate. Elena et al. [24,25] analyzed two different types of models with two strains of infection in the prey population. Their mathematical findings display no coexistence with both disease strains. Similarly, Roman et al. [26] considered two strain infections in predators, Bosica et al. [27] considered two strain infections in prey with mutualism but could not establish coexistence of equilibrium with both disease strains. However, our work shows the existence of an interior equilibrium point $E^{*}(S^{*},I_1^{*},I_2^{*},P^{*})$ where both the infections coexist. In Figure 3 we detected a supercritical hopf bifurcation about $E^{*}$ for prey herd shape $k=0.722467$. Gupta and Dubey [12] observed a resembling situation for $k=0.6621$ in the case of one infection strain. Moreover, in Figure 2 our research reveals that for a particular herd shape at $k=0.950245$ prey infected with strain 1 dies and other populations settle to a stable limit cycle. Further investigations demonstrate changes in population on increasing the mortality rate($\delta$) of predators due to prey herd. It is observed in Figure 6, that at $\delta =0.0296$ system undergoes supercritical Hopf bifurcation and a stable limit cycle bifurcates. In addition to this, it is seen in Figure 8, for high values of $\delta$ (say$=0.04$) prey infected with strain 1 dies out. Further, in the model we have detected significant bifurcation phenomena, consisting of the generalized Hopf bifurcation, the cusp point, and the limit point of cycles (LPC) curve. However, a detailed analysis and extensive investigation of these bifurcations, including their biological inferences and further understanding of the system’s dynamics, will be attended in future work. These factors hold significant prospects in understanding the complex interplay between spread of infection and population dynamics and will form the basis of later studies.