These authors contributed equally to this work.

[2]Kwadwo Antwi-Fordjour These authors contributed equally to this work.

These authors contributed equally to this work.

[1]Department of Mathematics, Iowa State University, Ames, 50011, IA, USA

2]Department of Mathematics and Computer Science, Samford University,
Birmingham, 35229, AL, USA

# Exploring unique dynamics in a predator-prey model with generalist predator
and group defence in prey

Vaibhava Srivastava<EMAIL_ADDRESS><EMAIL_ADDRESS>Rana D. Parshad
<EMAIL_ADDRESS>*  [

###### Abstract

In the current manuscript we consider a predator-prey model where the predator
is modeled as a generalist using a modified Leslie-Gower scheme, and the prey
exhibits group defence via a generalised response. We show that the model
could exhibit finite time blow-up, contrary to the current literature (Eur.
Phys. J. Plus 137, 28). We also propose a new concept via which the predator
population blows up in finite time while the prey population _quenches_ in
finite time. The blow-up and quenching times are proved to be one and the
same. Our analysis is complemented by numerical findings. This includes a
numerical description of the basin of attraction for large data blow-up
solutions, as well as several rich bifurcations leading to multiple limit
cycles, both in co-dimension one and two. Lastly, we posit a delayed version
of the model with globally existing solutions for any initial data.

###### keywords:

Predator-prey model, Finite time blowup, Quenching, Basin of attraction

##  1 Introduction

The current work aims to comment on the dynamics of the predator-prey model
introduced in [1]. This model takes the following form,

$\left\\{\begin{array}[]{l}\displaystyle\frac{dX}{dt}=RX\left(1-\displaystyle\frac{X(t-\tau)}{K}\right)-\dfrac{MXY}{X^{p}+C}\vspace{2ex}\\\
\displaystyle\frac{dY}{dt}=\left(D-\dfrac{E}{X+A}\right)Y^{2}\end{array}\right.$
(1)

This system of equations represents the dynamics of an interacting predator
and prey species, where $X(t)$ and $Y(t)$ are a prey and predator population,
respectively. The prey population obeys logistic growth, with a delay term
representing gestation effect [2]. The prey population’s functional response
is non-standard, taking the form $\frac{MX}{X^{p}+C}$, where $p\geq 1$. If
$p=1$, this is the classical Holling type II response. If $p=2$, it is the
Holling type IV response. If $1<p<2$, it is called the Cosner response [3].
The predator population is modeled according to the modified Leslie-Gower
scheme [4]. The other parameters used in the model have the following
interpretations: $R$ and $K$ are the intrinsic growth rate and carrying
capacity of the prey population, $M$ is the maximum predation rate, $C$
represents the protection provided to the prey by the environment, $D$ is the
reproduction rate of the generalist predator by sexual reproduction, $E$ is
the maximum rate of death of predator population, and $A$ represents the
residual loss of the predator population due to the scarcity of prey species
$X$.

In the event that there is no delay, or $\tau=0$, the model is a subsystem of
the well-known Upadhyay-Rai food chain model that models the interaction
between three trophic groups, a prey species $x$ is depredated on by a middle
predator species $y$ that in turn is depredated on by a top predator species
$z$, [4]. The top predator ($Y(t)$ here) is assumed to be a generalist. That
is, it can switch its food source in the absence of its favorite food $X$. The
growth of this top predator $Y$ is modeled in a nonstandard way via the
modified Leslie–Gower scheme. Here the $Y$ population grows due to sexual
reproduction - essentially modeled as a direct product between the males and
females in the population, $\frac{Y}{2}\times\frac{Y}{2}=\frac{1}{4}Y^{2}$,
which is controlled by a density-dependent, growth rate
$\left(D-\frac{E}{X+A}\right)$. This term is positive if the food source, $X$,
is above a certain threshold, and in this case, the top predator population
grows, and if below this threshold, it decays.

The Upadhyay-Rai model has a long-standing history in the mathematical ecology
literature. It is the first model to show chaotic dynamics when the top
predator is a generalist - and the second after the Hastings-Powell model [5],
to show chaos in a three-species food chain model. Boundedness for the model
was first shown in [6]. However, in 2015, the Upadhyay-Rai model was shown to
possess solutions that blow up in finite time for sufficiently large initial
data [7], and then subsequently even for small initial data [8]. Since then,
many works have appeared in the literature that has attempted to either dampen
the blow-up dynamic present in the basic model via ecological mechanisms such
as time delay, fear effects, or refuge [9, 10], while others have attempted to
show that blow-up does not occur altogether, via mechanisms, such as toxins,
Allee effects, and fractional derivatives [11, 12].

The following results concerning the properties of system (1) are proved in
[1].

######  Lemma 1 (Lemma 3.1 [1]).

The positive quadrant $\mathbb{R}^{2}$ is invariant for system (1).

######  Theorem 1 (Theorem 3.2 [1]).

The solution of system (1), originating in $\mathbb{R}^{2}_{+}$ are bounded,
provided the conditions

$\displaystyle\mu<\dfrac{DX^{*}}{(X^{*})^{p}+C}\quad\&\quad
D-\dfrac{E}{X^{*}+A}<\dfrac{D}{M}\Big{(}\dfrac{X^{*}}{Y^{*}}\Big{)}^{2}$ (2)

holds, where $(X^{*},Y^{*})$ is the interior equilibrium point and
$\mu=\min\\{M,E\\}.$

In the current manuscript, we show the following:

  * •

Solutions to (1) are unbounded. In fact, they can blow up in finite time for
sufficiently large data, even under the parametric restrictions of Theorem 3.2
[1]. This is shown via Theorem 3.

  * •

The positive quadrant $\mathbb{R}^{2}$ is not invariant for system (1). This
is shown via Corollary 1.

  * •

Introducing delay (motivated by gestation period) in the predator yields
global in-time existence. This is shown via Theorem 6.

  * •

A new concept of quenching in finite time is introduced for the prey species.
Subsequently, we proved that the predator’s blow-up time and the prey’s
quenching time are one and the same. For the non-delayed model, this is shown
via Theorem 4 and numerically validated through Remark 4 and Figure 3.
Moreover, for the delayed model, it is proposed via Theorem 7 and numerically
validated through Remark 5 and Figure 5.

  * •

Several numerical simulations are performed to corroborate all of our
analytical findings. See Figs. [1,2,3,4,5,6,7,8,9].

  * •

It is numerically validated that the linear feedback control proposed in [1]
does not prevent the blow-up for both delayed and non-delayed models. This is
numerically shown via Remarks 6 and 7, and Figs [6,8].

  * •

In [1], the asymptotic stability of the interior equilibrium for the delayed
linear feedback control model (under some parametric restrictions) is proposed
and proved. This is disproved by providing a counterexample. This is
numerically shown via Remark 7 and Figure 8.

  * •

The basin of attraction is thoroughly discussed, supported by several theorems
(refer to Theorems 8 and 9) and illustrated through Example 1.

  * •

In [1], the authors observed the presence of multiple concurrent limit cycles.
Expanding on their findings, we propose several conjectures (refer to
Conjectures 3, 4, and 5) and provide numerical validations (refer to Remark 8,
and Figures 10 and 11) regarding the occurrence of two concurrent limit
cycles. These validations are obtained by systematically varying either one or
two parameters.

  * •

A detailed discussion is carried out with applications of our results to
ecological scenarios. Furthermore, several open questions are proposed and
discussed.

##  2 Preliminaries

Given a system of ODE, depending on the non-linearity in the equations, one
might not expect a solution to always exist globally in time. In particular,
solutions of some ODE may blow up in finite time [13]. Recall,

######  Definition 1.

(Finite time blow-up for ODE) We say that a solution $u(t)$, of a given ODE,
with suitable initial conditions, blows up at a finite time if

$\underset{t\rightarrow
T^{\ast}<\infty}{\lim}\left|u\left(t\right)\right|=+\infty,$

where $T^{\ast}<\infty$ is the blow-up time.

We need the following alternative (See [14, 15, 16, 17, 18]).

######  Proposition 2.

The two-component system (1) admits a unique local in time, classical solution
$\left(X,Y\right)$ on an interval $[0,T_{\max}]$, and either
(i) The solution is bounded on $[0,T_{\max})$ and it is global ( i.e.
$T_{\max}=+\infty$).
(ii) Or

$\underset{t\nearrow T_{\max}}{\lim}\max\hskip
1.4457pt\left|X(t)\right|+\left|Y(t)\right|=+\infty,$ (3)

in this case, the solution is not global, and we say that it blows up in
finite time $T_{\max}$, or it ceases to exist, where $T_{\max}<\infty$ denotes
the eventual blowing-up time.

######  Definition 2.

(Finite time quenching for ODE) We say that a solution $u(t)$ of a given ODE,
with suitable initial conditions, quenches in finite time if

$\underset{t\rightarrow
T^{\ast\ast}<\infty}{\lim}\left|\frac{du\left(t\right)}{dt}\right|=+\infty,$

where $T^{\ast\ast}<\infty$ is the quenching time.

The idea of quenching and it’s several applications can be found in [19, 20].

Figure 1: Time series simulation of the non-delayed predator-prey dynamics
given by the mathematical model in (1). The simulation used specific parameter
values: $R=1$, $K=1$, $M=1.2$, $p=2$, $D=0.5$, $E=0.2$, $A=0.2$, and $C=0.3$,
with initial data $[X(0),Y(0)]=[78,30]$. The parameter choices ensure that the
model satisfies the requirements for a blow-up, as stated by Theorem 1.

Figure 2:  (A) Dynamics of non-delayed predator-prey dynamics given by the
mathematical model in (1) when certain initial conditions are applied. The
colored regions indicate two possible outcomes: the red region represents a
non-blowup case, while the blue region represents a blow-up case if we choose
initial conditions from that region. The model was simulated using specific
parameter values: $R=1$, $K=1$, $M=1.2$, $p=2$, $D=0.5$, $E=0.2$, $A=0.2$, and
$C=0.3$. (B) Shape-preserving interpolation of the boundary of the region of
attraction. The blow-up tolerance parameter is assumed as $10^{8}.$

######  Theorem 3.

Consider the model (1). Then $Y(t)$ blows up in finite time, that is

$\mathop{\lim}_{t\rightarrow T^{*}<\infty}|Y(t)|\rightarrow\infty,$ (4)

as long as the initial data $(X_{0},Y_{0})$ is sufficiently large.

###### Proof.

Consider (1), with positive initial conditions $(X_{0},Y_{0})$. By integrating
the equation for the predator in time, we obtain

$Y=\frac{1}{\dfrac{1}{Y_{0}}-\displaystyle\int_{0}^{t}\left(D-\frac{E}{X(s)+A}\right)ds}.$

Now consider the continuous function:

$\psi\left(t\right)=\frac{1}{Y_{0}}-\int_{0}^{t}\left(D-\frac{E}{X(s)+A}\right)ds.$
(5)

Our goal is to show $\psi(T)=0$, at some time $T>0$. Now we know from
Proposition 2 since solutions are classical and so continuous $\forall\hskip
3.6854ptT$, where $t<T_{max}$ and $T<T_{max}$, where $T_{max}$ is the maximal
interval of existence of the local solutions, that for a given $\delta>0$,
$\delta<T_{max}$, we can choose $X_{0}$ sufficiently large such that

$X+A>\dfrac{E}{D},\ \forall\hskip 3.61371ptt\in[0,\delta),$ (6)

and so

$\frac{E}{X+A}<D\ \forall\hskip 3.61371ptt\in[0,\delta).$ (7)

This implies

$\int_{0}^{t}\left(\frac{E}{X+A}\right)ds<\int_{0}^{t}Dds=Dt,\ \text{\ for all
}t\in[0,\delta).$

This entails, that for some $T^{*}<\delta$,

$\displaystyle\psi\left(T^{*}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{Y_{0}}-\int_{0}^{T}\left(D-\frac{E}{X(s)+A}\right)ds$
$\displaystyle=$
$\displaystyle\frac{1}{Y_{0}}+\int_{0}^{T}\left(\frac{E}{X(s)+A}\right)ds-DT$
$\displaystyle<$
$\displaystyle\frac{1}{Y_{0}}+\int_{0}^{\delta}\left(\frac{E}{X_{0}}\right)ds-D\delta$
$\displaystyle=$ $\displaystyle\frac{1}{Y_{0}}+\left(\frac{\delta
E}{X_{0}}\right)-D\delta\leq 0.$

The last inequality follows by choosing $(X_{0},Y_{0})$ appropriately large.
Thus $\psi\left(T^{*}\right)\leq 0$, since $\psi\left(0\right)=Y_{0}>0$, a
simple application of the intermediate value theorem tells us that there
exists a $T^{**}<T^{*}$, where $\psi\left(T^{**}\right)=0$, and so $Y$ has
blown-up at the finite time $T^{**}$, by a standard comparison argument [15].
∎

A corollary to this is,

######  Corollary 1.

The positive quadrant $\mathbb{R}^{2}$ is not invariant for system (1).

######  Remark 1.

We next provide a counter-example to the boundedness claim via Theorem 3.2
[1]. Pick $R=1,K=1,M=1.2,p=2,D=0.5,E=0.2,A=0.2$ and $C=0.3$, then $X^{*}=0.2$
and $Y^{*}=0.226667$. Let’s calculate the fractions involved in parametric
restrictions given by (2),

$0.2=\min\\{1.2,0.2\\}<\dfrac{0.5\times 0.2}{(0.2)^{2}+0.3}=0.294118,$

and

$0=0.5-\dfrac{0.2}{0.2+0.2}<\dfrac{0.5}{1.2}\Big{(}\dfrac{0.2}{0.226667}\Big{)}^{2}=0.324585$

This choice of parameter satisfies the constraints of Theorem 1, but from the
time-series plot Figure 1, the predator $Y$ blows up to infinity in finite
time $t=T^{*}=6.789603\times 10^{-2}.$

######  Remark 2.

For the parameters: $R=1,K=1,M=1.2,p=2,D=0.5,E=0.2,A=0.2$ and $C=0.3$, the
parameters satisfies the parametric restriction given in Theorem 8:

$0.3=C>\frac{\left(\frac{E}{D}-A\right)^{p}(ADp+AD+DKp-Ep-E)}{E-AD}=0.28.$

The interior equilibrium $(X^{*},Y^{*})=(0.2,0.226667)$ is locally asymptotic
stable.

######  Remark 3.

Note, via a simple comparison to the problem $\frac{dY}{dt}=DY^{2}$, the blow-
up time $T^{**}$ for (1), has a lower estimate given by,
$\frac{1}{DY_{0}}<T^{*}$. Also via standard theory, $T_{max}<T^{**}$. Here
$T_{max}$ was the maximal interval of existence assumed on $Y(t)$, so that we
could make the formal estimates.

Figure 3: Time series simulation of the predator-prey dynamics given by the
mathematical model in (1). The simulation used specific parameter values:
$R=1$, $K=1$, $M=1.2$, $p=2$, $D=0.5$, $E=0.2$, $A=0.2$, and $C=0.3$, with
initial data $[X(0),Y(0)]=[78,30]$.

######  Remark 4.

Let’s numerically validate that predator-prey dynamics given by the
mathematical model in (1) quenches in finite time. Consider the parameters as
$R=1,K=1,M=1.2,p=2,D=0.5,E=0.2,A=0.2$ and $C=0.3$. From the simulation, the
blow up time $T^{*}=6.789603\times 10^{-2}$, and $X(T^{*})=2.2434$ and
$Y(T^{*})=1.328\times 10^{16}$. Hence,

$\displaystyle\frac{dX}{dt}=-6.7078\times
10^{15}\quad\&\quad\displaystyle\frac{dY}{dt}=7.3834\times 10^{31}.$

This motivates the following Theorem:

######  Theorem 4.

Consider the system (1), for any choice of parameters and initial conditions
chosen sufficiently large, the prey population $X(t)$ will quench at a finite
time $T^{*}$, which is also a blow-up time for predator population $Y(t)$,
i.e.,

$\underset{t\rightarrow
T^{\ast}<\infty}{\lim}\left|\frac{dX\left(t\right)}{dt}\right|=+\infty\quad\&\quad\mathop{\lim}_{t\rightarrow
T^{*}<\infty}|Y(t)|\rightarrow\infty,$

This is proved via a series of two lemmas.

######  Lemma 2.

Consider the system (1), for any choice of parameters and initial conditions
chosen sufficiently large, the prey population $X(t)$ will quench at a finite
time $T^{**}$, which is greater than or equal to the blow-up time $T^{*}$ for
predator population $Y(t)$.

###### Proof.

From the Theorem 3, we have the $Y(t)$ blows up in finite time for
sufficiently large initial data $(X_{0},Y_{0})$, that is

$\mathop{\lim}_{t\rightarrow T^{*}<\infty}|Y(t)|\rightarrow\infty.$

We now claim that $T^{\ast\ast}$ the finite-time quenching for prey population
is greater than or equal to the blow-up time $T^{*}$ for predator population
$Y(t)$, i.e., $T^{\ast\ast}\geq T^{\ast}$. To claim this it suffices that
$\underset{t\rightarrow T^{\ast}<\infty}{\lim}X(t)\neq 0.$ Let’s on contrary
assume that

$\mathop{\lim}_{t\rightarrow T^{*}<\infty}X(t)=0\implies\text{From Proposition
}\ref{prop.exis.blow.2.2},\text{ we have}\quad X(t)\equiv 0,\quad\forall\hskip
3.61371pt\Omega_{\epsilon}=(T^{\ast}-\epsilon,T^{\ast}).$

Consider, the $Y$ nullclines in $\Omega_{\epsilon}$, we have

$\frac{dY}{dt}=\left(D-\dfrac{E}{A}\right)Y^{2}<0.$

Therefore, $Y$ is decreasing in $\Omega_{\epsilon}$, which contradicts the
fact that $Y(t)$ blows up in finite time.

Hence,

$\mathop{\lim}_{t\rightarrow T^{*}<\infty}X(t)=\widehat{x}(\neq
0)\quad\&\quad\mathop{\lim}_{t\rightarrow T^{*}<\infty}Y(t)=+\infty.$

Therefore, we get

$\underset{t\rightarrow
T^{\ast}<\infty}{\lim}\left|\frac{dX\left(t\right)}{dt}\right|=\underset{t\rightarrow
T^{\ast}<\infty}{\lim}\Big{\\{}RX\left(1-\displaystyle\frac{X(t)}{K}\right)-\dfrac{MXY}{X^{p}+C}\Big{\\}}=+\infty.$

Hence, we have that $T^{\ast\ast}$ the finite-time quenching for prey
population is greater than or equal to the blow-up time $T^{*}$ for predator
population $Y(t)$. ∎

Next we show,

######  Lemma 3.

Consider the system (1), for any choice of parameters and initial conditions
chosen sufficiently large, the prey population $X(t)$ will quench at a finite
time $T^{**}$, which is less than or equal to the blow-up time $T^{*}$ for
predator population $Y(t)$.

###### Proof.

We prove this via contradiction. Assume to the contrary that the quenching
time of the prey, $T^{**}>T^{*}$, where $T^{*}$ is the blow-up time of the
predator. Since $T^{\ast\ast}$ is the finite quenching time, we have
$\underset{t\rightarrow
T^{\ast\ast}<\infty}{\lim}\left|\frac{dX\left(t\right)}{dt}\right|=+\infty$.
From Proposition 2, and using the fact that the prey population is bounded
from above, we have

$\underset{t\rightarrow
T^{\ast\ast}<\infty}{\lim}RX\left(1-\displaystyle\frac{X(t)}{K}\right)<\infty.$

Upon using this information in the prey equation equation, and passing the
limit, we get

$\underset{t\rightarrow
T^{\ast\ast}<\infty}{\lim}\dfrac{XY}{X^{p}+C}=+\infty\quad\implies\underset{t\rightarrow
T^{\ast\ast}<\infty}{\lim}\quad Y(t)=+\infty.$

This gives us a contradiction, because if $T^{**}>T^{*}$, where we assume that
$T^{\ast}$ is the finite time of blowup, then the solution will cease to exist
past $T^{\ast}$.

This shows that the prey population $X(t)$ will quench at a finite time
$T^{**}$, which is less than or equal to the blow-up time $T^{*}$ for predator
population $Y(t)$. ∎

The Lemmas 2 and 3, used in conjunction prove Theorem 4, giving the
equivalence of $T^{**}$ and $T^{*}$.

##  3 Delayed ODE

Adding the delay parameter $\tau$, to the system (1), the delayed model is
given by

$\left\\{\begin{array}[]{l}\displaystyle\frac{dX}{dt}=X\left(1-\displaystyle
X(t-\tau)\right)-\dfrac{MXY}{X^{p}+C}\vspace{2ex}\\\
\displaystyle\frac{dY}{dt}=\left(D-\dfrac{E}{X+A}\right)Y^{2}\end{array}\right.$
(9)

subject to initial conditions

$\displaystyle X(\theta)$ $\displaystyle=\phi_{1}(\theta)>0,$ $\displaystyle
Y(\theta)$ $\displaystyle=\phi_{2}(\theta)>0,$ $\displaystyle\theta$
$\displaystyle\in[-\tau,0);\quad\phi_{i}(0)>0,\quad i=1,2.$

######  Theorem 5.

Consider the model (9). Then

$Y(t)$ blows up in finite time, that is

$\mathop{\lim}_{t\rightarrow T^{*}<\infty}|Y(t)|\rightarrow\infty,$ (10)

as long as the initial data $(X_{0},Y_{0})$ is sufficiently large.

###### Proof.

The proof follows similarly as in the proof of Theorem 3.

∎

Figure 4: Time series simulation of the delayed predator-prey dynamics given
by the mathematical model in (9). The simulation used specific parameter
values: $R=1$, $K=1$, $M=1.2$, $p=2$, $D=0.5$, $E=0.2$, $A=0.2,\tau=1$, and
$C=0.3$, with initial data $[X(0),Y(0)]=[100,100]$. The parameter choices
ensure that the model satisfies the requirements for a blow-up, as stated by
Theorem 6.

We now derive a delayed version of (1) s.t. there is global in-time existence
of solutions for all initial conditions. We add the delay parameter $\tau$, to
the system (1), by considering a gestation effect in the predator,

$\left\\{\begin{array}[]{l}\displaystyle\frac{dX}{dt}=X\left(1-\displaystyle
X\right)-\dfrac{MXY}{X^{p}+C}\vspace{2ex}\\\
\displaystyle\frac{dY}{dt}=DY(t-\tau)^{2}-\dfrac{E}{X+A}Y^{2}\end{array}\right.$
(11)

subject to initial conditions

$\displaystyle X(\theta)$ $\displaystyle=\phi_{1}(\theta)>0,$ $\displaystyle
Y(\theta)$ $\displaystyle=\phi_{2}(\theta)>0,$ $\displaystyle\theta$
$\displaystyle\in[-\tau,0);\quad\phi_{i}(0)>0,\quad i=1,2.$

######  Theorem 6.

Consider the system (11), for any choice of a parameter such that
$\delta_{1}>0$, and any initial condition the predator population will remain
bounded for all time.

###### Proof.

We proceed by contradiction. Assume the predator population blows up at a
finite time $t=T^{*}$. Let us integrate the predator equation in the time
interval $[T^{*}-1,T^{*}]$. Thus we obtain

$\displaystyle
Y(T^{*})=D\int^{T^{*}}_{T^{*}-1}Y(t-\tau)^{2}dt-\int^{T^{*}}_{T^{*}-1}\dfrac{E}{X+A}Y^{2}dt+Y(T^{*}-1)$
$\displaystyle<D\int^{T^{*}}_{T^{*}-1}Y(t-\tau)^{2}+Y(T^{*}-1)$
$\displaystyle=D\int^{T^{*}-\tau}_{T^{*}-1-\tau}Y(t)^{2}dt+Y(T^{*}-1)$

this follows via the positivity of solutions and a simple change of time
variable. Now using standard estimates, we obtain,

$Y(T^{*})\leq D\left(\sup_{T^{*}-1-\tau\leq t\leq
T^{*}-\tau}Y(t)\right)[(T^{*}-\tau)-(T^{*}-\tau-1)]$ (13)

Via the assumption that blow-up occurs at $t=T^{*}$, we have that

$Y(T^{*})=\infty\leq D\left(\sup_{T^{*}-1-\tau\leq t\leq
T^{*}-\tau}Y(t)\right)<\infty$ (14)

Thus we have a contradiction. This proves the theorem.

∎

Figure 5: Time series simulation of the predator-prey dynamics given by the
mathematical model in (9). The simulation used specific parameter values:
$R=1$, $K=1$, $M=1.2$, $p=2$, $D=0.5$, $E=0.2$, $A=0.2,\tau=1$, and $C=0.3$,
with initial data $[X(0),Y(0)]=[100,100]$.

######  Remark 5.

Let’s numerically validate that predator-prey dynamics given by the
mathematical model in (1) quenches in finite time. Consider the parameters as:
$R=1$, $K=1$, $M=1.2$, $p=2$, $D=0.5$, $E=0.2$, $A=0.2,\tau=1$, and $C=0.3$.
From the simulation, the blow up time $T^{*}=2.029316\times 10^{-2}$, and
$X(T^{*})=5.4949$ and $Y(T^{*})=7.086\times 10^{15}$. Hence,

$\displaystyle\frac{dX}{dt}=-1.532\times
10^{15}\quad\&\quad\displaystyle\frac{dY}{dt}=2.3347\times 10^{31}.$

This motivated us to present a Theorem:

######  Theorem 7.

Consider the system (9), for any choice of a parameter such that the initial
condition is sufficiently large, then the prey population $X(t)$ will quench
in finite time $T^{*}$, which is also a blow-up time for predator population
$Y(t)$, i.e.,

$\underset{t\rightarrow
T^{\ast}<\infty}{\lim}\left|\frac{dX\left(t\right)}{dt}\right|=+\infty\quad\&\quad\mathop{\lim}_{t\rightarrow
T^{*}<\infty}|Y(t)|\rightarrow\infty,$

###### Proof.

The proof follows similarly as in the proof of Theorem 4.

∎

##  4 Feedback Control

This section shows that the linear feedback model proposed in [1] to control
the system instability is not enough to prevent blow-up for both the delayed
and non-delayed models. Consider the linear feedback model proposed in [1].

$\left\\{\begin{array}[]{l}\displaystyle\frac{dX}{dt}=X\left(1-X(t-\tau)\right)-\dfrac{MXY}{X^{p}+C}-u(X-X^{*})\vspace{2ex}\\\
\displaystyle\frac{dY}{dt}=\left(D-\dfrac{E}{X+A}\right)Y^{2}-u(Y-Y^{*})\end{array}\right.$
(15)

where $(X^{*},Y^{*})$ are the interior equilibrium point of system (15) when
linear feedback control parameter $u$ is set to be zero.

The authors in [1] justified that the interior equilibrium for non-delayed
(15) is stable if the parameters satisfy

$\displaystyle u>\dfrac{1}{2}\Big{(}A_{11}-X^{*}\Big{)}\quad\&\quad
u^{2}-(A_{11}-X^{*})u+A_{12}A_{21}>0,$ (16)

where

$\displaystyle A_{11}=\dfrac{MpY^{*}(X^{*})^{p}}{((X^{*})^{p}+C)^{2}},\quad
A_{12}=\dfrac{mX^{*}}{(X^{*})^{p}+C},\quad\&\quad
A_{21}=\dfrac{D^{2}(Y^{*})^{2}}{E}.$ (17)

######  Remark 6.

Let’s consider a few parameters that satisfy the parametric restriction given
by (16), i.e.,

$M=1.2,\quad p=2,\quad C=0.3,\quad D=0.4,\quad E=0.2,\quad A=0.2\quad\&\quad
u=0.02.$

Figure 6 shows that the solution blows up in a finite time.

Moreover, the authors in [1] also argued about asymptotic stability of
interior for delayed (15) is stable if the parameters for $\lambda=i\omega$
simultaneously satisfies

$\begin{array}[]{l}u^{2}-A_{11}u+A_{12}A_{21}+X^{*}u\cos(\omega_{0}\tau)-\omega_{0}^{2}+X^{*}\omega_{0}\sin(\omega_{0}\tau)=0,\\\
(2u-A_{11})\omega_{0}+X^{*}\omega_{0}\cos(\omega_{0}\tau)-X^{*}u\sin(\omega_{0}\tau)>0,\end{array}$
(18)

where $A_{11},A_{12},A_{21},$ and $A_{22}$ are given by (17).

######  Remark 7.

Let’s consider a few parameters that satisfy the parametric restriction given
by (18), i.e.,

$M=1.2,\quad p=2,\quad C=0.3,\quad D=0.4,\quad E=0.2,\quad
A=0.2,\quad\tau=2\quad\&\quad u=0.06.$

From Figure 8, it is evident that the solution blows up in finite time.

Figure 6: Time series simulation of the non-delayed predator-prey dynamics
given by the mathematical model in (15). The simulation used specific
parameter values: $M=1.2$, $p=2$, $D=0.4$, $E=0.2$, $A=0.2,u=.02,\tau=0$ and
$C=0.3$, with initial data $[X(0),Y(0)]=[100,100]$.

Figure 7: (A) Dynamics of non-delayed predator-prey dynamics with feedback
control given by the mathematical model in (15) when certain initial
conditions are applied. The colored regions indicate two possible outcomes:
the red region represents a non-blowup case, while the blue region represents
a blow-up case if we choose initial conditions from that region. The
simulation used specific parameter values: $M=1.2$, $p=2$, $D=0.4$, $E=0.2$,
$A=0.2,u=.02,\tau=0$ and $C=0.3$ (B) Shape-preserving interpolation of the
boundary of the region of attraction. The blow-up tolerance parameter is
assumed as $10^{8}.$

Figure 8: Time series simulation of the delayed predator-prey dynamics given
by the mathematical model in (15). The simulation used specific parameter
values: $M=1.2$, $p=2$, $D=0.4$, $E=0.2$, $A=0.2,u=.06,\tau=2$ and $C=0.3$,,
with initial data $[X(0),Y(0)]=[200,200]$.

##  5 The blow-up boundary

In this section we ascertain the curve in the phase space that delineates
initial conditions that lead to blow-up in finite time versus those that are
attracted to equilibrium or limit cycles, hence stay bounded - this is the
_blow-up boundary_.

WLOG assumes $R=K=1$. To prove that system (1) and (9) blow up in finite time,
we will use the blow-up techniques mentioned in [21, 22].

######  Conjecture 1.

Consider the system (1), for any choice of parameter such that $\delta_{1}>0$,
and the initial condition satisfies the largeness conditions

$\displaystyle\ln\Big{(}\dfrac{|X(0)|}{\frac{E}{D-\delta_{1}}-A}\Big{)}>\frac{1}{\delta_{1}|Y_{0}|},$
(19)

then the predator population will explode to infinity in finite time, i.e.,

$\lim_{t\to T^{*}<\infty}||Y||\to\infty,$

where the blow-up time is given by

$T^{*}\leq\frac{1}{\delta_{1}|Y_{0}|}$

Let’s analyze this conjecture by considering a few cases:

Case I: Assume $D>\frac{E}{A}$. Hence the result follows trivially, since then
the predator nullcline would be of the form

$\displaystyle\frac{dY}{dt}\geq\widehat{C}Y^{2},$

where $\widehat{C}>0$, and hence the predator blows up to infinity in finite
time.

Case II Assume $\underset{X(t)}{\min}\\{D-\frac{E}{X+A}\\}=\delta_{1}>0$, then
again we have the finite time blow-up for predator population since

$\displaystyle\frac{dY}{dt}\geq\delta_{1}Y^{2}.$

Case II: Let’s work on non-obivous case, i.e., when
$D-\frac{E}{X+A}>\delta_{1}>0$. For this case, we need to make sure the
$D-\frac{E}{X+A}$ has the same sign, i.e.,

$\displaystyle D-\frac{E}{X+A}>\delta_{1}>0\quad\forall\hskip 3.61371ptX$

Equivalently,

$\displaystyle X>\frac{E}{D-\delta_{1}}-A>0$ (20)

Let’s consider the prey nullcline, and using the positivity and logistic
comparison arguments, we have $|X(t)|\leq 1\hskip 3.61371pt\forall\hskip
3.61371ptt$. On using these estimates in prey nullclines, we have

$\displaystyle\frac{dX}{dt}=X\left(1-\displaystyle
X\right)-\dfrac{MXY}{X^{p}+C}\quad\implies\quad X\geq X(0)e^{-t}.$

Using the lower bound for prey population $X$ in (20), we have

$\displaystyle|X(0)|e^{-t}>\frac{E}{D-\delta_{1}}-A\quad\implies\quad\ln\Big{(}\dfrac{|X(0)|}{\frac{E}{D-\delta_{1}}-A}\Big{)}>t$

Moreover, we know that:

$\displaystyle\frac{dY}{dt}=\delta_{1}Y^{2},$

blows up in $t=T^{*}\leq\frac{1}{\delta_{1}|Y_{0}|}.$ Hence, we have

$\displaystyle\ln\Big{(}\dfrac{|X(0)|}{\frac{E}{D-\delta_{1}}-A}\Big{)}>\frac{1}{\delta_{1}|Y_{0}|}\quad\implies\quad|Y_{0}|\ln\Big{(}\dfrac{|X(0)|}{\frac{E}{D-\delta_{1}}-A}\Big{)}>\dfrac{1}{\delta_{1}}.$

Therefore under the largeness restriction on the initial data given by (19),
the predator population will blow up to infinity in finite time.

######  Conjecture 2.

Consider the system (9), for any choice of a parameter such that
$\delta_{1}>0$, and the initial condition satisfies the largeness conditions
given by (19), then the predator population will explode to infinity in finite
time.

Consider the prey nullcline defined for the delayed model (9), and make use of
the logistics comparison argument, which entails that $X(t-\tau)\leq 1\hskip
3.61371pt\forall\hskip 3.61371ptt$. Hence, we have

$\displaystyle\displaystyle\frac{dX}{dt}$
$\displaystyle=X\left(1-\displaystyle
X(t-\tau)\right)-\dfrac{MXY}{X^{p}+C}\geq-X(t)\quad$
$\displaystyle\implies\quad X$ $\displaystyle\geq X(0)e^{-t}\quad\forall\hskip
3.61371ptt$

The lower-bound on prey population gives the following inequality,

$\displaystyle-\dfrac{1}{X+A}$ $\displaystyle\geq-\dfrac{1}{X(0)e^{-t}+A}$

On substituting this in predator nullclines, we have

$\displaystyle\displaystyle\frac{dY}{dt}$
$\displaystyle\geq\Big{(}D-\dfrac{1}{X(0)e^{-t}+A}\Big{)}Y^{2}$

So, if we have

$D-\dfrac{1}{X(0)e^{-t}+A}>\delta_{1}>0\quad\implies\Big{(}\dfrac{|X(0)|}{\frac{E}{D-\delta_{1}}-A}\Big{)}>t.$

Moreover, we know that:

$\displaystyle\frac{dY}{dt}=\delta_{1}Y^{2},$

blows up in $t=T^{*}\leq\frac{1}{\delta_{1}|Y_{0}|}.$

$\displaystyle\ln\Big{(}\dfrac{|X(0)|}{\frac{E}{D-\delta_{1}}-A}\Big{)}>\frac{1}{\delta_{1}|Y_{0}|}\quad\implies\quad|Y_{0}|\ln\Big{(}\dfrac{|X(0)|}{\frac{E}{D-\delta_{1}}-A}\Big{)}>\dfrac{1}{\delta_{1}}.$

Given the condition on the initial data as stated in (19), the predator
population will experience blow-up, meaning it will grow to infinity within a
finite time.

Figure 9: Approximation of blow-up boundary with the help of nonlinear
function $\phi_{1}(x)$ for non-delayed predator-prey dynamics with and without
feedback control. The red dots in the figure represent points on the blow-up
boundary. In panel (A), the parameters for the non-delayed model without
feedback control are the same as those used in Figure 2. Panel (B) parameters
are identical to those in Figure 7.

##  6 The basin of attraction

The basin of attraction for the model under consideration is fairly large.
That is initial data much larger than carrying capacity can be attracted to
equilibrium. However, despite the non-delayed model being a two-species one
(thus chaos is precluded) the structure of the basin of attraction can be
richer than a unique interior equilibrium or limit cycle dynamics. We provide
a few details next.

We first provide a summary for the stability of the interior equilibrium point
in the following theorem.

######  Theorem 8.

The interior equilibrium of model (1) is locally asymptotically if
$C>\frac{\left(\frac{E}{D}-A\right)^{p}(ADp+AD+DKp-Ep-E)}{E-AD}$.

###### Proof.

The proof follows directly from the linear stability analysis in the appendix.
∎

Next, we discuss the critical point of the prey’s protection by the
environment ($C$) around the interior equilibrium point where changes in
stability occur.

######  Theorem 9.

The model (1) experiences a Hopf bifucation around the interior equilibrium
point $E_{2}$ when $C$ crosses some critical value $C_{H}$ where
$C=C_{H}=\frac{\left(\frac{E}{D}-A\right)^{p}(ADp+AD+DKp-Ep-E)}{E-AD}$.

###### Proof.

The proof is routine and thus omitted for brevity. ∎

We shall provide an example to corroborate with the explicit parameter
dependent critical value $C_{H}$ from Theorem 9.

######  Example 1.

Consider the parameters $M=1.2,p=2,D=0.4,E=0.2$, and $A=0.2$ as used in [1].
Then

$\displaystyle C_{H}$
$\displaystyle=\frac{\left(\frac{E}{D}-A\right)^{p}(ADp+AD+DKp-Ep-E)}{E-AD}$
$\displaystyle=\dfrac{(0.5-0.2)^{2}(0.16+0.08+0.8-0.4-0.2)}{0.2-0.08}$
$\displaystyle=\dfrac{0.0393}{0.12}$ $\displaystyle=0.33$

The value obtained above corroborates with the $C_{H}$ value computed with a
computer algebra system as depicted in [1] (Figure 5(a) on page 13).

From this it appears if the interior equilibrium is a nodal sink its basin of
attraction can be fairly large. However if it is a spiral sink, this may not
be so. The interior equilibrium could loose stability via a Hopf bifurcation
as seen in Theorem 9 to yield a limit cycle. Thus fairly large initial
conditions could now be attracted to this limit cycle. The situation can be
more interesting and even richer dynamics could abound. We numerically see the
possibilities of multiple concurrent limit cycles. This was also observed in
[1]. Thus we conjecture that the occurrence of two concurrent limit cycles is
possible,

######  Conjecture 3.

Consider (1). Then there exist single parameter regimes, under which the
interior equilibrium point can bifurcate to produce at least two limit cycles.

This is also possible via varying two parameters.

######  Conjecture 4.

Consider (1). Then there exist two parameter regimes, under which a Bautin
bifurcation is possible.

Further these could collide and annihilate each other,

######  Conjecture 5.

Consider (1). Then there exist parameter regimes, under which a cyclical fold
bifurcation is possible.

Figure 10: The presence of two limit cycles is illustrated in the figures.
Left Figure: The stable interior equilibrium point $(0.2277,0.2264)$ is
denoted by a solid black circle, while the initial conditions $(0.274,0.228)$,
$(0.265,0.24)$, and $(0.20,0.192)$ are represented by solid red, black, and
blue squares, respectively. Right Figure: Period of cycles versus the
parameter $D$. Limit Point Cycle (LPC) observed at $D=0.4674$ and $D=0.4681$.
The chosen parameter values for this illustration are as follows: $R=1$,
$K=1$, $M=1.2$, $p=2$, $D=0.4676$, $E=0.2$, $A=0.2$, and $C=0.3$.

Figure 11: Two parameter bifurcation plots showing Bautin bifurcation. Left
Figure: $(D,A)=(0.3302,0.3777)$. Right Figure: $(D,C)=(0.4576,0.3055)$. Bottom
Figure: $(D,K)=(0.4699,0.9910)$. The chosen parameter values for this
illustration are as follows: $R=1$, $K=1$, $M=1.2$, $p=2$, $E=0.2$, $A=0.2$,
and $C=0.3$.

######  Remark 8.

In Figure 10, the initial conditions emanating from the solid blue and black
squares exhibit convergence towards an outer stable limit cycle, whereas the
initial conditions from black and red squares display movement away from an
inner unstable limit cycle. To visualize the limit point cycle (LPC) in an
alternate way, we show a curve with limit points clearly indicating the
existence of two limit cycles for $D<0.4681$. Additionally, the normal form
coefficient of the LPC derived with the aid of MATCONT is $-4.01\times
10^{-2}\neq 0$, thus the limit cycle has a fold.

In these situations the basin structure would be richer. The basin for the
interior (stable) equilibrium would be small, and larger data would be
attracted to the outer limit cycle. Thus bi-stable behavior, between a stable
equilibrium and a stable limit cycle is possible, with very large data leading
to blow-up solutions.

##  7 Numerical Simulations

We performed a time series simulation of the mathematical model given by (1)
and (9) using specific parameter values: $R=1$, $K=1$, $M=1.2$, $p=2$,
$D=0.5$, $E=0.2$, $A=0.2$, and $C=0.3$.

For the non-delay model (1), with initial data of $[X(0),Y(0)]=[78,30]$, we
observed that the predator population experiences blow-up to infinity in a
finite time of $T^{*}\approx 6.789603\times 10^{-2}$. We ensured that the
chosen parameter values satisfied the requirements for a blow-up, as
stipulated by Theorem 1. The numerical time series simulation was generated
using `MATLAB`’s standard solver for ordinary differential equations, ode45.

For the delayed model with the same parameter values as the non-delay model
and an additional delay parameter $\tau=1$, we found that the predator
population explodes to infinity in a finite time of $T^{*}\approx
2.029316\times 10^{-2}$ for constant positive initial data of $[100,100]$ (See
Figure (4)). We ensured that the selected parameter values met the
requirements for a blow-up, as prescribed by Theorem 5. The time series
simulation was produced using dde23, which is `MATLAB`’s solver for delay
differential equations (DDEs) with constant delays.

We conducted numerical investigations to analyze the region of attraction or
where the solution is bounded. Figure 2 displays the results, where colored
regions indicate two possible outcomes. The red region represents a non-blowup
case, whereas the blue region represents a blow-up case if we choose initial
conditions from respective regions. The simulations were conducted using the
non-delayed model (1) with specific parameter values: $R=1$, $K=1$, $M=1.2$,
$p=2$, $D=0.5$, $E=0.2$, $A=0.2$, and $C=0.3$. We performed the simulations
$9\times 10^{4}$ times, picking an initial data point from a mesh grid where
$0\leq X(0),Y(0)\leq 100$.

##  8 Conclusion

The current manuscript aims to provide commentary on the dynamics of the
predator-prey model introduced in [1]. The authors establish several stability
and boundedness results for both the delayed and non-delayed models under
certain parameter restrictions. However, our findings reveal that the
solutions to (1) are actually unbounded and can even exhibit finite-time blow-
up for sufficiently large data, contrary to the specified parameter
restrictions of Theorem 3.2 [1]. We demonstrate this through Theorem 3.
Consequently, the positive quadrant $\mathbb{R}^{2}$ is not invariant for
system (1), as shown in Corollary 1.

In addition, the authors introduce a delay in the predator equation, motivated
by the gestation period, which ensures global in-time existence, as
demonstrated in Theorem 6.

Furthermore, we introduce the concept of quenching in finite time for the prey
species. Subsequently, we propose that the predator’s blow-up time and the
prey’s quenching time are equivalent. For the non-delayed model, this is
presented in Theorem 4 and validated numerically in Remark 4 and Figure 3.
Similarly, for the delayed model, this is shown via Theorem 7 and numerically
validated in Remark 5 and Figure 5. To validate our analytical findings, we
conduct several numerical simulations, as depicted in Figs
[1,2,3,4,5,6,7,8,9]. To the best of our knowledge, this is the first instance
of _simultaneous_ blow-up and quenching reported in the predator-prey
literature. It is also one of the first instances where a smooth RHS can lead
to quenching. To this end the finite time blow-up dynamic of the predator is
exploited. The idea of using blow-up solutions to mimic the explosive increase
of invasive species has been considered in the literature [9]. This has been
seen to fit population patterns where the explosive increase of invasive
Burmese Pythons in the Florida Everglades, has resulted in the sharp decline
of mammalian prey [23]. To this end the simultaneous blow-up and quenching
models introduced in the current work may be more apt, as these would not only
catch the sharp population increase in predator populations, but the sharp
decrease in their prey populations as well. Furthermore, as quenching requires
blow-up in the time derivative (while the state variable stays bounded and
positive) such a model could depict the sharp decline of prey to very low (but
still positive) levels, such that they fall below almost any sensitive (and
endangered) levels. To this end fitting some of these data sets via the
current models would make for an interesting future endeavor.

Furthermore, we provide numerical evidence that the linear feedback control
proposed in [1] is inadequate in preventing blow-ups for both delayed and non-
delayed models. This is demonstrated numerically in Remarks 6 and 7, as well
as Figs [6,8]. In [1], the authors proposed and proved the asymptotic
stability of the interior solution for the delayed linear feedback control
model under specific parameter restrictions. However, we disprove their claim
by providing a counterexample, as illustrated numerically in Remark 7 and
Figure 8.

The discussion on the basin of attraction is comprehensive, incorporating the
support of various theorems (refer to Theorems 8 and 9). To provide a
practical understanding, we present an illustrative example (refer to Example
1). In [1], the presence of multiple concurrent limit cycles was initially
observed by the authors. Expanding upon their findings, we present several
conjectures (refer to Conjectures 3, 4, and 5). To validate these conjectures
numerically, we conducted simulations and obtained compelling results (refer
to Remark 8, and Figures 10 and 11). These simulations involved systematic
variations of either one or two parameters, providing evidence for the
occurrence of two concurrent limit cycles. Note, the occurrence of multiple
limit cycles for both Gauss type as well as Holling-Tanner type predator prey
models has been shown in the literature [24, 25]. However, these models all
possess globally existing bounded solutions, for all positive initial
conditions and across all parameter regimes. However, this is not the case in
the current class of models, where the top predator is modeled via a modified
Leslie-Gower scheme. In this setting we know via Theorem 3 that finite time
blow-up is possible. Thus the dynamics of the system are richer than the basic
Gauss type as well as Holling-Tanner type predator prey models. Herein one
could it would be pertinent to ask the maximum radius of the outer (stable)
limit cycle.

Ecologically speaking one could see bi-stability between the interior
equilibrium and the stable outer limit cycle, see Fig. 10, where the predator
and prey populations would go between stable levels to cyclical populations.
Note further one could see blow-up solutions if initial conditions were very
far from the outer cycle. To this end one could conjecture a the possible
existence of a third larger cycle encircling the outer stable one - and thus
possible changes to the structure of the blow-up boundary seen in Fig. 2.
Currently the blow-up boundary is approximated numerically, see Fig 9, using a
function of the form $\phi(x)=\frac{ax}{bx-c}$, which follows from estimates
in Theorem 3. Later analysis shows fits of the form
$\phi(x)=\frac{1}{b\log(cx)}$ might also be feasible. The exact sufficient
(and necessary) conditions to compute the blow-up boundary for these class of
model systems remains open.

Overall, one must use extreme caution in deriving global in time existence
results, in predator-prey systems when the top predator is modeled via the
modified Leslie-Gower scheme. It is more apt to derive restrictions on the
initial data such that solutions remain bounded. Also, one can look at the
various ecological mechanisms that could prevent blow-up solutions, for
initial conditions that would indeed lead to blow-up without these mechanisms.

## Conflict of interest

The authors declare there is no conflict of interest in this paper.

## References

  *   * Patra et al. [2022] Patra, R.R., Kundu, S., Maitra, S.: Effect of delay and control on a predator–prey ecosystem with generalist predator and group defence in the prey species. The European Physical Journal Plus 137(1), 28 (2022)
  * Kuang [1993] Kuang, Y.: Delay Differential Equations: with Applications in Population Dynamics. Academic press, (1993)
  * Cosner et al. [1999] Cosner, C., DeAngelis, D.L., Ault, J.S., Olson, D.B.: Effects of spatial grouping on the functional response of predators. Theoretical population biology 56(1), 65–75 (1999)
  * Upadhyay et al. [1998] Upadhyay, R.K., Iyengar, S., Rai, V.: Chaos: an ecological reality? International Journal of Bifurcation and Chaos 8(06), 1325–1333 (1998)
  * Hastings and Powell [1991] Hastings, A., Powell, T.: Chaos in a three-species food chain. Ecology 72(3), 896–903 (1991)
  * Aziz-Alaoui [2002] Aziz-Alaoui, M.: Study of a leslie–gower-type tritrophic population model. Chaos, Solitons & Fractals 14(8), 1275–1293 (2002)
  * Parshad et al. [2015] Parshad, R.D., Kumari, N., Kouachi, S.: A remark on “study of a leslie–gower-type tritrophic population model”[chaos, solitons and fractals 14 (2002) 1275–1293]. Chaos, Solitons & Fractals 71, 22–28 (2015)
  * Parshad et al. [2017] Parshad, R.D., Quansah, E., Beauregard, M.A., Kouachi, S.: On “small” data blow-up in a three species food chain model. Computers & Mathematics with Applications 73(4), 576–587 (2017)
  * Parshad et al. [2016] Parshad, R.D., Quansah, E., Black, K., Beauregard, M.: Biological control via “ecological” damping: an approach that attenuates non-target effects. Mathematical biosciences 273, 23–44 (2016)
  * Verma et al. [2021] Verma, H., Antwi-Fordjour, K., Hossain, M., Pal, N., Parshad, R.D., Mathur, P.: A “double” fear effect in a tri-trophic food chain model. The European Physical Journal Plus 136, 1–17 (2021)
  * Misra and Raveendra Babu [2016] Misra, O., Raveendra Babu, A.: Mathematical study of a leslie–gower-type tritrophic population model in a polluted environment. Modeling earth systems and environment 2, 1–11 (2016)
  * Debnath et al. [2020] Debnath, S., Ghosh, U., Sarkar, S.: Global dynamics of a tritrophic food chain model subject to the allee effects in the prey population with sexually reproductive generalized-type top predator. Computational and Mathematical Methods 2(2), 1079 (2020)
  * Quittner and Souplet [2007] Quittner, P., Souplet, P.: Superlinear Parabolic Problems: Blow-up, Global Existence and Steady States. Birkhäuser Basel, (2007)
  * Friedman [2008] Friedman, A.: Partial Differential Equations of Parabolic Type. Courier Dover Publications, (2008)
  * Henry [2006] Henry, D.: Geometric Theory of Semilinear Parabolic Equations vol. 840. Springer, (2006)
  * Pazy [2012] Pazy, A.: Semigroups of Linear Operators and Applications to Partial Differential Equations vol. 44. Springer, (2012)
  * Smoller [2012] Smoller, J.: Shock Waves and Reaction—diffusion Equations vol. 258. Springer, (2012)
  * Rothe [1984] Rothe, F.: Global-solutions of reaction-diffusion systems. Lecture notes in mathematics 1072, 1–214 (1984)
  * Levine [1983] Levine, H.A.: The quenching of solutions of linear parabolic and hyperbolic equations with nonlinear boundary conditions. SIAM journal on mathematical analysis 14(6), 1139–1153 (1983)
  * Beauregard and Selcuk [2022] Beauregard, M.A., Selcuk, B.: Quenching estimates for a non-newtonian filtration equation with singular boundary conditions. Tamkang Journal of Mathematics (2022)
  * Parshad et al. [2017] Parshad, R.D., Bhowmick, S., Quansah, E., Agrawal, R., Upadhyay, R.K.: Finite time blow-up in a delayed diffusive population model with competitive interference. International Journal of Nonlinear Sciences and Numerical Simulation 18(5), 435–450 (2017)
  * Parshad et al. [2016] Parshad, R.D., Bhowmick, S., Quansah, E., Basheer, A., Upadhyay, R.K.: Predator interference effects on biological control: The “paradox” of the generalist predator revisited. Communications in Nonlinear Science and Numerical Simulation 39, 169–184 (2016)
  * Dorcas et al. [2012] Dorcas, M.E., Willson, J.D., Reed, R.N., Snow, R.W., Rochford, M.R., Miller, M.A., Meshaka Jr, W.E., Andreadis, P.T., Mazzotti, F.J., Romagosa, C.M., et al.: Severe mammal declines coincide with proliferation of invasive burmese pythons in everglades national park. Proceedings of the National Academy of Sciences 109(7), 2418–2422 (2012)
  * Hofbauer and So [1989] Hofbauer, J., So, J.W.H.: Multiple limit cycles for predator-prey models (1989)
  * Sáez and González-Olivares [1999] Sáez, E., González-Olivares, E.: Dynamics of a predator-prey model. SIAM Journal on Applied Mathematics 59(5), 1867–1878 (1999)

##  Appendix A Stability Analysis of Non-Delayed Model

In this section, we shall investigate the existence of biologically feasible
equilibria of the model (1). Also, we will analyze the dynamics of the model
around these equilibria.

###  A.1 Existence of Equilibrium Points

From model (1), there are three biologically feasible equilibria, namely

  * (i)

The extinction (or population-free) equilibrium point is given by
${E_{0}}=(0,0)$.

  * (ii)

The predator-free (or axial) equilibrium point is given by ${E_{1}}=(K,0)$.

  * (iii)

The unique interior equilibrium (or coexistence) point is given by

${E_{2}}=(X^{*},Y^{*})=\left(\frac{E-AD}{D},-\frac{\frac{CR(E-AD)}{D}-KR\left(\frac{E-AD}{D}\right)^{p}+R\left(\frac{E-AD}{D}\right)^{p+1}-CKR}{KM}\right),$

where $X^{*}$ and $Y^{*}$ are the positive solutions of the prey and predator
nullclines:

$\displaystyle 0$
$\displaystyle=R\left(1-\displaystyle\frac{X}{K}\right)-\dfrac{MY}{X^{p}+C}\vspace{2ex}$
$\displaystyle 0$ $\displaystyle=\left(D-\dfrac{E}{X+A}\right).$

###  A.2 Local Stability Guidelines

In this subsection, we shall investigate the possible parametric restrictions
for the local stability of the model dynamics around each equilibrium point
with the aid of the Jacobian matrix. The Jacobian matrix is given by

$J=\left[\begin{array}[]{cc}\frac{MY\left((p-1)X^{p}-C\right)}{\left(C+X^{p}\right)^{2}}-\frac{2RX}{K}+R&-\frac{MX}{C+X^{p}}\\\
\frac{EY^{2}}{(A+X)^{2}}&2Y\left(D-\frac{E}{A+X}\right)\\\
\end{array}\right].$

Next, we will investigate the possible stability criteria around each
equilibrium point.

####  A.2.1 Extinction Equilibrium

The Jacobian matrix evaluated at the extinction equilibrium ($J_{E_{0}}$)
becomes

$J_{E_{0}}=\left[\begin{array}[]{cc}R&0\\\ 0&0\\\ \end{array}\right],$

where

$\det(J_{E_{0}})=0\quad\text{and}\quad\operatorname{tr}(J_{E_{0}})=R.$ (21)

Clearly, the extinction equilibrium point is a saddle based on equation (21).

####  A.2.2 Predator-free Equilibrium

The Jacobian matrix evaluated at the predator-free equilibrium ($J_{E_{1}}$)
becomes

$J_{E_{1}}=\left[\begin{array}[]{cc}-R&-\frac{MK}{C+K^{p}}\\\ 0&0\\\
\end{array}\right],$

where

$\det(J_{E_{1}})=0\quad\text{and}\quad\operatorname{tr}(J_{E_{1}})=-R.$ (22)

The stability of the predator-free equilibrium cannot be determined from the
information given in equation (22) via linearization. Here, the equilibrium
point is said to be degenerate. The stability of the predator-free equilibrium
can be determined by using the center manifold theory.

####  A.2.3 Interior Equilibrium

The Jacobian matrix evaluated at the interior equilibrium ($J_{E_{2}}$)
becomes

$J_{E_{2}}=\left[\begin{array}[]{cc}J^{*}_{11}&J^{*}_{12}\\\
J^{*}_{21}&J^{*}_{22}\\\ \end{array}\right],$

where

$\displaystyle J^{*}_{11}$
$\displaystyle=\frac{CR(AD-E)+R\left(\frac{E}{D}-A\right)^{p}(AD(p+1)+DKp-E(p+1))}{DK\left(\left(\frac{E}{D}-A\right)^{p}+C\right)}$
$\displaystyle J^{*}_{12}$
$\displaystyle=-\frac{M(E-AD)}{D\left(\left(\frac{E}{D}-A\right)^{p}+C\right)}<0$
$\displaystyle J^{*}_{21}$
$\displaystyle=\frac{R^{2}(E-D(A+K))^{2}\left(\left(\frac{E}{D}-A\right)^{p}+C\right)^{2}}{EK^{2}M^{2}}>0$
$\displaystyle J^{*}_{22}$ $\displaystyle=0.$

Now, we note that

$\det
J_{E_{2}}=-J^{*}_{12}J^{*}_{21}=\frac{R^{2}(E-AD)(AD+DK-E)^{2}\left(\left(\frac{E}{D}-A\right)^{p}+C\right)}{DEK^{2}M}>0$
(23)

and

$\operatorname{tr}J_{E_{2}}=J^{*}_{11}=\frac{CR(AD-E)+R\left(\frac{E}{D}-A\right)^{p}(AD(p+1)+DKp-E(p+1))}{DK\left(\left(\frac{E}{D}-A\right)^{p}+C\right)}.$
(24)

The determinant of the interior equilibrium is always positive from equation
(23) whiles the trace can assume either a positive or a negative number.